From f195ebbb1ebd14fbdaa0c003421b026d7079d582 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 20 Nov 2018 18:45:57 +0900 Subject: [PATCH 01/23] Launching jet uncertainty evaluator, and modifying the object selection code for adjusting this efficiently --- analysis/interface/topObjectSelection.h | 15 +++ analysis/src/topObjectSelection.cc | 20 ++- nanoAOD/plugins/JetUncertaintyEvaluator.cc | 142 +++++++++++++++++++++ nanoAOD/plugins/JetUncertaintyEvaluator.h | 67 ++++++++++ nanoAOD/prod/run2_2016MC_NANO.py | 5 +- nanoAOD/python/JetMetUncertainty_cff.py | 13 ++ nanoAOD/python/nano_cff.py | 12 +- 7 files changed, 263 insertions(+), 11 deletions(-) create mode 100644 nanoAOD/plugins/JetUncertaintyEvaluator.cc create mode 100644 nanoAOD/plugins/JetUncertaintyEvaluator.h create mode 100644 nanoAOD/python/JetMetUncertainty_cff.py diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index b19f2cf4..90a0c283 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -76,6 +76,21 @@ class topObjectSelection : public nanoBase virtual bool additionalConditionForJet(UInt_t nIdx) {return true;}; virtual bool additionalConditionForBJet(UInt_t nIdx) {return true;}; + // In uncertainty study we need to switch the kinematic variables of jets + // The following variables are for this switch + // In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on. + // In default, these are same as the original ones, but when a user wants to study systematic uncertainty + // so that he/she needs to switch them to the evaluated ones, + // just touching them in anlalyser class will be okay, and this is for it. + virtual void GetJetMassPt(UInt_t nIdx, + Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi) + { + fJetMass = Jet_mass[ nIdx ]; + fJetPt = Jet_pt[ nIdx ]; + fJetEta = Jet_eta[ nIdx ]; + fJetPhi = Jet_phi[ nIdx ]; + } + public: std::vector muonSelection(); std::vector elecSelection(); diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index db377532..61f6187a 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -122,11 +122,14 @@ vector topObjectSelection::jetSelection() { vector jets; //float Jet_SF_CSV[19] = {1.0,}; for (UInt_t i = 0; i < nJet; ++i){ - if (Jet_pt[i] < cut_JetPt) continue; - if (std::abs(Jet_eta[i]) > cut_JetEta) continue; + Float_t fJetMass, fJetPt, fJetEta, fJetPhi; + GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi); + + if (fJetPt < cut_JetPt) continue; + if (std::abs(fJetEta) > cut_JetEta) continue; if (Jet_jetId[i] < cut_JetID) continue; TLorentzVector mom; - mom.SetPtEtaPhiM(Jet_pt[i], Jet_eta[i], Jet_phi[i], Jet_mass[i]); + mom.SetPtEtaPhiM(fJetPt, fJetEta, fJetPhi, fJetMass); bool hasOverLap = false; for (auto lep : recoleps){ if (mom.TLorentzVector::DeltaR(lep) < cut_JetConeSizeOverlap) hasOverLap = true; @@ -141,7 +144,7 @@ vector topObjectSelection::jetSelection() { BTagEntry::JetFlavor JF = BTagEntry::FLAV_UDSG; if (abs(Jet_hadronFlavour[i]) == 5) JF = BTagEntry::FLAV_B; //else if (abs(Jet_hadronFlavour[i]) == 4) JF = BTagEntry::FLAV_C; - auto bjetSF = m_btagSF.eval_auto_bounds("central", JF , Jet_eta[i], Jet_pt[i], Jet_btagCSVV2[i]); + auto bjetSF = m_btagSF.eval_auto_bounds("central", JF , fJetEta, fJetPt, Jet_btagCSVV2[i]); b_btagweight *= bjetSF; for (UInt_t iu = 0; iu < 19; iu++) { // Jet_SF_CSV[iu] *= m_btagSF.getSF(jet, Jet_btagCSVV2[i], Jet_hadronFlavour[i], iu); @@ -156,12 +159,15 @@ vector topObjectSelection::jetSelection() { vector topObjectSelection::bjetSelection() { vector bjets; for (UInt_t i = 0; i < nJet; ++i ) { - if (Jet_pt[i] < cut_BJetPt) continue; - if (std::abs(Jet_eta[i]) > cut_BJetEta) continue; + Float_t fJetMass, fJetPt, fJetEta, fJetPhi; + GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi); + + if (fJetPt < cut_BJetPt) continue; + if (std::abs(fJetEta) > cut_BJetEta) continue; if (Jet_jetId[i] < cut_BJetID) continue; if (cut_BJetTypeBTag[i] < cut_BJetBTagCut) continue; TLorentzVector mom; - mom.SetPtEtaPhiM(Jet_pt[i], Jet_eta[i], Jet_phi[i], Jet_mass[i]); + mom.SetPtEtaPhiM(fJetPt, fJetEta, fJetPhi, fJetMass); bool hasOverLap = false; for (auto lep : recoleps) { if (mom.TLorentzVector::DeltaR(lep) < cut_BJetConeSizeOverlap) hasOverLap = true; diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.cc b/nanoAOD/plugins/JetUncertaintyEvaluator.cc new file mode 100644 index 00000000..eec18723 --- /dev/null +++ b/nanoAOD/plugins/JetUncertaintyEvaluator.cc @@ -0,0 +1,142 @@ +#include "JetMetUncertainty.h" + + +//#define debugMode +using namespace edm; +using namespace std; + + +JetMetUncertainty::JetMetUncertainty(const edm::ParameterSet &iConfig) : + src_(consumes>(iConfig.getParameter("src"))), + rhoToken_(consumes(iConfig.getParameter("rho"))), + payloadName_(iConfig.getParameter("payloadName")), + jetResFilePath_(edm::FileInPath(iConfig.getParameter("jetResFile")).fullPath()), + jetResSFFilePath_(edm::FileInPath(iConfig.getParameter("jetResSFFile")).fullPath()) +{ + produces(); +} + + +void JetMetUncertainty::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + + +void JetMetUncertainty::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) { + edm::Service rng; + rng_ = &rng->getEngine(lumi.index()); +} + + +void JetMetUncertainty::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + runOnMC_ = !iEvent.isRealData(); + + Handle jets; + iEvent.getByToken(src_, jets); + + // The following codes are from: + // https://github.com/vallot/CATTools/blob/cat80x/CatProducer/plugins/CATJetProducer.cc + edm::Handle rhoHandle; + iEvent.getByToken(rhoToken_, rhoHandle); + const double rho = *rhoHandle; + + if ( !payloadName_.empty() ) { + // temp measure - payloadName should be AK4PFchs, but PHYS14_25_V2 does not have uncertainty + edm::ESHandle JetCorParColl; + iSetup.get().get(payloadName_, JetCorParColl); + JetCorrectorParameters const & JetCorPar = ( *JetCorParColl )[ "Uncertainty" ]; + + jecUnc = new JetCorrectionUncertainty(JetCorPar); + } + + JME::JetResolution jetResObj; + JME::JetResolutionScaleFactor jetResSFObj; + + if ( runOnMC_ ) { + jetResObj = JME::JetResolution(jetResFilePath_); + jetResSFObj = JME::JetResolutionScaleFactor(jetResSFFilePath_); + } + + std::vector vecJetJESUp; + std::vector vecJetJESDn; + std::vector vecJetJERUp; + std::vector vecJetJERDn; + + for ( unsigned int i = 0 ; i < jets->size() ; i++ ) { + const auto & jet = ( *jets )[ i ]; + + double jetPt = jet.pt(); + double jetEta = jet.eta(); + + // Computing JEC uncertainty + float fJESUp = 1, fJESDn = 1; + + if ( !payloadName_.empty() ) { + jecUnc->setJetEta(jetEta); + jecUnc->setJetPt(jetPt); // here you must use the CORRECTED jet pt + fJESUp = 1 + jecUnc->getUncertainty(true); + + jecUnc->setJetEta(jetEta); + jecUnc->setJetPt(jetPt); // here you must use the CORRECTED jet pt + fJESDn = 1 - jecUnc->getUncertainty(false); + } + + // Computing JER uncertainty + float fJERUp = 1, fJERDn = 1; + + if ( runOnMC_ ) { + // adding genJet + auto genJet = jet.genJetFwdRef(); + + JME::JetParameters jetPars = {{JME::Binning::JetPt, jetPt}, + {JME::Binning::JetEta, jetEta}, + {JME::Binning::Rho, rho}}; + const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. + const double cJERUp = jetResSFObj.getScaleFactor(jetPars, Variation::UP); + const double cJERDn = jetResSFObj.getScaleFactor(jetPars, Variation::DOWN); + + // JER - apply scaling method if matched genJet is found, + // apply gaussian smearing method if unmatched + if ( genJet.isNonnull() && deltaR(genJet->p4(), jet.p4()) < 0.2 + && std::abs(genJet->pt() - jetPt) < jetRes * 3 * jetPt ) + { + const double genJetPt = genJet->pt(); + const double dPt = jetPt - genJetPt; + + fJERUp = std::max(0., (genJetPt + dPt * cJERUp) / jetPt); + fJERDn = std::max(0., (genJetPt + dPt * cJERDn) / jetPt); + } else { + const double smear = CLHEP::RandGaussQ::shoot(rng_); + + fJERUp = ( cJERUp <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJERUp * cJERUp - 1) ); + fJERDn = ( cJERDn <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJERDn * cJERDn - 1) ); + } + } + + vecJetJESUp.push_back(fJESUp); + vecJetJESDn.push_back(fJESDn); + vecJetJERUp.push_back(fJERUp); + vecJetJERDn.push_back(fJERDn); + } + + auto outJetCor = std::make_unique(jets->size(), "Jet", false, true); + //outJetCor->setDoc("Jet energy scale factor uncertainty"); + + outJetCor->addColumn("jes_up", vecJetJESUp, "Scale factor for JES up", + nanoaod::FlatTable::FloatColumn); + outJetCor->addColumn("jes_dn", vecJetJESDn, "Scale factor for JES down", + nanoaod::FlatTable::FloatColumn); + + outJetCor->addColumn("jer_up", vecJetJERUp, "Scale factor for JER up", + nanoaod::FlatTable::FloatColumn); + outJetCor->addColumn("jer_dn", vecJetJERDn, "Scale factor for JER down", + nanoaod::FlatTable::FloatColumn); + + iEvent.put(move(outJetCor)); +} + +DEFINE_FWK_MODULE(JetMetUncertainty); + + diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.h b/nanoAOD/plugins/JetUncertaintyEvaluator.h new file mode 100644 index 00000000..c1d4423a --- /dev/null +++ b/nanoAOD/plugins/JetUncertaintyEvaluator.h @@ -0,0 +1,67 @@ +#ifndef JetMetUncertainty_H +#define JetMetUncertainty_H + +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Utilities/interface/RandomNumberGenerator.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" + +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "DataFormats/Common/interface/ValueMap.h" + +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" +#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" +#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" + +#include "JetMETCorrections/Modules/interface/JetResolution.h" + +#include "CLHEP/Random/RandomEngine.h" +#include "CLHEP/Random/RandGaussQ.h" + +//#define debugMode + +using namespace edm; +using namespace std; + + +class JetMetUncertainty: public edm::stream::EDProducer<> { +public: + explicit JetMetUncertainty(const edm::ParameterSet &iConfig); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) override; + void produce( edm::Event&, const edm::EventSetup& ) override; + + edm::EDGetTokenT src_; + + edm::EDGetTokenT rhoToken_; + std::string payloadName_; + const std::string jetResFilePath_, jetResSFFilePath_; + + bool runOnMC_; + + JetCorrectionUncertainty *jecUnc; + CLHEP::HepRandomEngine* rng_; +}; + +#endif // JetMetUncertainty_H + + diff --git a/nanoAOD/prod/run2_2016MC_NANO.py b/nanoAOD/prod/run2_2016MC_NANO.py index bb5fd40c..d2c569fd 100644 --- a/nanoAOD/prod/run2_2016MC_NANO.py +++ b/nanoAOD/prod/run2_2016MC_NANO.py @@ -22,13 +22,14 @@ process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(-1) + input = cms.untracked.int32(100) ) # Input source process.source = cms.Source("PoolSource", # fileNames = cms.untracked.vstring('file:10D3267C-6FBD-E611-AED8-70106F48BC1E.root'), - fileNames = cms.untracked.vstring('/store/user/jlee/tsW_13TeV_PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v4_miniAOD/miniAOD_226.root'), +# fileNames = cms.untracked.vstring('/store/user/jlee/tsW_13TeV_PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v4_miniAOD/miniAOD_226.root'), + fileNames = cms.untracked.vstring('/store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/50000/0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'), secondaryFileNames = cms.untracked.vstring() ) diff --git a/nanoAOD/python/JetMetUncertainty_cff.py b/nanoAOD/python/JetMetUncertainty_cff.py new file mode 100644 index 00000000..cd72887e --- /dev/null +++ b/nanoAOD/python/JetMetUncertainty_cff.py @@ -0,0 +1,13 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * + +jetUncEvaluatorTable = cms.EDProducer("JetMetUncertainty", + src = cms.InputTag("linkedObjects","jets"), + rho = cms.InputTag('fixedGridRhoFastjetAll'), + payloadName = cms.string('AK4PFchs'), + jetResFile = cms.string("nano/nanoAOD/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"), + jetResSFFile = cms.string("nano/nanoAOD/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt"), +) + +jetUncEvaluatorTables = cms.Sequence(jetUncEvaluatorTable) + diff --git a/nanoAOD/python/nano_cff.py b/nanoAOD/python/nano_cff.py index 7532b698..e6d95a18 100644 --- a/nanoAOD/python/nano_cff.py +++ b/nanoAOD/python/nano_cff.py @@ -22,9 +22,17 @@ def customise(process, doHadron=True, fastSim=False): # DATA process.NANOAODoutput.fileName = fileName jecFile = 'Summer16_07Aug2017All_V10_DATA' - + + process.RandomNumberGeneratorService.jetUncEvaluatorTable = cms.PSet( + engineName = cms.untracked.string('TRandom3'), + initialSeed = cms.untracked.uint32(1), + ) + customiseMuons(process) + process.load('nano.nanoAOD.JetMetUncertainty_cff') + process.nanoAOD_step += process.jetUncEvaluatorTable + process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") if doHadron: process.load('nano.nanoAOD.jetMLID_cff') @@ -65,7 +73,7 @@ def customise(process, doHadron=True, fastSim=False): ) process.es_prefer_jec = cms.ESPrefer("PoolDBESSource","jec") print "JEC based on", process.jec.connect - + process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(1000) process.MessageLogger.cerr.FwkSummary.reportEvery = cms.untracked.int32(1000) From eef611d2cf53ee2ccf1b42dbc57cd0b67045f7f7 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 20 Nov 2018 18:51:23 +0900 Subject: [PATCH 02/23] Sorry, I kept modifications for test... rolled back them --- nanoAOD/prod/run2_2016MC_NANO.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/nanoAOD/prod/run2_2016MC_NANO.py b/nanoAOD/prod/run2_2016MC_NANO.py index d2c569fd..bb5fd40c 100644 --- a/nanoAOD/prod/run2_2016MC_NANO.py +++ b/nanoAOD/prod/run2_2016MC_NANO.py @@ -22,14 +22,13 @@ process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(100) + input = cms.untracked.int32(-1) ) # Input source process.source = cms.Source("PoolSource", # fileNames = cms.untracked.vstring('file:10D3267C-6FBD-E611-AED8-70106F48BC1E.root'), -# fileNames = cms.untracked.vstring('/store/user/jlee/tsW_13TeV_PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v4_miniAOD/miniAOD_226.root'), - fileNames = cms.untracked.vstring('/store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/50000/0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'), + fileNames = cms.untracked.vstring('/store/user/jlee/tsW_13TeV_PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v4_miniAOD/miniAOD_226.root'), secondaryFileNames = cms.untracked.vstring() ) From 43d20ebd0e1939d372685ae1b3331aeb18defd4b Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 20 Nov 2018 19:56:02 +0900 Subject: [PATCH 03/23] Fixed naming --- nanoAOD/plugins/JetUncertaintyEvaluator.cc | 12 ++++++------ nanoAOD/plugins/JetUncertaintyEvaluator.h | 10 +++++----- nanoAOD/python/JetMetUncertainty_cff.py | 2 +- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.cc b/nanoAOD/plugins/JetUncertaintyEvaluator.cc index eec18723..89c479c7 100644 --- a/nanoAOD/plugins/JetUncertaintyEvaluator.cc +++ b/nanoAOD/plugins/JetUncertaintyEvaluator.cc @@ -1,4 +1,4 @@ -#include "JetMetUncertainty.h" +#include "JetUncertaintyEvaluator.h" //#define debugMode @@ -6,7 +6,7 @@ using namespace edm; using namespace std; -JetMetUncertainty::JetMetUncertainty(const edm::ParameterSet &iConfig) : +JetUncertaintyEvaluator::JetUncertaintyEvaluator(const edm::ParameterSet &iConfig) : src_(consumes>(iConfig.getParameter("src"))), rhoToken_(consumes(iConfig.getParameter("rho"))), payloadName_(iConfig.getParameter("payloadName")), @@ -17,20 +17,20 @@ JetMetUncertainty::JetMetUncertainty(const edm::ParameterSet &iConfig) : } -void JetMetUncertainty::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { +void JetUncertaintyEvaluator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.setUnknown(); descriptions.addDefault(desc); } -void JetMetUncertainty::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) { +void JetUncertaintyEvaluator::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) { edm::Service rng; rng_ = &rng->getEngine(lumi.index()); } -void JetMetUncertainty::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { +void JetUncertaintyEvaluator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { runOnMC_ = !iEvent.isRealData(); Handle jets; @@ -137,6 +137,6 @@ void JetMetUncertainty::produce(edm::Event &iEvent, const edm::EventSetup &iSetu iEvent.put(move(outJetCor)); } -DEFINE_FWK_MODULE(JetMetUncertainty); +DEFINE_FWK_MODULE(JetUncertaintyEvaluator); diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.h b/nanoAOD/plugins/JetUncertaintyEvaluator.h index c1d4423a..bce99b5a 100644 --- a/nanoAOD/plugins/JetUncertaintyEvaluator.h +++ b/nanoAOD/plugins/JetUncertaintyEvaluator.h @@ -1,5 +1,5 @@ -#ifndef JetMetUncertainty_H -#define JetMetUncertainty_H +#ifndef JetUncertaintyEvaluator_H +#define JetUncertaintyEvaluator_H #include @@ -41,9 +41,9 @@ using namespace edm; using namespace std; -class JetMetUncertainty: public edm::stream::EDProducer<> { +class JetUncertaintyEvaluator: public edm::stream::EDProducer<> { public: - explicit JetMetUncertainty(const edm::ParameterSet &iConfig); + explicit JetUncertaintyEvaluator(const edm::ParameterSet &iConfig); static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); private: @@ -62,6 +62,6 @@ class JetMetUncertainty: public edm::stream::EDProducer<> { CLHEP::HepRandomEngine* rng_; }; -#endif // JetMetUncertainty_H +#endif // JetUncertaintyEvaluater_H diff --git a/nanoAOD/python/JetMetUncertainty_cff.py b/nanoAOD/python/JetMetUncertainty_cff.py index cd72887e..71e87cf3 100644 --- a/nanoAOD/python/JetMetUncertainty_cff.py +++ b/nanoAOD/python/JetMetUncertainty_cff.py @@ -1,7 +1,7 @@ import FWCore.ParameterSet.Config as cms from PhysicsTools.NanoAOD.common_cff import * -jetUncEvaluatorTable = cms.EDProducer("JetMetUncertainty", +jetUncEvaluatorTable = cms.EDProducer("JetUncertaintyEvaluator", src = cms.InputTag("linkedObjects","jets"), rho = cms.InputTag('fixedGridRhoFastjetAll'), payloadName = cms.string('AK4PFchs'), From 4aeabb7344658ef2830890e45d9eeb57149ba4c0 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 5 Dec 2018 13:56:45 +0900 Subject: [PATCH 04/23] Evaluation of JER/JES uncertainty on the fly --- analysis/BuildFile.xml | 2 + analysis/interface/topEventSelectionDL.h | 4 +- analysis/interface/topEventSelectionSL.h | 4 +- analysis/interface/topObjectSelection.h | 56 ++++++++--- analysis/src/topEventSelectionDL.cc | 4 +- analysis/src/topEventSelectionSL.cc | 4 +- analysis/src/topObjectSelection.cc | 123 ++++++++++++++++++++++- 7 files changed, 172 insertions(+), 25 deletions(-) diff --git a/analysis/BuildFile.xml b/analysis/BuildFile.xml index 024e91ec..0f59fad4 100644 --- a/analysis/BuildFile.xml +++ b/analysis/BuildFile.xml @@ -1,4 +1,6 @@ + + diff --git a/analysis/interface/topEventSelectionDL.h b/analysis/interface/topEventSelectionDL.h index f1314edf..a65cd94f 100644 --- a/analysis/interface/topEventSelectionDL.h +++ b/analysis/interface/topEventSelectionDL.h @@ -52,8 +52,8 @@ class topEventSelectionDL : public topObjectSelection virtual int EventSelection(); - topEventSelectionDL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t dl = true, Bool_t sle = false, Bool_t slm = false); - topEventSelectionDL(TTree *tree=0, Bool_t isMC=false, Bool_t dl=true, Bool_t sle=false, Bool_t slm=false) : topEventSelectionDL(tree, 0, 0, isMC, dl, sle, slm) {} + topEventSelectionDL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t dl = true, Bool_t sle = false, Bool_t slm = false, UInt_t unFlag = 0); + topEventSelectionDL(TTree *tree=0, Bool_t isMC=false, Bool_t dl=true, Bool_t sle=false, Bool_t slm=false, UInt_t unFlag = 0) : topEventSelectionDL(tree, 0, 0, isMC, dl, sle, slm, unFlag) {} ~topEventSelectionDL(); virtual void Loop() = 0; diff --git a/analysis/interface/topEventSelectionSL.h b/analysis/interface/topEventSelectionSL.h index a7c80bcb..618c5d31 100644 --- a/analysis/interface/topEventSelectionSL.h +++ b/analysis/interface/topEventSelectionSL.h @@ -51,8 +51,8 @@ class topEventSelectionSL : public topObjectSelection virtual Bool_t TrigForEl() {return HLT_Ele27_WPTight_Gsf;} virtual Bool_t TrigForMu() {return HLT_IsoTkMu24 || HLT_IsoMu24;} - topEventSelectionSL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t sle = false, Bool_t slm = false); - topEventSelectionSL(TTree *tree=0, Bool_t isMC=false, Bool_t sle=false, Bool_t slm=false) : topEventSelectionSL(tree, 0, 0, isMC, sle, slm) {} + topEventSelectionSL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t sle = false, Bool_t slm = false, UInt_t unFlag = 0); + topEventSelectionSL(TTree *tree=0, Bool_t isMC=false, Bool_t sle=false, Bool_t slm=false, UInt_t unFlag = 0) : topEventSelectionSL(tree, 0, 0, isMC, sle, slm, unFlag) {} ~topEventSelectionSL(); virtual void Loop() = 0; }; diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 9a7aac4d..c019f981 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -5,6 +5,12 @@ #include "nano/external/interface/TopTriggerSF.h" //#include "nano/external/interface/TTbarModeDefs.h" +#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" +#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" +#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" +#include "JetMETCorrections/Modules/interface/JetResolution.h" +#include + class topObjectSelection : public nanoBase { protected: @@ -14,7 +20,16 @@ class topObjectSelection : public nanoBase Float_t b_maxBDiscr_nonb; + UInt_t m_unFlag; + + JetCorrectionUncertainty *jecUnc; + + JME::JetResolution jetResObj; + JME::JetResolutionScaleFactor jetResSFObj; + TRandom3 *rndEngine; + public: + Float_t Jet_pt_jer_nom[35]; // YOU MUST SET UP ALL IN THE BELOW!!! // (SetCutValues() will force you to do it) Float_t cut_ElectronPt; @@ -73,8 +88,11 @@ class topObjectSelection : public nanoBase return Muon_isPFcand[ nIdx ] && ( Muon_globalMu[ nIdx ] || Muon_trackerMu[ nIdx ] ); }; virtual bool additionalConditionForGenJet(UInt_t nIdx) {return true;}; - virtual bool additionalConditionForJet(UInt_t nIdx) {return true;}; - virtual bool additionalConditionForBJet(UInt_t nIdx) {return true;}; + + virtual bool additionalConditionForJet(UInt_t nIdx, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, Float_t &fJetMass) + {return true;}; + virtual bool additionalConditionForBJet(UInt_t nIdx, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, Float_t &fJetMass) + {return true;}; // In uncertainty study we need to switch the kinematic variables of jets // The following variables are for this switch @@ -83,13 +101,9 @@ class topObjectSelection : public nanoBase // so that he/she needs to switch them to the evaluated ones, // just touching them in anlalyser class will be okay, and this is for it. virtual void GetJetMassPt(UInt_t nIdx, - Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi) - { - fJetMass = Jet_mass[ nIdx ]; - fJetPt = Jet_pt[ nIdx ]; - fJetEta = Jet_eta[ nIdx ]; - fJetPhi = Jet_phi[ nIdx ]; - } + Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi); + + Int_t GetMatchGenJet(UInt_t nIdxJet); public: std::vector muonSelection(); @@ -102,15 +116,33 @@ class topObjectSelection : public nanoBase std::vector genJetSelection(); - topObjectSelection(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false); - topObjectSelection(TTree *tree=0, Bool_t isMC=false) : topObjectSelection(tree, 0, 0, isMC) {} - ~topObjectSelection() {} + topObjectSelection(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, UInt_t unFlag = 0); + topObjectSelection(TTree *tree=0, Bool_t isMC=false, UInt_t unFlag = 0) : + topObjectSelection(tree, 0, 0, isMC, unFlag) {} + ~topObjectSelection() { + if ( jecUnc != NULL ) delete jecUnc; + if ( rndEngine != NULL ) delete rndEngine; + } // In this function you need to set all the cut conditions in the above // If you do not set this function up (so that you didn't set the cuts), the compiler will deny your code, // so you can be noticed that you forgot the setting up. // And you don't need to run this function indivisually; it will be run in the creator of this class. virtual int SetCutValues() = 0; + + enum { + OptBit_JER_Up = 0, + OptBit_JER_Dn, + OptBit_JES_Up, + OptBit_JES_Dn + }; + + enum { + OptFlag_JER_Up = ( 1 << OptBit_JER_Up ), + OptFlag_JER_Dn = ( 1 << OptBit_JER_Dn ), + OptFlag_JES_Up = ( 1 << OptBit_JES_Up ), + OptFlag_JES_Dn = ( 1 << OptBit_JES_Dn ) + }; }; #endif diff --git a/analysis/src/topEventSelectionDL.cc b/analysis/src/topEventSelectionDL.cc index 0864f128..9eea06ee 100644 --- a/analysis/src/topEventSelectionDL.cc +++ b/analysis/src/topEventSelectionDL.cc @@ -2,8 +2,8 @@ using std::vector; -topEventSelectionDL::topEventSelectionDL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t dl, Bool_t sle, Bool_t slm) : - topObjectSelection(tree, had, hadTruth, isMC), +topEventSelectionDL::topEventSelectionDL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t dl, Bool_t sle, Bool_t slm, UInt_t unFlag) : + topObjectSelection(tree, had, hadTruth, isMC, unFlag), m_isDL(dl), m_isSL_e(sle), m_isSL_m(slm) { SetCutValues(); diff --git a/analysis/src/topEventSelectionSL.cc b/analysis/src/topEventSelectionSL.cc index c78cd5bf..f7a1da6e 100644 --- a/analysis/src/topEventSelectionSL.cc +++ b/analysis/src/topEventSelectionSL.cc @@ -2,8 +2,8 @@ using std::vector; -topEventSelectionSL::topEventSelectionSL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t sle, Bool_t slm) : - topObjectSelection(tree, had, hadTruth, isMC), +topEventSelectionSL::topEventSelectionSL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t sle, Bool_t slm, UInt_t unFlag) : + topObjectSelection(tree, had, hadTruth, isMC, unFlag), h_nevents(0), h_genweights(0), h_cutFlow(0), diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index 8660b47e..42637c61 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -2,9 +2,32 @@ using std::vector; -topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) : - nanoBase(tree, had, hadTruth, isMC) -{} +topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, UInt_t unFlag) : + nanoBase(tree, had, hadTruth, isMC), m_unFlag(unFlag), jecUnc(NULL), rndEngine(NULL) +{ + if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn | OptFlag_JER_Up | OptFlag_JER_Dn ) ) != 0 ) { + std::string env = getenv("CMSSW_BASE"); + + std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jetunc/" + "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; + jetResSFObj = JME::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); + + std::string strPathJetResObj = env + "/src/nano/analysis/data/jetunc/" + "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; + jetResObj = JME::JetResolution(strPathJetResObj.c_str()); + + rndEngine = new TRandom3(12345); + + if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { + std::string strPathJecUnc = env + "/src/nano/analysis/data/jetunc/" + //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt"; + "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; + + JetCorrectorParameters JetCorPar(strPathJecUnc, "Total"); + jecUnc = new JetCorrectionUncertainty(JetCorPar); + } + } +} vector topObjectSelection::elecSelection() { vector elecs; @@ -139,7 +162,7 @@ vector topObjectSelection::jetSelection() { if (mom.TLorentzVector::DeltaR(lep) < cut_JetConeSizeOverlap) hasOverLap = true; } if (hasOverLap) continue; - if ( !additionalConditionForJet(i) ) continue; + if ( !additionalConditionForJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue; auto jet = TParticle(); jet.SetMomentum(mom); @@ -204,7 +227,7 @@ vector topObjectSelection::bjetSelection() { if (mom.TLorentzVector::DeltaR(lep) < cut_BJetConeSizeOverlap) hasOverLap = true; } if (hasOverLap) continue; - if ( !additionalConditionForBJet(i) ) continue; + if ( !additionalConditionForBJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue; auto bjet = TParticle(); bjet.SetMomentum(mom); @@ -213,3 +236,93 @@ vector topObjectSelection::bjetSelection() { } return bjets; } + + +// In uncertainty study we need to switch the kinematic variables of jets +// The following variables are for this switch +// In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on. +// In default, these are same as the original ones, but when a user wants to study systematic uncertainty +// so that he/she needs to switch them to the evaluated ones, +// just touching them in anlalyser class will be okay, and this is for it. +void topObjectSelection::GetJetMassPt(UInt_t nIdx, + Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi) +{ + Float_t fCorrFactor = 1.0; + + fJetMass = Jet_mass[ nIdx ]; + fJetPt = Jet_pt[ nIdx ]; + fJetEta = Jet_eta[ nIdx ]; + fJetPhi = Jet_phi[ nIdx ]; + + if ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { + // Evaluating the central part cJER of the factor + JME::JetParameters jetPars = {{JME::Binning::JetPt, fJetPt}, + {JME::Binning::JetEta, fJetEta}, + {JME::Binning::Rho, fixedGridRhoFastjetAll}}; + + // We need corresponding genJet + //Int_t nIdxGen = Jet_genJetIdx[ nIdx ]; + Int_t nIdxGen = GetMatchGenJet(nIdx); + const double genJetPt = GenJet_pt[ nIdxGen ]; + + const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. + const double cJER = jetResSFObj.getScaleFactor(jetPars, + ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn ) ) == 0 ? Variation::NOMINAL : + ( ( m_unFlag & OptFlag_JER_Up ) != 0 ? Variation::UP : Variation::DOWN ) )); + + // JER (nominal and up and down) - apply scaling method if matched genJet is found, + // apply gaussian smearing method if unmatched + /*if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool + std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) + fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt);*/ + if ( nIdxGen >= 0 ) { + fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt; + } else { + // Different from the postprocess tool, but it might not be important + const double smear = rndEngine->Gaus(0, 1); + fCorrFactor = ( cJER <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJER * cJER - 1) ); + } + + if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { // JES + // The evaluator needs pT and eta of the current jet + jecUnc->setJetPt(fCorrFactor * fJetPt); + jecUnc->setJetEta(fJetEta); + + if ( ( m_unFlag & OptFlag_JES_Up ) != 0 ) { + fCorrFactor *= 1 + jecUnc->getUncertainty(true); + } else { + fCorrFactor *= 1 - jecUnc->getUncertainty(true); + } + } + } + + fJetMass *= fCorrFactor; + fJetPt *= fCorrFactor; +} + + +Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet) { + UInt_t i; + + Float_t dRFound = 999; + UInt_t nIdxFound = 999; + + for ( i = 0 ; i < nGenJet ; i++ ) { + Float_t dDEta = Jet_eta[ nIdxJet ] - GenJet_eta[ i ]; + Float_t dDPhi = Jet_phi[ nIdxJet ] - GenJet_phi[ i ]; + + while ( dDPhi > 3.141592 ) dDPhi -= 2 * 3.141592; + while ( dDPhi < -3.141592 ) dDPhi += 2 * 3.141592; + + Float_t dR = sqrt(dDEta * dDEta + dDPhi * dDPhi); + + if ( dRFound > dR ) { + dRFound = dR; + nIdxFound = i; + } + } + + return ( dRFound < 0.4 ? (Int_t)nIdxFound : -1 ); +} + + From a78925c9f4a3701c89af25243c252bdcd27a4a90 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Sun, 16 Dec 2018 23:59:09 +0900 Subject: [PATCH 05/23] For nominal --- analysis/interface/topObjectSelection.h | 1 - analysis/src/topObjectSelection.cc | 13 +++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index c019f981..32058189 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -29,7 +29,6 @@ class topObjectSelection : public nanoBase TRandom3 *rndEngine; public: - Float_t Jet_pt_jer_nom[35]; // YOU MUST SET UP ALL IN THE BELOW!!! // (SetCutValues() will force you to do it) Float_t cut_ElectronPt; diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index 42637c61..d00859ad 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -254,7 +254,8 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, fJetEta = Jet_eta[ nIdx ]; fJetPhi = Jet_phi[ nIdx ]; - if ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { + //if ( m_isMC && ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) + if ( m_isMC ) { // Evaluating the central part cJER of the factor JME::JetParameters jetPars = {{JME::Binning::JetPt, fJetPt}, {JME::Binning::JetEta, fJetEta}, @@ -272,11 +273,11 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, // JER (nominal and up and down) - apply scaling method if matched genJet is found, // apply gaussian smearing method if unmatched - /*if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool - std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) - fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt);*/ - if ( nIdxGen >= 0 ) { - fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt; + if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool + std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) { + fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt); + /*if ( nIdxGen >= 0 ) + fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt;*/ } else { // Different from the postprocess tool, but it might not be important const double smear = rndEngine->Gaus(0, 1); From 9886bd92434b2ad5f8d1fe06e7d41c1bfdac6e98 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Mon, 17 Dec 2018 15:09:16 +0900 Subject: [PATCH 06/23] Sync to the code in CMSSW --- analysis/interface/topObjectSelection.h | 5 +++- analysis/src/topObjectSelection.cc | 35 +++++++++++++++---------- 2 files changed, 25 insertions(+), 15 deletions(-) diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 32058189..48171042 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -72,6 +72,9 @@ class topObjectSelection : public nanoBase Float_t cut_BJetConeSizeOverlap; Float_t *cut_BJetTypeBTag; // For example, set it as cut_BJetTypeBTag = Jet_btagCSVV2; Float_t cut_BJetBTagCut; + + Float_t m_fDRcone_JER; + Float_t m_fResFactorMathcer; public: // Tip: If you want to use your own additional cut with the existing cut, @@ -102,7 +105,7 @@ class topObjectSelection : public nanoBase virtual void GetJetMassPt(UInt_t nIdx, Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi); - Int_t GetMatchGenJet(UInt_t nIdxJet); + Int_t GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution); public: std::vector muonSelection(); diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index d00859ad..42d5f4bd 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -5,7 +5,8 @@ using std::vector; topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, UInt_t unFlag) : nanoBase(tree, had, hadTruth, isMC), m_unFlag(unFlag), jecUnc(NULL), rndEngine(NULL) { - if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn | OptFlag_JER_Up | OptFlag_JER_Dn ) ) != 0 ) { + //if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn | OptFlag_JER_Up | OptFlag_JER_Dn ) ) != 0 ) + if ( m_isMC ) { std::string env = getenv("CMSSW_BASE"); std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jetunc/" @@ -27,6 +28,9 @@ topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, jecUnc = new JetCorrectionUncertainty(JetCorPar); } } + + m_fDRcone_JER = 0.4; // For AK4 jets + m_fResFactorMathcer = 3; // According to https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetResolution } vector topObjectSelection::elecSelection() { @@ -261,23 +265,23 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, {JME::Binning::JetEta, fJetEta}, {JME::Binning::Rho, fixedGridRhoFastjetAll}}; - // We need corresponding genJet - //Int_t nIdxGen = Jet_genJetIdx[ nIdx ]; - Int_t nIdxGen = GetMatchGenJet(nIdx); - const double genJetPt = GenJet_pt[ nIdxGen ]; - const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. const double cJER = jetResSFObj.getScaleFactor(jetPars, ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn ) ) == 0 ? Variation::NOMINAL : ( ( m_unFlag & OptFlag_JER_Up ) != 0 ? Variation::UP : Variation::DOWN ) )); + // We need corresponding genJet + //Int_t nIdxGen = Jet_genJetIdx[ nIdx ]; + Int_t nIdxGen = GetMatchGenJet(nIdx, fJetPt * jetRes); + const double genJetPt = GenJet_pt[ nIdxGen ]; + // JER (nominal and up and down) - apply scaling method if matched genJet is found, // apply gaussian smearing method if unmatched - if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool - std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) { - fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt); - /*if ( nIdxGen >= 0 ) - fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt;*/ + /*if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool + std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) + fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt);*/ + if ( nIdxGen >= 0 ) { + fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt; } else { // Different from the postprocess tool, but it might not be important const double smear = rndEngine->Gaus(0, 1); @@ -302,10 +306,10 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, } -Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet) { +Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { UInt_t i; - Float_t dRFound = 999; + Float_t dRFound = m_fDRcone_JER; UInt_t nIdxFound = 999; for ( i = 0 ; i < nGenJet ; i++ ) { @@ -317,13 +321,16 @@ Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet) { Float_t dR = sqrt(dDEta * dDEta + dDPhi * dDPhi); + if ( dR >= m_fDRcone_JER * 0.5 ) continue; if ( dRFound > dR ) { + if ( std::abs(Jet_pt[ nIdxJet ] - GenJet_pt[ i ]) >= m_fResFactorMathcer * fResolution ) continue; + dRFound = dR; nIdxFound = i; } } - return ( dRFound < 0.4 ? (Int_t)nIdxFound : -1 ); + return ( dRFound < m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); } From 0c745421afc5c6907be55dbb3513ac9048b764e9 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Mon, 17 Dec 2018 19:29:03 +0900 Subject: [PATCH 07/23] Changed the paths for JER/JES text files to a better place --- analysis/src/topObjectSelection.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index 42d5f4bd..d1760542 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -9,18 +9,18 @@ topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, if ( m_isMC ) { std::string env = getenv("CMSSW_BASE"); - std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jetunc/" + std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jer/" "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; jetResSFObj = JME::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); - std::string strPathJetResObj = env + "/src/nano/analysis/data/jetunc/" + std::string strPathJetResObj = env + "/src/nano/analysis/data/jer/" "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; jetResObj = JME::JetResolution(strPathJetResObj.c_str()); rndEngine = new TRandom3(12345); if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { - std::string strPathJecUnc = env + "/src/nano/analysis/data/jetunc/" + std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt"; "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; From 0c2144310921834340e60200229aa4f048bd8ed8 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 10:18:45 +0900 Subject: [PATCH 08/23] Using deltaR provided in CMSSW --- analysis/src/topObjectSelection.cc | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index d1760542..f61a7935 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -1,3 +1,4 @@ +#include "DataFormats/Math/interface/deltaR.h" #include "nano/analysis/interface/topObjectSelection.h" using std::vector; @@ -283,7 +284,6 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, if ( nIdxGen >= 0 ) { fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt; } else { - // Different from the postprocess tool, but it might not be important const double smear = rndEngine->Gaus(0, 1); fCorrFactor = ( cJER <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJER * cJER - 1) ); } @@ -309,17 +309,12 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { UInt_t i; + Float_t dR; Float_t dRFound = m_fDRcone_JER; UInt_t nIdxFound = 999; for ( i = 0 ; i < nGenJet ; i++ ) { - Float_t dDEta = Jet_eta[ nIdxJet ] - GenJet_eta[ i ]; - Float_t dDPhi = Jet_phi[ nIdxJet ] - GenJet_phi[ i ]; - - while ( dDPhi > 3.141592 ) dDPhi -= 2 * 3.141592; - while ( dDPhi < -3.141592 ) dDPhi += 2 * 3.141592; - - Float_t dR = sqrt(dDEta * dDEta + dDPhi * dDPhi); + dR = deltaR(Jet_eta[ nIdxJet ], Jet_phi[ nIdxJet ], GenJet_eta[ i ], GenJet_phi[ i ]); if ( dR >= m_fDRcone_JER * 0.5 ) continue; if ( dRFound > dR ) { @@ -330,7 +325,7 @@ Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { } } - return ( dRFound < m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); + return ( dRFound < 0.75 * m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); } From d4bce5b548a55cdbde0573fe1e9acccdbe6ab417 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 11:25:07 +0900 Subject: [PATCH 09/23] I forgot to remove the previous evaluator in nanoAOD production... --- nanoAOD/python/JetMetUncertainty_cff.py | 13 ------------- nanoAOD/python/nano_cff.py | 12 ++---------- 2 files changed, 2 insertions(+), 23 deletions(-) delete mode 100644 nanoAOD/python/JetMetUncertainty_cff.py diff --git a/nanoAOD/python/JetMetUncertainty_cff.py b/nanoAOD/python/JetMetUncertainty_cff.py deleted file mode 100644 index 71e87cf3..00000000 --- a/nanoAOD/python/JetMetUncertainty_cff.py +++ /dev/null @@ -1,13 +0,0 @@ -import FWCore.ParameterSet.Config as cms -from PhysicsTools.NanoAOD.common_cff import * - -jetUncEvaluatorTable = cms.EDProducer("JetUncertaintyEvaluator", - src = cms.InputTag("linkedObjects","jets"), - rho = cms.InputTag('fixedGridRhoFastjetAll'), - payloadName = cms.string('AK4PFchs'), - jetResFile = cms.string("nano/nanoAOD/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"), - jetResSFFile = cms.string("nano/nanoAOD/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt"), -) - -jetUncEvaluatorTables = cms.Sequence(jetUncEvaluatorTable) - diff --git a/nanoAOD/python/nano_cff.py b/nanoAOD/python/nano_cff.py index e6d95a18..7532b698 100644 --- a/nanoAOD/python/nano_cff.py +++ b/nanoAOD/python/nano_cff.py @@ -22,17 +22,9 @@ def customise(process, doHadron=True, fastSim=False): # DATA process.NANOAODoutput.fileName = fileName jecFile = 'Summer16_07Aug2017All_V10_DATA' - - process.RandomNumberGeneratorService.jetUncEvaluatorTable = cms.PSet( - engineName = cms.untracked.string('TRandom3'), - initialSeed = cms.untracked.uint32(1), - ) - + customiseMuons(process) - process.load('nano.nanoAOD.JetMetUncertainty_cff') - process.nanoAOD_step += process.jetUncEvaluatorTable - process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") if doHadron: process.load('nano.nanoAOD.jetMLID_cff') @@ -73,7 +65,7 @@ def customise(process, doHadron=True, fastSim=False): ) process.es_prefer_jec = cms.ESPrefer("PoolDBESSource","jec") print "JEC based on", process.jec.connect - + process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(1000) process.MessageLogger.cerr.FwkSummary.reportEvery = cms.untracked.int32(1000) From 4810b9ce212b0de9896d38b6cfdcfc294c45eb31 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 11:36:18 +0900 Subject: [PATCH 10/23] Ah, plugin... --- nanoAOD/plugins/JetUncertaintyEvaluator.cc | 142 --------------------- nanoAOD/plugins/JetUncertaintyEvaluator.h | 67 ---------- 2 files changed, 209 deletions(-) delete mode 100644 nanoAOD/plugins/JetUncertaintyEvaluator.cc delete mode 100644 nanoAOD/plugins/JetUncertaintyEvaluator.h diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.cc b/nanoAOD/plugins/JetUncertaintyEvaluator.cc deleted file mode 100644 index 89c479c7..00000000 --- a/nanoAOD/plugins/JetUncertaintyEvaluator.cc +++ /dev/null @@ -1,142 +0,0 @@ -#include "JetUncertaintyEvaluator.h" - - -//#define debugMode -using namespace edm; -using namespace std; - - -JetUncertaintyEvaluator::JetUncertaintyEvaluator(const edm::ParameterSet &iConfig) : - src_(consumes>(iConfig.getParameter("src"))), - rhoToken_(consumes(iConfig.getParameter("rho"))), - payloadName_(iConfig.getParameter("payloadName")), - jetResFilePath_(edm::FileInPath(iConfig.getParameter("jetResFile")).fullPath()), - jetResSFFilePath_(edm::FileInPath(iConfig.getParameter("jetResSFFile")).fullPath()) -{ - produces(); -} - - -void JetUncertaintyEvaluator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; - desc.setUnknown(); - descriptions.addDefault(desc); -} - - -void JetUncertaintyEvaluator::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) { - edm::Service rng; - rng_ = &rng->getEngine(lumi.index()); -} - - -void JetUncertaintyEvaluator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { - runOnMC_ = !iEvent.isRealData(); - - Handle jets; - iEvent.getByToken(src_, jets); - - // The following codes are from: - // https://github.com/vallot/CATTools/blob/cat80x/CatProducer/plugins/CATJetProducer.cc - edm::Handle rhoHandle; - iEvent.getByToken(rhoToken_, rhoHandle); - const double rho = *rhoHandle; - - if ( !payloadName_.empty() ) { - // temp measure - payloadName should be AK4PFchs, but PHYS14_25_V2 does not have uncertainty - edm::ESHandle JetCorParColl; - iSetup.get().get(payloadName_, JetCorParColl); - JetCorrectorParameters const & JetCorPar = ( *JetCorParColl )[ "Uncertainty" ]; - - jecUnc = new JetCorrectionUncertainty(JetCorPar); - } - - JME::JetResolution jetResObj; - JME::JetResolutionScaleFactor jetResSFObj; - - if ( runOnMC_ ) { - jetResObj = JME::JetResolution(jetResFilePath_); - jetResSFObj = JME::JetResolutionScaleFactor(jetResSFFilePath_); - } - - std::vector vecJetJESUp; - std::vector vecJetJESDn; - std::vector vecJetJERUp; - std::vector vecJetJERDn; - - for ( unsigned int i = 0 ; i < jets->size() ; i++ ) { - const auto & jet = ( *jets )[ i ]; - - double jetPt = jet.pt(); - double jetEta = jet.eta(); - - // Computing JEC uncertainty - float fJESUp = 1, fJESDn = 1; - - if ( !payloadName_.empty() ) { - jecUnc->setJetEta(jetEta); - jecUnc->setJetPt(jetPt); // here you must use the CORRECTED jet pt - fJESUp = 1 + jecUnc->getUncertainty(true); - - jecUnc->setJetEta(jetEta); - jecUnc->setJetPt(jetPt); // here you must use the CORRECTED jet pt - fJESDn = 1 - jecUnc->getUncertainty(false); - } - - // Computing JER uncertainty - float fJERUp = 1, fJERDn = 1; - - if ( runOnMC_ ) { - // adding genJet - auto genJet = jet.genJetFwdRef(); - - JME::JetParameters jetPars = {{JME::Binning::JetPt, jetPt}, - {JME::Binning::JetEta, jetEta}, - {JME::Binning::Rho, rho}}; - const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. - const double cJERUp = jetResSFObj.getScaleFactor(jetPars, Variation::UP); - const double cJERDn = jetResSFObj.getScaleFactor(jetPars, Variation::DOWN); - - // JER - apply scaling method if matched genJet is found, - // apply gaussian smearing method if unmatched - if ( genJet.isNonnull() && deltaR(genJet->p4(), jet.p4()) < 0.2 - && std::abs(genJet->pt() - jetPt) < jetRes * 3 * jetPt ) - { - const double genJetPt = genJet->pt(); - const double dPt = jetPt - genJetPt; - - fJERUp = std::max(0., (genJetPt + dPt * cJERUp) / jetPt); - fJERDn = std::max(0., (genJetPt + dPt * cJERDn) / jetPt); - } else { - const double smear = CLHEP::RandGaussQ::shoot(rng_); - - fJERUp = ( cJERUp <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJERUp * cJERUp - 1) ); - fJERDn = ( cJERDn <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJERDn * cJERDn - 1) ); - } - } - - vecJetJESUp.push_back(fJESUp); - vecJetJESDn.push_back(fJESDn); - vecJetJERUp.push_back(fJERUp); - vecJetJERDn.push_back(fJERDn); - } - - auto outJetCor = std::make_unique(jets->size(), "Jet", false, true); - //outJetCor->setDoc("Jet energy scale factor uncertainty"); - - outJetCor->addColumn("jes_up", vecJetJESUp, "Scale factor for JES up", - nanoaod::FlatTable::FloatColumn); - outJetCor->addColumn("jes_dn", vecJetJESDn, "Scale factor for JES down", - nanoaod::FlatTable::FloatColumn); - - outJetCor->addColumn("jer_up", vecJetJERUp, "Scale factor for JER up", - nanoaod::FlatTable::FloatColumn); - outJetCor->addColumn("jer_dn", vecJetJERDn, "Scale factor for JER down", - nanoaod::FlatTable::FloatColumn); - - iEvent.put(move(outJetCor)); -} - -DEFINE_FWK_MODULE(JetUncertaintyEvaluator); - - diff --git a/nanoAOD/plugins/JetUncertaintyEvaluator.h b/nanoAOD/plugins/JetUncertaintyEvaluator.h deleted file mode 100644 index bce99b5a..00000000 --- a/nanoAOD/plugins/JetUncertaintyEvaluator.h +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef JetUncertaintyEvaluator_H -#define JetUncertaintyEvaluator_H - -#include - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/stream/EDProducer.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/Utilities/interface/RandomNumberGenerator.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "FWCore/Utilities/interface/StreamID.h" -#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" -#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" - -#include "CommonTools/Utils/interface/StringCutObjectSelector.h" - -#include "DataFormats/NanoAOD/interface/FlatTable.h" -#include "DataFormats/Common/interface/ValueMap.h" - -#include "DataFormats/PatCandidates/interface/Jet.h" -#include "DataFormats/PatCandidates/interface/PackedCandidate.h" - -#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" -#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" -#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" - -#include "JetMETCorrections/Modules/interface/JetResolution.h" - -#include "CLHEP/Random/RandomEngine.h" -#include "CLHEP/Random/RandGaussQ.h" - -//#define debugMode - -using namespace edm; -using namespace std; - - -class JetUncertaintyEvaluator: public edm::stream::EDProducer<> { -public: - explicit JetUncertaintyEvaluator(const edm::ParameterSet &iConfig); - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - -private: - void beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup&) override; - void produce( edm::Event&, const edm::EventSetup& ) override; - - edm::EDGetTokenT src_; - - edm::EDGetTokenT rhoToken_; - std::string payloadName_; - const std::string jetResFilePath_, jetResSFFilePath_; - - bool runOnMC_; - - JetCorrectionUncertainty *jecUnc; - CLHEP::HepRandomEngine* rng_; -}; - -#endif // JetUncertaintyEvaluater_H - - From 137084c39080c861e481fba090568f4946c42ecc Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 21:12:53 +0900 Subject: [PATCH 11/23] Ported the JER/JES codes in CMSSW to be standalone --- analysis/interface/topObjectSelection.h | 10 +- analysis/src/topObjectSelection.cc | 9 +- external/interface/JetCorrectionUncertainty.h | 72 ++ external/interface/JetCorrectorParameters.h | 271 ++++++++ .../interface/JetCorrectorParametersHelper.h | 63 ++ external/interface/JetResolution.h | 88 +++ external/interface/JetResolutionObject.h | 295 ++++++++ external/src/JetCorrectionUncertainty.cc | 304 +++++++++ external/src/JetCorrectorParameters.cc | 633 ++++++++++++++++++ external/src/JetCorrectorParametersHelper.cc | 210 ++++++ external/src/JetResolution.cc | 83 +++ external/src/JetResolutionObject.cc | 440 ++++++++++++ 12 files changed, 2472 insertions(+), 6 deletions(-) create mode 100644 external/interface/JetCorrectionUncertainty.h create mode 100644 external/interface/JetCorrectorParameters.h create mode 100644 external/interface/JetCorrectorParametersHelper.h create mode 100644 external/interface/JetResolution.h create mode 100644 external/interface/JetResolutionObject.h create mode 100644 external/src/JetCorrectionUncertainty.cc create mode 100644 external/src/JetCorrectorParameters.cc create mode 100644 external/src/JetCorrectorParametersHelper.cc create mode 100644 external/src/JetResolution.cc create mode 100644 external/src/JetResolutionObject.cc diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 48171042..2e8a609c 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -5,10 +5,12 @@ #include "nano/external/interface/TopTriggerSF.h" //#include "nano/external/interface/TTbarModeDefs.h" -#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" -#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" -#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" -#include "JetMETCorrections/Modules/interface/JetResolution.h" +//#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" +//#include "JetMETCorrections/Modules/interface/JetResolution.h" +#include "nano/external/interface/JetCorrectionUncertainty.h" +#include "nano/external/interface/JetCorrectorParameters.h" +#include "nano/external/interface/JetResolution.h" #include class topObjectSelection : public nanoBase diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index f61a7935..d0adfac2 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -309,12 +309,17 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { UInt_t i; - Float_t dR; + Float_t dEta, dPhi, dR; Float_t dRFound = m_fDRcone_JER; UInt_t nIdxFound = 999; for ( i = 0 ; i < nGenJet ; i++ ) { - dR = deltaR(Jet_eta[ nIdxJet ], Jet_phi[ nIdxJet ], GenJet_eta[ i ], GenJet_phi[ i ]); + dEta = Jet_eta[ nIdxJet ] - GenJet_eta[ i ]; + dPhi = std::abs(Jet_eta[ nIdxJet ] - GenJet_eta[ i ]); + while ( dPhi > 3.14159265359 ) dPhi -= 2 * 3.14159265359; + + //dR = deltaR(Jet_eta[ nIdxJet ], Jet_phi[ nIdxJet ], GenJet_eta[ i ], GenJet_phi[ i ]); + dR = sqrt(dEta * dEta + dPhi * dPhi); if ( dR >= m_fDRcone_JER * 0.5 ) continue; if ( dRFound > dR ) { diff --git a/external/interface/JetCorrectionUncertainty.h b/external/interface/JetCorrectionUncertainty.h new file mode 100644 index 00000000..c34a6952 --- /dev/null +++ b/external/interface/JetCorrectionUncertainty.h @@ -0,0 +1,72 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h +// in 4d7b6ef commit +// Copied by Byeonghak Ko +// + + +#ifndef JetCorrectionUncertainty_h +#define JetCorrectionUncertainty_h + +//#include "CondFormats/Serialization/interface/Serializable.h" // MODIFIED + + +// MODIFIED +#define COND_SERIALIZABLE +#define COND_TRANSIENT + + +#include +#include +class SimpleJetCorrectionUncertainty; +class JetCorrectorParameters; + +class JetCorrectionUncertainty +{ + public: + JetCorrectionUncertainty(); + JetCorrectionUncertainty(const std::string& fDataFile); + JetCorrectionUncertainty(const JetCorrectorParameters& fParameters); + ~JetCorrectionUncertainty(); + + void setParameters (const std::string& fDataFile); + void setJetEta (float fEta); + void setJetPt (float fPt); + void setJetE (float fE); + void setJetPhi (float fE); + void setJetEMF (float fEMF); + void setLepPx (float fLepPx); + void setLepPy (float fLepPy); + void setLepPz (float fLepPz); + void setAddLepToJet (bool fAddLepToJet) {mAddLepToJet = fAddLepToJet;} + float getUncertainty(bool fDirection); + + private: + JetCorrectionUncertainty(const JetCorrectionUncertainty&) = delete; + JetCorrectionUncertainty& operator= (const JetCorrectionUncertainty&) = delete; + std::vector fillVector(const std::vector& fNames); + float getPtRel(); + //---- Member Data --------- + float mJetE; + float mJetEta; + float mJetPt; + float mJetPhi; + float mJetEMF; + float mLepPx; + float mLepPy; + float mLepPz; + bool mAddLepToJet; + bool mIsJetEset; + bool mIsJetPtset; + bool mIsJetPhiset; + bool mIsJetEtaset; + bool mIsJetEMFset; + bool mIsLepPxset; + bool mIsLepPyset; + bool mIsLepPzset; + SimpleJetCorrectionUncertainty* mUncertainty; +}; + +#endif + diff --git a/external/interface/JetCorrectorParameters.h b/external/interface/JetCorrectorParameters.h new file mode 100644 index 00000000..ded2b78e --- /dev/null +++ b/external/interface/JetCorrectorParameters.h @@ -0,0 +1,271 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetCorrectorParameters.h +// in 01870b1 commit +// Copied by Byeonghak Ko +// + + +// +// Original Author: Fedor Ratnikov Nov 9, 2007 +// $Id: JetCorrectorParameters.h,v 1.15 2012/03/01 18:24:52 srappocc Exp $ +// +// Generic parameters for Jet corrections +// +#ifndef JetCorrectorParameters_h +#define JetCorrectorParameters_h + + +//#include "CondFormats/Serialization/interface/Serializable.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/Utilities.h" // MODIFIED +#include "nano/external/interface/Utilities.h" // MODIFIED + + +// MODIFIED +#define COND_SERIALIZABLE +#define COND_TRANSIENT + + +#include +#include +#include +#include +#include +#include +#include +//#include "FWCore/Utilities/interface/Exception.h" // MODIFIED +//#include "FWCore/MessageLogger/interface/MessageLogger.h" // MODIFIED + +class JetCorrectorParametersHelper; + +class JetCorrectorParameters +{ + //---------------- JetCorrectorParameters class ---------------- + //-- Encapsulates all the information of the parametrization --- + public: + //---------------- Definitions class --------------------------- + //-- Global iformation about the parametrization is kept here -- + class Definitions + { + public: + //-------- Constructors -------------- + Definitions() : mIsResponse(false) {} + Definitions(const std::vector& fBinVar, const std::vector& fParVar, const std::string& fFormula, bool fIsResponse); + Definitions(const std::string& fLine); + //-------- Member functions ---------- + unsigned nBinVar() const {return mBinVar.size(); } + unsigned nParVar() const {return mParVar.size(); } + std::vector parVar() const {return mParVar; } + std::vector binVar() const {return mBinVar; } + std::string parVar(unsigned fIndex) const {return mParVar[fIndex];} + std::string binVar(unsigned fIndex) const {return mBinVar[fIndex];} + std::string formula() const {return mFormula; } + std::string level() const {return mLevel; } + bool isResponse() const {return mIsResponse; } + private: + //-------- Member variables ---------- + bool mIsResponse; + std::string mLevel; + std::string mFormula; + std::vector mParVar; + std::vector mBinVar; + + COND_SERIALIZABLE; + }; + //---------------- Record class -------------------------------- + //-- Each Record holds the properties of a bin ----------------- + class Record + { + public: + //-------- Constructors -------------- + Record() : mNvar(0),mMin(0),mMax(0) {} + Record(unsigned fNvar, const std::vector& fXMin, const std::vector& fXMax, const std::vector& fParameters) : mNvar(fNvar),mMin(fXMin),mMax(fXMax),mParameters(fParameters) {} + Record(const std::string& fLine, unsigned fNvar); + //-------- Member functions ---------- + unsigned nVar() const {return mNvar; } + float xMin(unsigned fVar) const {return mMin[fVar]; } + float xMax(unsigned fVar) const {return mMax[fVar]; } + float xMiddle(unsigned fVar) const {return 0.5*(xMin(fVar)+xMax(fVar));} + float parameter(unsigned fIndex) const {return mParameters[fIndex]; } + std::vector parameters() const {return mParameters; } + unsigned nParameters() const {return mParameters.size(); } + bool operator< (const Record& other) const + { + if (xMin(0) < other.xMin(0)) return true; + if (xMin(0) > other.xMin(0)) return false; + if (xMin(1) < other.xMin(1)) return true; + if (xMin(1) > other.xMin(1)) return false; + return (xMin(2) < other.xMin(2)); + } + private: + //-------- Member variables ---------- + unsigned mNvar; + std::vector mMin; + std::vector mMax; + std::vector mParameters; + + COND_SERIALIZABLE; + }; + + //-------- Constructors -------------- + JetCorrectorParameters() { valid_ = false;} + JetCorrectorParameters(const std::string& fFile, const std::string& fSection = ""); + JetCorrectorParameters(const JetCorrectorParameters::Definitions& fDefinitions, + const std::vector& fRecords) + : mDefinitions(fDefinitions),mRecords(fRecords) { valid_ = true;} + //-------- Member functions ---------- + const Record& record(unsigned fBin) const {return mRecords[fBin]; } + const Definitions& definitions() const {return mDefinitions; } + unsigned size() const {return mRecords.size();} + unsigned size(unsigned fVar) const; + int binIndex(const std::vector& fX) const; + int binIndexN(const std::vector& fX) const; + int neighbourBin(unsigned fIndex, unsigned fVar, bool fNext) const; + std::vector binCenters(unsigned fVar) const; + void printScreen() const; + void printFile(const std::string& fFileName) const; + bool isValid() const { return valid_; } + void init(); + + static const int MAX_SIZE_DIMENSIONALITY = 3 COND_TRANSIENT; + + private: + //-------- Member variables ---------- + JetCorrectorParameters::Definitions mDefinitions; + std::vector mRecords; + bool valid_; /// is this a valid set? + + std::shared_ptr helper COND_TRANSIENT; + + COND_SERIALIZABLE; +}; +std::ostream& operator<<(std::ostream& out, const JetCorrectorParameters::Record& fBin); + + +class JetCorrectorParametersCollection { + //---------------- JetCorrectorParametersCollection class ---------------- + //-- Adds several JetCorrectorParameters together by algorithm type --- + //-- to reduce the number of payloads in the Database --- + //-- NB: The enum is listed in "human-logical" order, but the actual + // enum value is in order of appearance when people thought of them. + public: + enum Level_t { L1Offset=0, + L1JPTOffset=7, + L1FastJet = 10, + L2Relative=1, + L3Absolute=2, + L2L3Residual=8, + L4EMF=3, + L5Flavor=4, + L6UE=5, + L7Parton=6, + Uncertainty=9, + UncertaintyAbsolute=11, + UncertaintyHighPtExtra=12, + UncertaintySinglePionECAL=13, + UncertaintySinglePionHCAL=27, + UncertaintyFlavor=14, + UncertaintyTime=15, + UncertaintyRelativeJEREC1=16, + UncertaintyRelativeJEREC2=17, + UncertaintyRelativeJERHF=18, + UncertaintyRelativePtEC1=28, + UncertaintyRelativePtEC2=29, + UncertaintyRelativePtHF=30, + UncertaintyRelativeStatEC2=19, + UncertaintyRelativeStatHF=20, + UncertaintyRelativeFSR=21, + UncertaintyRelativeSample=31, + UncertaintyPileUpDataMC=22, + UncertaintyPileUpOOT=23, + UncertaintyPileUpPtBB=24, + UncertaintyPileUpPtEC=32, + UncertaintyPileUpPtHF=33, + UncertaintyPileUpBias=25, + UncertaintyPileUpJetRate=26, + L1RC=34, + L1Residual=35, + UncertaintyAux3=36, + UncertaintyAux4=37, + N_LEVELS=38 + }; + + + enum L5_Species_t {L5_bJ=0,L5_cJ,L5_qJ,L5_gJ,L5_bT,L5_cT,L5_qT,L5_gT,N_L5_SPECIES}; + enum L7_Species_t {L7_gJ=0,L7_qJ,L7_cJ,L7_bJ,L7_jJ,L7_qT,L7_cT,L7_bT,L7_jT,N_L7_SPECIES}; + typedef int key_type; + typedef std::string label_type; + typedef JetCorrectorParameters value_type; + typedef std::pair pair_type; + typedef std::vector collection_type; + + + // Constructor... initialize all three vectors to zero + JetCorrectorParametersCollection() { corrections_.clear(); correctionsL5_.clear(); correctionsL7_.clear(); } + + // Add a JetCorrectorParameter object, possibly with flavor. + void push_back( key_type i, value_type const & j, label_type const & flav = "" ); + + // Access the JetCorrectorParameter via the key k. + // key_type is hashed to deal with the three collections + JetCorrectorParameters const & operator[]( key_type k ) const; + + // Access the JetCorrectorParameter via a string. + // Will find the hashed value for the label, and call via that + // operator. + JetCorrectorParameters const & operator[]( std::string const & label ) const { + return operator[]( findKey(label) ); + } + + // Get a list of valid keys. These will contain hashed keys + // that are aware of all three collections. + void validKeys(std::vector & keys ) const; + + + + // Helper method to find all of the sections in a given + // parameters file + static void getSections(std::string inputFile, + std::vector & outputs ); + + // Find the L5 bin for hashing + static key_type getL5Bin( std::string const & flav ); + // Find the L7 bin for hashing + static key_type getL7Bin( std::string const & flav ); + // Check if this is an L5 hashed value + static bool isL5( key_type k ); + // Check if this is an L7 hashed value + static bool isL7( key_type k ); + + static std::string findLabel( key_type k ); + static std::string findL5Flavor( key_type k ); + static std::string findL7Parton( key_type k ); + + protected: + + // Find the key corresponding to each label + key_type findKey( std::string const & label ) const; + + collection_type corrections_; + collection_type correctionsL5_; + collection_type correctionsL7_; + + collection_type& getCorrections() {return corrections_;} + collection_type& getCorrectionsL5() {return correctionsL5_;} + collection_type& getCorrectionsL7() {return correctionsL7_;} + + friend struct JetCorrectorParametersInitializeTransients; + + COND_SERIALIZABLE; + +}; + +struct JetCorrectorParametersInitializeTransients { + void operator()(JetCorrectorParametersCollection& jcpc) { + for (auto & ptype : jcpc.getCorrections()) {ptype.second.init();} + for (auto & ptype : jcpc.getCorrectionsL5()) {ptype.second.init();} + for (auto & ptype : jcpc.getCorrectionsL7()) {ptype.second.init();} + } +}; + +#endif diff --git a/external/interface/JetCorrectorParametersHelper.h b/external/interface/JetCorrectorParametersHelper.h new file mode 100644 index 00000000..136e3ff3 --- /dev/null +++ b/external/interface/JetCorrectorParametersHelper.h @@ -0,0 +1,63 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetCorrectorParametersHelper.h +// in 01870b1 commit +// Copied by Byeonghak Ko +// + + +// +// Original Author: Fedor Ratnikov Nov 9, 2007 +// $Id: JetCorrectorParameters.h,v 1.15 2012/03/01 18:24:52 srappocc Exp $ +// +// Generic parameters for Jet corrections +// +#ifndef JetCorrectorParametersHelper_h +#define JetCorrectorParametersHelper_h + +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/Utilities.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParameters.h" // MODIFIED +#include "nano/external/interface/Utilities.h" // MODIFIED + +#include +#include +#include +#include +#include +#include +#include // MODIFIED +//#include "FWCore/Utilities/interface/Exception.h" +//#include "FWCore/MessageLogger/interface/MessageLogger.h" + + +//---------- JetCorrectorParametersHelper class ---------------- +//-- The helper is used to find the correct Record to access --- +class JetCorrectorParametersHelper +{ + public: + //-------- Member functions ---------- + unsigned size() const {return mIndexMap.size();} + void initTransientMaps(); + void init(const JetCorrectorParameters::Definitions& mDefinitions, + const std::vector& mRecords); + void checkMiddleBinUniformity(const std::vector& mRecords) const; + void binIndexChecks(unsigned N, const std::vector& fX) const; + bool binBoundChecks(unsigned dim, const float& value, const float& min, const float& max) const; + int binIndexN(const std::vector& fX, const std::vector& mRecords) const; + + using tuple_type = typename generate_tuple_type::type; + using tuple_type_Nm1 = typename generate_tuple_type::type; + private: + //-------- Member variables ---------- + // Stores the lower and upper bounds of the bins for each binned dimension + std::vector > mBinBoundaries; + // Maps a set of lower bounds for N binned dimensions to the index in mRecords + std::unordered_map mIndexMap; + // Maps a set of lower bounds for the first N-1 dimensions to the range of lower bound indices mBinBoundaries for the N dimension + std::unordered_map > mMap; + // The number of binned dimensions as given by the JetCorrectorParameters::Definitions class + unsigned SIZE; +}; + +#endif diff --git a/external/interface/JetResolution.h b/external/interface/JetResolution.h new file mode 100644 index 00000000..7ad735b8 --- /dev/null +++ b/external/interface/JetResolution.h @@ -0,0 +1,88 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/JetMETCorrections/Modules/interface/JetResolution.h +// in a4c8a25 commit +// Copied by Byeonghak Ko +// + + +#ifndef JetResolution_h +#define JetResolution_h + + +#define STANDALONE // MODIFIED + + +// If you want to use the JER code in standalone mode, you'll need to create a new define named 'STANDALONE'. If you use gcc for compiling, you'll need to add +// -DSTANDALONE to the command line +// In standalone mode, no reference to CMSSW exists, so the only way to retrieve resolutions and scale factors are from text files. + +//#include // MODIFIED +#include "nano/external/interface/JetResolutionObject.h" // MODIFIED + +#ifndef STANDALONE +namespace edm { + class EventSetup; +} +#endif + + +namespace JME { + class JetResolution { + public: + JetResolution(const std::string& filename); + JetResolution(const JetResolutionObject& object); + JetResolution() { + // Empty + } + +#ifndef STANDALONE + static const JetResolution get(const edm::EventSetup&, const std::string&); +#endif + + float getResolution(const JetParameters& parameters) const; + + void dump() const { + m_object->dump(); + } + + // Advanced usage + const JetResolutionObject* getResolutionObject() const { + return m_object.get(); + } + + private: + std::shared_ptr m_object; + }; + + class JetResolutionScaleFactor { + public: + JetResolutionScaleFactor(const std::string& filename); + JetResolutionScaleFactor(const JetResolutionObject& object); + JetResolutionScaleFactor() { + // Empty + } + +#ifndef STANDALONE + static const JetResolutionScaleFactor get(const edm::EventSetup&, const std::string&); +#endif + + float getScaleFactor(const JetParameters& parameters, Variation variation = Variation::NOMINAL) const; + + void dump() const { + m_object->dump(); + } + + // Advanced usage + const JetResolutionObject* getResolutionObject() const { + return m_object.get(); + } + + private: + std::shared_ptr m_object; + }; + +}; + +#endif + diff --git a/external/interface/JetResolutionObject.h b/external/interface/JetResolutionObject.h new file mode 100644 index 00000000..aad1f3c3 --- /dev/null +++ b/external/interface/JetResolutionObject.h @@ -0,0 +1,295 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetResolutionObject.h +// in 28c2e9f commit +// Copied by Byeonghak Ko +// + + +#ifndef JetResolutionObject_h +#define JetResolutionObject_h + +#define STANDALONE // MODIFIED - For standalone of nanoAOD analysers + +// If you want to use the JER code in standalone mode, you'll need to create a new define named 'STANDALONE'. If you use gcc for compiling, you'll need to add +// -DSTANDALONE to the command line +// In standalone mode, no reference to CMSSW exists, so the only way to retrieve resolutions and scale factors are from text files. + +#ifndef STANDALONE +#include "CondFormats/Serialization/interface/Serializable.h" +#else +// Create no-op definitions of CMSSW macro +#define COND_SERIALIZABLE +#define COND_TRANSIENT +#endif + +#include +#include +#include +#include +#include +#include + +#ifndef STANDALONE +#include "CommonTools/Utils/interface/FormulaEvaluator.h" +#else +#include +#endif + +enum class Variation { + NOMINAL = 0, + DOWN = 1, + UP = 2 +}; + +template +T clip(const T& n, const T& lower, const T& upper) { + return std::max(lower, std::min(n, upper)); +} + +namespace JME { + template + struct bimap { + typedef std::unordered_map left_type; + typedef std::unordered_map right_type; + + left_type left; + right_type right; + + bimap(std::initializer_list l) { + for (auto& v: l) { + left.insert(v); + right.insert(typename right_type::value_type(v.second, v.first)); + } + } + + bimap() { + // Empty + } + + bimap(bimap&& rhs) { + left = std::move(rhs.left); + right = std::move(rhs.right); + } + }; + + enum class Binning { + JetPt = 0, + JetEta, + JetAbsEta, + JetE, + JetArea, + Mu, + Rho, + NPV, + }; + +}; + +// Hash function for Binning enum class +namespace std { + template<> + struct hash { + typedef JME::Binning argument_type; + typedef std::size_t result_type; + + hash int_hash; + + result_type operator()(argument_type const& s) const { + return int_hash(static_cast(s)); + } + }; +}; + +namespace JME { + + class JetParameters { + public: + typedef std::unordered_map value_type; + + JetParameters() = default; + JetParameters(JetParameters&& rhs); + JetParameters(std::initializer_list init); + + + JetParameters& setJetPt(float pt); + JetParameters& setJetEta(float eta); + JetParameters& setJetE(float e); + JetParameters& setJetArea(float area); + JetParameters& setMu(float mu); + JetParameters& setRho(float rho); + JetParameters& setNPV(float npv); + JetParameters& set(const Binning& bin, float value); + JetParameters& set(const typename value_type::value_type& value); + + static const bimap binning_to_string; + + std::vector createVector(const std::vector& binning) const; + + private: + value_type m_values; + }; + + + class JetResolutionObject { + + public: + + struct Range { + float min; + float max; + + Range() { + // Empty + } + + Range(float min, float max) { + this->min = min; + this->max = max; + } + + bool is_inside(float value) const { + return (value >= min) && (value <= max); + } + + COND_SERIALIZABLE; + }; + + class Definition { + + public: + Definition() { + // Empty + } + + Definition(const std::string& definition); + + const std::vector& getBinsName() const { + return m_bins_name; + } + + const std::vector& getBins() const { + return m_bins; + } + + std::string getBinName(size_t bin) const { + return m_bins_name[bin]; + } + + size_t nBins() const { + return m_bins_name.size(); + } + + const std::vector& getVariablesName() const { + return m_variables_name; + } + + const std::vector& getVariables() const { + return m_variables; + } + + std::string getVariableName(size_t variable) const { + return m_variables_name[variable]; + } + + size_t nVariables() const { + return m_variables.size(); + } + + std::string getFormulaString() const { + return m_formula_str; + } + +#ifndef STANDALONE + const reco::FormulaEvaluator* getFormula() const { + return m_formula.get(); + } +#else + TFormula const * getFormula() const { + return m_formula.get(); + } +#endif + void init(); + + private: + std::vector m_bins_name; + std::vector m_variables_name; + std::string m_formula_str; + +#ifndef STANDALONE + std::shared_ptr m_formula COND_TRANSIENT; +#else + std::shared_ptr m_formula COND_TRANSIENT; +#endif + std::vector m_bins COND_TRANSIENT; + std::vector m_variables COND_TRANSIENT; + + COND_SERIALIZABLE; + }; + + class Record { + public: + Record() { + // Empty + } + + Record(const std::string& record, const Definition& def); + + const std::vector& getBinsRange() const { + return m_bins_range; + } + + const std::vector& getVariablesRange() const { + return m_variables_range; + } + + const std::vector& getParametersValues() const { + return m_parameters_values; + } + + size_t nVariables() const { + return m_variables_range.size(); + } + + size_t nParameters() const { + return m_parameters_values.size(); + } + + private: + std::vector m_bins_range; + std::vector m_variables_range; + std::vector m_parameters_values; + + COND_SERIALIZABLE; + }; + + public: + JetResolutionObject(const std::string& filename); + JetResolutionObject(const JetResolutionObject& filename); + JetResolutionObject(); + + void dump() const; + void saveToFile(const std::string& file) const; + + const Record* getRecord(const JetParameters& bins) const; + float evaluateFormula(const Record& record, const JetParameters& variables) const; + + const std::vector& getRecords() const { + return m_records; + } + + const Definition& getDefinition() const { + return m_definition; + } + + private: + Definition m_definition; + std::vector m_records; + + bool m_valid = false; + + COND_SERIALIZABLE; + }; +}; + +#endif diff --git a/external/src/JetCorrectionUncertainty.cc b/external/src/JetCorrectionUncertainty.cc new file mode 100644 index 00000000..a9539cfe --- /dev/null +++ b/external/src/JetCorrectionUncertainty.cc @@ -0,0 +1,304 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/JetCorrectionUncertainty.cc +// in deb815b commit +// Copied by Byeonghak Ko +// + + +//#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" // MODIFIED +#include "nano/external/interface/JetCorrectionUncertainty.h" // MODIFIED +#include "nano/external/interface/SimpleJetCorrectionUncertainty.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParameters.h" // MODIFIED + +//#include "FWCore/MessageLogger/interface/MessageLogger.h" // MODIFIED + +//#include "Math/PtEtaPhiE4D.h" // MODIFIED +//#include "Math/Vector3D.h" // MODIFIED +//#include "Math/LorentzVector.h" // MODIFIED +#include // MODIFIED +#include // MODIFIED +#include +#include + + +// MODIFIED +namespace edm { + inline std::ostream &LogError(const char *szTitle) {return std::cerr << szTitle << ": ";} +}; + + +///////////////////////////////////////////////////////////////////////// +JetCorrectionUncertainty::JetCorrectionUncertainty () +{ + mJetEta = -9999; + mJetPt = -9999; + mJetPhi = -9999; + mJetE = -9999; + mJetEMF = -9999; + mLepPx = -9999; + mLepPy = -9999; + mLepPz = -9999; + mIsJetEset = false; + mIsJetPtset = false; + mIsJetPhiset = false; + mIsJetEtaset = false; + mIsJetEMFset = false; + mIsLepPxset = false; + mIsLepPyset = false; + mIsLepPzset = false; + mAddLepToJet = false; + mUncertainty = new SimpleJetCorrectionUncertainty(); +} +///////////////////////////////////////////////////////////////////////// +JetCorrectionUncertainty::JetCorrectionUncertainty(const std::string& fDataFile) +{ + mJetEta = -9999; + mJetPt = -9999; + mJetPhi = -9999; + mJetE = -9999; + mJetEMF = -9999; + mLepPx = -9999; + mLepPy = -9999; + mLepPz = -9999; + mIsJetEset = false; + mIsJetPtset = false; + mIsJetPhiset = false; + mIsJetEtaset = false; + mIsJetEMFset = false; + mIsLepPxset = false; + mIsLepPyset = false; + mIsLepPzset = false; + mAddLepToJet = false; + mUncertainty = new SimpleJetCorrectionUncertainty(fDataFile); +} +///////////////////////////////////////////////////////////////////////// +JetCorrectionUncertainty::JetCorrectionUncertainty(const JetCorrectorParameters& fParameters) +{ + mJetEta = -9999; + mJetPt = -9999; + mJetPhi = -9999; + mJetE = -9999; + mJetEMF = -9999; + mLepPx = -9999; + mLepPy = -9999; + mLepPz = -9999; + mIsJetEset = false; + mIsJetPtset = false; + mIsJetPhiset = false; + mIsJetEtaset = false; + mIsJetEMFset = false; + mIsLepPxset = false; + mIsLepPyset = false; + mIsLepPzset = false; + mAddLepToJet = false; + mUncertainty = new SimpleJetCorrectionUncertainty(fParameters); +} +///////////////////////////////////////////////////////////////////////// +JetCorrectionUncertainty::~JetCorrectionUncertainty () +{ + delete mUncertainty; +} +///////////////////////////////////////////////////////////////////////// +void JetCorrectionUncertainty::setParameters(const std::string& fDataFile) +{ + //---- delete the mParameters pointer before setting the new address --- + delete mUncertainty; + mUncertainty = new SimpleJetCorrectionUncertainty(fDataFile); +} +///////////////////////////////////////////////////////////////////////// +float JetCorrectionUncertainty::getUncertainty(bool fDirection) +{ + float result; + std::vector vx,vy; + vx = fillVector(mUncertainty->parameters().definitions().binVar()); + vy = fillVector(mUncertainty->parameters().definitions().parVar()); + result = mUncertainty->uncertainty(vx,vy[0],fDirection); + mIsJetEset = false; + mIsJetPtset = false; + mIsJetPhiset = false; + mIsJetEtaset = false; + mIsJetEMFset = false; + mIsLepPxset = false; + mIsLepPyset = false; + mIsLepPzset = false; + return result; +} +//------------------------------------------------------------------------ +//--- Reads the parameter names and fills a vector of floats ------------- +//------------------------------------------------------------------------ +std::vector JetCorrectionUncertainty::fillVector(const std::vector& fNames) +{ + std::vector result; + for(unsigned i=0;i > PtEtaPhiELorentzVector; + //typedef ROOT::Math::DisplacementVector3D > XYZVector; + //PtEtaPhiELorentzVector jet; + //XYZVector lep; + //jet.SetPt(mJetPt); + //jet.SetEta(mJetEta); + //jet.SetPhi(mJetPhi); + //jet.SetE(mJetE); + //lep.SetXYZ(mLepPx,mLepPy,mLepPz); + TLorentzVector jet; + TVector3 lep; + jet.SetPtEtaPhiE(mJetPt, mJetEta, mJetPhi, mJetE); + lep.SetXYZ(mLepPx,mLepPy,mLepPz); + float lj_x = (mAddLepToJet) ? lep.X()+jet.Px() : jet.Px(); + float lj_y = (mAddLepToJet) ? lep.Y()+jet.Py() : jet.Py(); + float lj_z = (mAddLepToJet) ? lep.Z()+jet.Pz() : jet.Pz(); + // absolute values squared + float lj2 = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z; + float pTrel2 = -999.0; + if (lj2 > 0) { + float lep2 = lep.X()*lep.X()+lep.Y()*lep.Y()+lep.Z()*lep.Z(); + // projection vec(mu) to lepjet axis + float lepXlj = lep.X()*lj_x+lep.Y()*lj_y+lep.Z()*lj_z; + // absolute value squared and normalized + float pLrel2 = lepXlj*lepXlj/lj2; + // lep2 = pTrel2 + pLrel2 + pTrel2 = lep2-pLrel2; + } else + edm::LogError("JetCorrectionUncertainty")<<" not positive lepton-jet momentum: "< 0) ? std::sqrt(pTrel2) : 0.0; +} +//------------------------------------------------------------------------ +//--- Setters ------------------------------------------------------------ +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setJetEta(float fEta) +{ + mJetEta = fEta; + mIsJetEtaset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setJetPt(float fPt) +{ + mJetPt = fPt; + mIsJetPtset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setJetPhi(float fPhi) +{ + mJetPhi = fPhi; + mIsJetPhiset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setJetE(float fE) +{ + mJetE = fE; + mIsJetEset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setJetEMF(float fEMF) +{ + mJetEMF = fEMF; + mIsJetEMFset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setLepPx(float fPx) +{ + mLepPx = fPx; + mIsLepPxset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setLepPy(float fPy) +{ + mLepPy = fPy; + mIsLepPyset = true; +} +//------------------------------------------------------------------------ +void JetCorrectionUncertainty::setLepPz(float fPz) +{ + mLepPz = fPz; + mIsLepPzset = true; +} +//------------------------------------------------------------------------ diff --git a/external/src/JetCorrectorParameters.cc b/external/src/JetCorrectorParameters.cc new file mode 100644 index 00000000..7fbc4a92 --- /dev/null +++ b/external/src/JetCorrectorParameters.cc @@ -0,0 +1,633 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/JetCorrectorParameters.cc +// in 4d7b6ef commit +// Copied by Byeonghak Ko +// + + +// +// Original Author: Fedor Ratnikov Nov 9, 2007 +// $Id: JetCorrectorParameters.cc,v 1.20 2012/03/01 18:24:53 srappocc Exp $ +// +// Generic parameters for Jet corrections +// +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParametersHelper.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParameters.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParametersHelper.h" // MODIFIED +#include +#include +#include +#include +#include +#include +#include +#include + +//------------------------------------------------------------------------ +//--- JetCorrectorParameters::Definitions constructor -------------------- +//--- takes specific arguments for the member variables ------------------ +//------------------------------------------------------------------------ +JetCorrectorParameters::Definitions::Definitions(const std::vector& fBinVar, const std::vector& fParVar, const std::string& fFormula, bool fIsResponse) +{ + for(unsigned i=0;i tokens = getTokens(fLine); + if (!tokens.empty()) + { + if (tokens.size() < 6) + { + std::stringstream sserr; + sserr<<"(line "< tokens = getTokens(fLine); + if (!tokens.empty()) + { + if (tokens.size() < 3) + { + std::stringstream sserr; + sserr<<"(line "<(); + helper->init(mDefinitions,mRecords); +} +//------------------------------------------------------------------------ +//--- returns the index of the record defined by fX ---------------------- +//------------------------------------------------------------------------ +int JetCorrectorParameters::binIndex(const std::vector& fX) const +{ + int result = -1; + unsigned N = mDefinitions.nBinVar(); + if (N != fX.size()) + { + std::stringstream sserr; + sserr<<"# bin variables "<= record(i).xMin(j) && fX[j] < record(i).xMax(j)) + tmp+=1; + if (tmp==N) + { + result = i; + break; + } + } + return result; +} +//------------------------------------------------------------------------ +//--- returns the index of the record defined by fX (non-linear search) -- +//------------------------------------------------------------------------ +int JetCorrectorParameters::binIndexN(const std::vector& fX) const +{ + if(helper->size()>0) + { + return helper->binIndexN(fX,mRecords); + } + else + { + return -1; + } +} +//------------------------------------------------------------------------ +//--- returns the neighbouring bins of fIndex in the direction of fVar --- +//------------------------------------------------------------------------ +int JetCorrectorParameters::neighbourBin(unsigned fIndex, unsigned fVar, bool fNext) const +{ + int result = -1; + unsigned N = mDefinitions.nBinVar(); + if (fVar >= N) + { + std::stringstream sserr; + sserr<<"# of bin variables "<= mDefinitions.nBinVar()) + { + std::stringstream sserr; + sserr<<"requested bin variable index "< tmpMin && record(i).xMax(fVar) > tmpMax) + { + result++; + tmpMin = record(i).xMin(fVar); + tmpMax = record(i).xMax(fVar); + } + return result; +} +//------------------------------------------------------------------------ +//--- returns the vector of bin centers of fVar -------------------------- +//------------------------------------------------------------------------ +std::vector JetCorrectorParameters::binCenters(unsigned fVar) const +{ + std::vector result; + for (unsigned i = 0; i < size(); ++i) + result.push_back(record(i).xMiddle(fVar)); + return result; +} +//------------------------------------------------------------------------ +//--- prints parameters on screen ---------------------------------------- +//------------------------------------------------------------------------ +void JetCorrectorParameters::printScreen() const +{ + std::cout<<"--------------------------------------------"< labels_ = { + "L1Offset", + "L2Relative", + "L3Absolute", + "L4EMF", + "L5Flavor", + "L6UE", + "L7Parton", + "L1JPTOffset", + "L2L3Residual", + "Uncertainty", + "L1FastJet", + "UncertaintyAbsolute", + "UncertaintyHighPtExtra", + "UncertaintySinglePionECAL", + "UncertaintyFlavor", + "UncertaintyTime", + "UncertaintyRelativeJEREC1", + "UncertaintyRelativeJEREC2", + "UncertaintyRelativeJERHF", + "UncertaintyRelativeStatEC2", + "UncertaintyRelativeStatHF", + "UncertaintyRelativeFSR", + "UncertaintyPileUpDataMC", + "UncertaintyPileUpOOT", + "UncertaintyPileUpPtBB", + "UncertaintyPileUpBias", + "UncertaintyPileUpJetRate", + "UncertaintySinglePionHCAL", + "UncertaintyRelativePtEC1", + "UncertaintyRelativePtEC2", + "UncertaintyRelativePtHF", + "UncertaintyRelativeSample", + "UncertaintyPileUpPtEC", + "UncertaintyPileUpPtHF", + "L1RC", + "L1Residual", + "UncertaintyAux3", + "UncertaintyAux4", +}; + +const std::vector l5Flavors_ = { + "L5Flavor_bJ", + "L5Flavor_cJ", + "L5Flavor_qJ", + "L5Flavor_gJ", + "L5Flavor_bT", + "L5Flavor_cT", + "L5Flavor_qT", + "L5Flavor_gT" +}; + +const std::vector l7Partons_ = { + "L7Parton_gJ", + "L7Parton_qJ", + "L7Parton_cJ", + "L7Parton_bJ", + "L7Parton_jJ", + "L7Parton_qT", + "L7Parton_cT", + "L7Parton_bT", + "L7Parton_tT" +}; +} +std::string +JetCorrectorParametersCollection::findLabel( key_type k ){ + if ( isL5(k) ) return findL5Flavor(k); + else if ( isL7(k) ) return findL7Parton(k); + else return labels_[k]; +} + +std::string +JetCorrectorParametersCollection::findL5Flavor( key_type k ){ + if ( k == L5Flavor ) return labels_[L5Flavor]; + else + return l5Flavors_[k / 100 - 1]; +} + +std::string +JetCorrectorParametersCollection::findL7Parton( key_type k ){ + if ( k == L7Parton ) return labels_[L7Parton]; + else + return l7Partons_[k / 1000 - 1]; +} + +void JetCorrectorParametersCollection::getSections(std::string inputFile, + std::vector & outputs ) +{ + outputs.clear(); + std::ifstream input( inputFile.c_str() ); + while( !input.eof() ) { + char buff[10000]; + input.getline(buff,10000); + std::string in(buff); + if ( in[0] == '[' ) { + std::string tok = getSection(in); + if ( tok != "" ) { + outputs.push_back( tok ); + } + } + } + std::cout << "Found these sections for file: " << std::endl; + copy(outputs.begin(),outputs.end(), std::ostream_iterator(std::cout, "\n") ); +} + + +// Add a JetCorrectorParameter object, possibly with flavor. +void JetCorrectorParametersCollection::push_back( key_type i, value_type const & j, label_type const & flav) { + std::cout << "i = " << i << std::endl; + std::cout << "flav = " << flav << std::endl; + if ( isL5(i) ) { + std::cout << "This is L5, getL5Bin = " << getL5Bin(flav) << std::endl; + correctionsL5_.push_back( pair_type(getL5Bin(flav),j) ); + } + else if ( isL7(i) ) { + std::cout << "This is L7, getL7Bin = " << getL7Bin(flav) << std::endl; + correctionsL7_.push_back( pair_type(getL7Bin(flav),j) ); + } + else if ( flav == "" ) { + corrections_.push_back( pair_type(i,j) ); + } else { + std::cout << "***** NOT ADDING " << flav << ", corresponding position in JetCorrectorParameters is not found." << std::endl; + } +} + + +// Access the JetCorrectorParameter via the key k. +// key_type is hashed to deal with the three collections +JetCorrectorParameters const & JetCorrectorParametersCollection::operator[]( key_type k ) const { + collection_type::const_iterator ibegin, iend, i; + if ( isL5(k) ) { + ibegin = correctionsL5_.begin(); + iend = correctionsL5_.end(); + i = ibegin; + } else if ( isL7(k) ) { + ibegin = correctionsL7_.begin(); + iend = correctionsL7_.end(); + i = ibegin; + } else { + ibegin = corrections_.begin(); + iend = corrections_.end(); + i = ibegin; + } + for ( ; i != iend; ++i ) { + if ( k == i->first ) return i->second; + } + // MODIFIED + //throw cms::Exception("InvalidInput") << " cannot find key " << static_cast(k) + // << " in the JEC payload, this usually means you have to change the global tag" << std::endl; + std::string strException; + std::stringstream(strException) << "InvalidInput: " + << " cannot find key " << static_cast(k) + << " in the JEC payload, this usually means you have to change the global tag" << std::endl; + throw std::runtime_error(strException); +} + +// Get a list of valid keys. These will contain hashed keys +// that are aware of all three collections. +void JetCorrectorParametersCollection::validKeys(std::vector & keys ) const { + keys.clear(); + for ( collection_type::const_iterator ibegin = corrections_.begin(), + iend = corrections_.end(), i = ibegin; i != iend; ++i ) { + keys.push_back( i->first ); + } + for ( collection_type::const_iterator ibegin = correctionsL5_.begin(), + iend = correctionsL5_.end(), i = ibegin; i != iend; ++i ) { + keys.push_back( i->first ); + } + for ( collection_type::const_iterator ibegin = correctionsL7_.begin(), + iend = correctionsL7_.end(), i = ibegin; i != iend; ++i ) { + keys.push_back( i->first ); + } +} + + +// Find the L5 bin for hashing +JetCorrectorParametersCollection::key_type +JetCorrectorParametersCollection::getL5Bin( std::string const & flav ){ + std::vector::const_iterator found = + find( l5Flavors_.begin(), l5Flavors_.end(), flav ); + if ( found != l5Flavors_.end() ) { + return (found - l5Flavors_.begin() + 1) * 100; + } + else return L5Flavor; +} +// Find the L7 bin for hashing +JetCorrectorParametersCollection::key_type +JetCorrectorParametersCollection::getL7Bin( std::string const & flav ){ + std::vector::const_iterator found = + find( l7Partons_.begin(), l7Partons_.end(), flav ); + if ( found != l7Partons_.end() ) { + return (found - l7Partons_.begin() + 1) * 1000; + } + else return L7Parton; +} + +// Check if this is an L5 hashed value +bool JetCorrectorParametersCollection::isL5( key_type k ) { + return k == L5Flavor || + ( k / 100 > 0 && k / 1000 == 0 ); +} +// Check if this is an L7 hashed value +bool JetCorrectorParametersCollection::isL7( key_type k ) { + return k == L7Parton || + ( k / 1000 > 0 ); +} + + +// Find the key corresponding to each label +JetCorrectorParametersCollection::key_type +JetCorrectorParametersCollection::findKey( std::string const & label ) const { + + // First check L5 corrections + std::vector::const_iterator found1 = + find( l5Flavors_.begin(), l5Flavors_.end(), label ); + if ( found1 != l5Flavors_.end() ) { + return getL5Bin(label); + } + + // Next check L7 corrections + std::vector::const_iterator found2 = + find( l7Partons_.begin(), l7Partons_.end(), label ); + if ( found2 != l7Partons_.end() ) { + return getL7Bin(label); + } + + // Finally check the default corrections + std::vector::const_iterator found3 = + find( labels_.begin(), labels_.end(), label ); + if ( found3 != labels_.end() ) { + return static_cast(found3 - labels_.begin()); + } + + // Didn't find default corrections, throw exception + // MODIFIED + //throw cms::Exception("InvalidInput") << " Cannot find label " << label << std::endl; + std::string strException; + std::stringstream(strException) << "InvalidInput: " << " Cannot find label " << label << std::endl; + throw std::runtime_error(strException); +} + + +//#include "FWCore/Framework/interface/EventSetup.h" +//#include "FWCore/Framework/interface/ESHandle.h" +//#include "FWCore/Framework/interface/ModuleFactory.h" +//#include "FWCore/Utilities/interface/typelookup.h" // MODIFIED + +//TYPELOOKUP_DATA_REG(JetCorrectorParameters); // MODIFIED +//TYPELOOKUP_DATA_REG(JetCorrectorParametersCollection); // MODIFIED diff --git a/external/src/JetCorrectorParametersHelper.cc b/external/src/JetCorrectorParametersHelper.cc new file mode 100644 index 00000000..4c52376a --- /dev/null +++ b/external/src/JetCorrectorParametersHelper.cc @@ -0,0 +1,210 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/JetCorrectorParametersHelper.cc +// in 4d7b6ef commit +// Copied by Byeonghak Ko +// + + +// +// Original Author: Alexx Perloff Feb 22, 2017 +// $Id: JetCorrectorParameters.cc,v 1.20 2012/03/01 18:24:53 srappocc Exp $ +// +// Helper class for JetCorrectorParameters +// +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParametersHelper.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParametersHelper.h" // MODIFIED +#include +#include +#include +#include +#include + +//------------------------------------------------------------------------ +//--- JetCorrectorParameters::JetCorrectorParametersHelper initializer --- +//--- initializes the mBinMap for quick lookup of mRecords --------------- +//------------------------------------------------------------------------ +void JetCorrectorParametersHelper::initTransientMaps() +{ + mIndexMap.clear(); + mMap.clear(); + mBinBoundaries.clear(); + mBinBoundaries.assign(JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY,std::vector(0,0)); +} +void JetCorrectorParametersHelper::init(const JetCorrectorParameters::Definitions& mDefinitionsLocal, + const std::vector& mRecordsLocal) +{ + SIZE = mDefinitionsLocal.nBinVar(); + if (SIZE>JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY) + { + std::stringstream sserr; + sserr<<"The number of binned variables requested ("<1 && (i==nRec-1 || mRecordsLocal[i].xMin(j-1)!=mRecordsLocal[i+1].xMin(j-1))) { + end = i; + //mMap.emplace(gen_tuple([&](size_t k){return mRecordsLocal[i].xMin(k);}),std::make_pair(start,end)); + mMap.emplace(gen_tuple([&](size_t k){return (k([&](size_t k){return (ksecond; + std::stringstream sserr; + sserr<<"Duplicate binning in record found (existing index,current index)=(" + <2) checkMiddleBinUniformity(mRecordsLocal); +} + +void JetCorrectorParametersHelper::checkMiddleBinUniformity(const std::vector& mRecords) const +{ + unsigned N = SIZE-2; + size_t nRec = mRecords.size(); + std::vector fN(N,-1); + //The order of looping (records or dimensions) does not matter because you have to go through all of them once anyway + //Loop over each binned dimension that isn't the first or the last + for (unsigned idim=1; idim<=N; idim++) { + int nBoundaryPassed = 0; + if (fN[idim-1]==-1) fN[idim-1] = mBinBoundaries[idim].size(); + //Loop over the mRecords vector + for (unsigned iRecord=0; iRecord& fX) const +{ + if (N != fX.size()) + { + std::stringstream sserr; + sserr<<"The number of binned variables, "< max) return false; + else return true; +} +//------------------------------------------------------------------------ +//--- JetCorrectorParameters::JetCorrectorParametersHelper binIndexN ----- +//--- returns the index of the record defined by fX (non-linear search) -- +//------------------------------------------------------------------------ +int JetCorrectorParametersHelper::binIndexN(const std::vector& fX, + const std::vector& mRecords) const +{ + unsigned Nm1 = SIZE-1; + binIndexChecks(SIZE,fX); + + //Create a container for the indices + std::vector fN(SIZE,-1); + std::vector::const_iterator tmpIt; + + // make sure that fX are within the first and last boundaries of mBinBoundaries (other than last dimension) + for (unsigned idim=0; idim==0||idim indexBounds; + if(SIZE>1) + { + tuple_type_Nm1 to_find_Nm1 = gen_tuple([&](size_t i){return (i([&](size_t i){return (i +#include +#include + +#include +#include +#else +//#include "JetResolution.h" // MODIFIED +#include "nano/external/interface/JetResolution.h" // MODIFIED +#endif + +namespace JME { + + JetResolution::JetResolution(const std::string& filename) { + m_object = std::make_shared(filename); + } + + JetResolution::JetResolution(const JetResolutionObject& object) { + m_object = std::make_shared(object); + } + +#ifndef STANDALONE + const JetResolution JetResolution::get(const edm::EventSetup& setup, const std::string& label) { + edm::ESHandle handle; + setup.get().get(label, handle); + + return *handle; + } +#endif + + float JetResolution::getResolution(const JetParameters& parameters) const { + const JetResolutionObject::Record* record = m_object->getRecord(parameters); + if (! record) + return 1; + + return m_object->evaluateFormula(*record, parameters); + } + + JetResolutionScaleFactor::JetResolutionScaleFactor(const std::string& filename) { + m_object = std::make_shared(filename); + } + + JetResolutionScaleFactor::JetResolutionScaleFactor(const JetResolutionObject& object) { + m_object = std::make_shared(object); + } + +#ifndef STANDALONE + const JetResolutionScaleFactor JetResolutionScaleFactor::get(const edm::EventSetup& setup, const std::string& label) { + edm::ESHandle handle; + setup.get().get(label, handle); + + return *handle; + } +#endif + + float JetResolutionScaleFactor::getScaleFactor(const JetParameters& parameters, Variation variation/* = Variation::NOMINAL*/) const { + const JetResolutionObject::Record* record = m_object->getRecord(parameters); + if (! record) + return 1; + + const std::vector& parameters_values = record->getParametersValues(); + return parameters_values[static_cast(variation)]; + } + +} + +#ifndef STANDALONE +#include "FWCore/Utilities/interface/typelookup.h" +TYPELOOKUP_DATA_REG(JME::JetResolution); +TYPELOOKUP_DATA_REG(JME::JetResolutionScaleFactor); +#endif diff --git a/external/src/JetResolutionObject.cc b/external/src/JetResolutionObject.cc new file mode 100644 index 00000000..373251f7 --- /dev/null +++ b/external/src/JetResolutionObject.cc @@ -0,0 +1,440 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetResolutionObject.h +// in a2f9cbb commit +// Copied by Byeonghak Ko +// + + +#define STANDALONE // MODIFIED - For standalone of nanoAOD analysers + +#ifndef STANDALONE +#include +#include +#include + +#else +//#include "JetResolutionObject.h" // MODIFIED +//#include "Utilities.h" // MODIFIED +#include "nano/external/interface/JetResolutionObject.h" // MODIFIED +#include "nano/external/interface/Utilities.h" // MODIFIED +#include + +namespace edm { + namespace errors { + enum ErrorCode { + NotFound = 8026, + ConfigFileReadError = 7002, + UnimplementedFeature = 8011, + FileReadError = 8021 + }; + }; +}; +#endif + +#include +#include +#include +#include +#include + +namespace JME { + + std::string getDefinitionLine(const std::string& line) { + size_t first = line.find ('{'); + size_t last = line.find ('}'); + + if (first != std::string::npos && last != std::string::npos && first < last) + return std::string(line, first + 1, last - first - 1); + + return ""; + } + + void throwException(uint32_t code, const std::string& message) { +#ifndef STANDALONE + throw edm::Exception(static_cast(code), message); +#else + std::stringstream error; + error << message << " Error code: " << code; + throw std::runtime_error(error.str()); + +#endif + } + + const bimap JetParameters::binning_to_string = { + {Binning::JetPt, "JetPt"}, {Binning::JetEta, "JetEta"}, + {Binning::JetAbsEta, "JetAbsEta"}, {Binning::JetE, "JetE"}, + {Binning::JetArea, "JetArea"}, {Binning::Mu, "Mu"}, + {Binning::Rho, "Rho"}, {Binning::NPV, "NPV"} + }; + + JetParameters::JetParameters(JetParameters&& rhs) { + m_values = std::move(rhs.m_values); + } + + JetParameters::JetParameters(std::initializer_list init) { + for (auto& i: init) { + set(i.first, i.second); + } + } + + JetParameters& JetParameters::setJetPt(float pt) { + m_values[Binning::JetPt] = pt; + return *this; + } + + JetParameters& JetParameters::setJetEta(float eta) { + m_values[Binning::JetEta] = eta; + m_values[Binning::JetAbsEta] = fabs(eta); + return *this; + } + + JetParameters& JetParameters::setJetE(float e) { + m_values[Binning::JetE] = e; + return *this; + } + + JetParameters& JetParameters::setJetArea(float area) { + m_values[Binning::JetArea] = area; + return *this; + } + + JetParameters& JetParameters::setMu(float mu) { + m_values[Binning::Mu] = mu; + return *this; + } + + JetParameters& JetParameters::setNPV(float npv) { + m_values[Binning::NPV] = npv; + return *this; + } + + JetParameters& JetParameters::setRho(float rho) { + m_values[Binning::Rho] = rho; + return *this; + } + + JetParameters& JetParameters::set(const Binning& bin, float value) { + m_values.emplace(bin, value); + + // Special case for eta + if (bin == Binning::JetEta) { + m_values.emplace(Binning::JetAbsEta, fabs(value)); + } + + return *this; + } + + JetParameters& JetParameters::set(const typename value_type::value_type& value) { + set(value.first, value.second); + return *this; + } + + std::vector JetParameters::createVector(const std::vector& binning) const { + std::vector values; + for (const auto& bin: binning) { + const auto& it = m_values.find(bin); + if (it == m_values.cend()) { + throwException(edm::errors::NotFound, "JER parametrisation depends on '" + + JetParameters::binning_to_string.left.at(bin) + + "' but no value for this parameter has been specified. Please call the appropriate 'set' function of the JME::JetParameters object"); + } + + values.push_back(it->second); + } + + return values; + } + + + JetResolutionObject::Definition::Definition(const std::string& definition) { + + std::vector tokens = getTokens(definition); + + // We need at least 3 tokens + if (tokens.size() < 3) { + throwException(edm::errors::ConfigFileReadError, "Definition line needs at least three tokens. Please check file format."); + } + + size_t n_bins = std::stoul(tokens[0]); + + if (tokens.size() < (n_bins + 2)) { + throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check."); + } + + for (size_t i = 0; i < n_bins; i++) { + m_bins_name.push_back(tokens[i + 1]); + } + + size_t n_variables = std::stoul(tokens[n_bins + 1]); + + if (tokens.size() < (1 + n_bins + 1 + n_variables + 1)) { + throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check."); + } + + for (size_t i = 0; i < n_variables; i++) { + m_variables_name.push_back(tokens[n_bins + 2 + i]); + } + + m_formula_str = tokens[n_bins + n_variables + 2]; + + std::string formula_str_lower = m_formula_str; + std::transform(formula_str_lower.begin(), formula_str_lower.end(), formula_str_lower.begin(), ::tolower); + + if (formula_str_lower == "none") + m_formula_str = ""; + + init(); + } + + void JetResolutionObject::Definition::init() { + if (!m_formula_str.empty()) +#ifndef STANDALONE + m_formula = std::make_shared(m_formula_str); +#else + m_formula = std::make_shared("jet_resolution_formula", m_formula_str.c_str()); +#endif + for (const auto& bin: m_bins_name) { + const auto& b = JetParameters::binning_to_string.right.find(bin); + if (b == JetParameters::binning_to_string.right.cend()) { + throwException(edm::errors::UnimplementedFeature, "Bin name not supported: '" + bin + "'"); + } + m_bins.push_back(b->second); + } + + for (const auto& v: m_variables_name) { + const auto& var = JetParameters::binning_to_string.right.find(v); + if (var == JetParameters::binning_to_string.right.cend()) { + throwException(edm::errors::UnimplementedFeature, "Variable name not supported: '" + v + "'"); + } + m_variables.push_back(var->second); + } + } + + JetResolutionObject::Record::Record(const std::string& line, const Definition& def) { + + std::vector tokens = getTokens(line); + + if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1)) { + throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line); + } + + size_t pos = 0; + + for (size_t i = 0; i < def.nBins(); i++) { + Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1])); + pos += 2; + m_bins_range.push_back(r); + } + + size_t n_parameters = std::stoul(tokens[pos++]); + + if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1 + (n_parameters - def.nVariables() * 2))) { + throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line); + } + + for (size_t i = 0; i < def.nVariables(); i++) { + Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1])); + pos += 2; + m_variables_range.push_back(r); + n_parameters -= 2; + } + + for (size_t i = 0; i < n_parameters; i++) { + m_parameters_values.push_back(std::stof(tokens[pos++])); + } + } + + JetResolutionObject::JetResolutionObject(const std::string& filename) { + + // Parse file + std::ifstream f(filename); + + if (! f.good()) { + throwException(edm::errors::FileReadError, "Can't read input file '" + filename + "'"); + } + + for (std::string line; std::getline(f, line); ) { + if ((line.empty()) || (line[0] == '#')) + continue; + + std::string definition = getDefinitionLine(line); + + if (!definition.empty()) { + m_definition = Definition(definition); + } else { + m_records.push_back(Record(line, m_definition)); + } + } + + m_valid = true; + } + + JetResolutionObject::JetResolutionObject(const JetResolutionObject& object) { + m_definition = object.m_definition; + m_records = object.m_records; + m_valid = object.m_valid; + + m_definition.init(); + } + + JetResolutionObject::JetResolutionObject() { + // Empty + } + + + void JetResolutionObject::dump() const { + std::cout << "Definition: " << std::endl; + std::cout << " Number of binning variables: " << m_definition.nBins() << std::endl; + std::cout << " "; + for (const auto& bin: m_definition.getBinsName()) { + std::cout << bin << ", "; + } + std::cout << std::endl; + std::cout << " Number of variables: " << m_definition.nVariables() << std::endl; + std::cout << " "; + for (const auto& bin: m_definition.getVariablesName()) { + std::cout << bin << ", "; + } + std::cout << std::endl; + std::cout << " Formula: " << m_definition.getFormulaString() << std::endl; + + std::cout << std::endl << "Bin contents" << std::endl; + + for (const auto& record: m_records) { + std::cout << " Bins" << std::endl; + size_t index = 0; + for (const auto& bin: record.getBinsRange()) { + std::cout << " " << m_definition.getBinName(index) << " [" << bin.min << " - " << bin.max << "]" << std::endl; + index++; + } + + std::cout << " Variables" << std::endl; + index = 0; + for (const auto& r: record.getVariablesRange()) { + std::cout << " " << m_definition.getVariableName(index) << " [" << r.min << " - " << r.max << "] " << std::endl; + index++; + } + + std::cout << " Parameters" << std::endl; + index = 0; + for (const auto& par: record.getParametersValues()) { + std::cout << " Parameter #" << index << " = " << par << std::endl; + index++; + } + } + } + + void JetResolutionObject::saveToFile(const std::string& file) const { + + std::ofstream fout(file); + fout.setf(std::ios::right); + + // Definition + fout << "{" << m_definition.nBins(); + + for (auto& bin: m_definition.getBinsName()) + fout << " " << bin; + + fout << " " << m_definition.nVariables(); + + for (auto& var: m_definition.getVariablesName()) + fout << " " << var; + + fout << " " << (m_definition.getFormulaString().empty() ? "None" : m_definition.getFormulaString()) << " Resolution}" << std::endl; + + // Records + for (auto& record: m_records) { + for (auto& r: record.getBinsRange()) { + fout << std::left << std::setw(15) << r.min << std::setw(15) << r.max << std::setw(15); + } + fout << (record.nVariables() * 2 + record.nParameters()) << std::setw(15); + + for (auto& r: record.getVariablesRange()) { + fout << r.min << std::setw(15) << r.max << std::setw(15); + } + + for (auto& p: record.getParametersValues()) { + fout << p << std::setw(15); + } + + fout << std::endl << std::setw(0); + } + + } + + const JetResolutionObject::Record* JetResolutionObject::getRecord(const JetParameters& bins_parameters) const { + // Find record for bins + if (! m_valid) + return nullptr; + + // Create vector of bins value. Throw if some values are missing + std::vector bins = bins_parameters.createVector(m_definition.getBins()); + + // Iterate over all records, and find the one for which all bins are valid + const Record* good_record = nullptr; + for (const auto& record: m_records) { + + // Iterate over bins + size_t valid_bins = 0; + size_t current_bin = 0; + for (const auto& bin: record.getBinsRange()) { + if (bin.is_inside(bins[current_bin])) + valid_bins++; + + current_bin++; + } + + if (valid_bins == m_definition.nBins()) { + good_record = &record; + break; + } + } + + return good_record; + } + + float JetResolutionObject::evaluateFormula(const JetResolutionObject::Record& record, const JetParameters& variables_parameters) const { + + if (! m_valid) + return 1; + +#ifndef STANDALONE + const auto* formula = m_definition.getFormula(); +#else + // Set parameters + auto const* pFormula = m_definition.getFormula(); + if (! pFormula) + return 1; + auto formula = *pFormula; +#endif + // Create vector of variables value. Throw if some values are missing + std::vector variables = variables_parameters.createVector(m_definition.getVariables()); + + double variables_[4] = {0}; + for (size_t index = 0; index < m_definition.nVariables(); index++) { + variables_[index] = clip(variables[index], record.getVariablesRange()[index].min, record.getVariablesRange()[index].max); + } + const std::vector& parameters = record.getParametersValues(); + +#ifndef STANDALONE + //ArrayAdaptor only takes doubles + std::vector parametersD(parameters.begin(),parameters.end()); + return formula->evaluate( + reco::formula::ArrayAdaptor(variables_,m_definition.nVariables()), + reco::formula::ArrayAdaptor(parametersD.data(),parametersD.size()) + ); +#else + for (size_t index = 0; index < parameters.size(); index++) { + formula.SetParameter(index, parameters[index]); + } + + return formula.EvalPar(variables_); +#endif + } +} + +#ifndef STANDALONE +#include "FWCore/Utilities/interface/typelookup.h" +TYPELOOKUP_DATA_REG(JME::JetResolutionObject); +#endif From 752627db0218ecb75eaea0d79219c94a371d6c00 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 22:20:53 +0900 Subject: [PATCH 12/23] Fixed on deltaR --- analysis/interface/topObjectSelection.h | 3 --- analysis/src/topObjectSelection.cc | 16 +++++++--------- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 2e8a609c..9063aed2 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -5,9 +5,6 @@ #include "nano/external/interface/TopTriggerSF.h" //#include "nano/external/interface/TTbarModeDefs.h" -//#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" -//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" -//#include "JetMETCorrections/Modules/interface/JetResolution.h" #include "nano/external/interface/JetCorrectionUncertainty.h" #include "nano/external/interface/JetCorrectorParameters.h" #include "nano/external/interface/JetResolution.h" diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index d0adfac2..ebd8e7b6 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -1,4 +1,3 @@ -#include "DataFormats/Math/interface/deltaR.h" #include "nano/analysis/interface/topObjectSelection.h" using std::vector; @@ -309,19 +308,18 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { UInt_t i; - Float_t dEta, dPhi, dR; - Float_t dRFound = m_fDRcone_JER; + double dEta, dPhi, dR; + double dRFound = m_fDRcone_JER; UInt_t nIdxFound = 999; for ( i = 0 ; i < nGenJet ; i++ ) { dEta = Jet_eta[ nIdxJet ] - GenJet_eta[ i ]; - dPhi = std::abs(Jet_eta[ nIdxJet ] - GenJet_eta[ i ]); - while ( dPhi > 3.14159265359 ) dPhi -= 2 * 3.14159265359; + dPhi = std::abs(Jet_phi[ nIdxJet ] - GenJet_phi[ i ]); + if ( dPhi > (double)M_PI ) dPhi -= (double)( 2 * M_PI ); - //dR = deltaR(Jet_eta[ nIdxJet ], Jet_phi[ nIdxJet ], GenJet_eta[ i ], GenJet_phi[ i ]); - dR = sqrt(dEta * dEta + dPhi * dPhi); + dR = std::sqrt(dEta * dEta + dPhi * dPhi); - if ( dR >= m_fDRcone_JER * 0.5 ) continue; + if ( dR >= (double)m_fDRcone_JER * 0.5 ) continue; if ( dRFound > dR ) { if ( std::abs(Jet_pt[ nIdxJet ] - GenJet_pt[ i ]) >= m_fResFactorMathcer * fResolution ) continue; @@ -330,7 +328,7 @@ Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { } } - return ( dRFound < 0.75 * m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); + return ( dRFound < 0.75 * (double)m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); } From d7fa1c0287e53828dab7f1d9f2bb9c704d06b3a6 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 22:54:28 +0900 Subject: [PATCH 13/23] Forgot some additional files... --- .../SimpleJetCorrectionUncertainty.h | 44 +++ external/interface/Utilities.h | 324 ++++++++++++++++++ .../src/SimpleJetCorrectionUncertainty.cc | 139 ++++++++ 3 files changed, 507 insertions(+) create mode 100644 external/interface/SimpleJetCorrectionUncertainty.h create mode 100644 external/interface/Utilities.h create mode 100644 external/src/SimpleJetCorrectionUncertainty.cc diff --git a/external/interface/SimpleJetCorrectionUncertainty.h b/external/interface/SimpleJetCorrectionUncertainty.h new file mode 100644 index 00000000..9be9feee --- /dev/null +++ b/external/interface/SimpleJetCorrectionUncertainty.h @@ -0,0 +1,44 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h +// in 4d7b6ef commit +// Copied by Byeonghak Ko +// + + +#ifndef SimpleJetCorrectionUncertainty_h +#define SimpleJetCorrectionUncertainty_h + +//#include "CondFormats/Serialization/interface/Serializable.h" // MODIFIED + + +// MODIFIED +#define COND_SERIALIZABLE +#define COND_TRANSIENT + + +#include +#include +class JetCorrectorParameters; + +class SimpleJetCorrectionUncertainty +{ + public: + SimpleJetCorrectionUncertainty(); + SimpleJetCorrectionUncertainty(const std::string& fDataFile); + SimpleJetCorrectionUncertainty(const JetCorrectorParameters& fParameters); + ~SimpleJetCorrectionUncertainty(); + const JetCorrectorParameters& parameters() const {return *mParameters;} + float uncertainty(const std::vector& fX, float fY, bool fDirection) const; + + private: + SimpleJetCorrectionUncertainty(const SimpleJetCorrectionUncertainty&) = delete; + SimpleJetCorrectionUncertainty& operator= (const SimpleJetCorrectionUncertainty&) = delete; + int findBin(const std::vector& v, float x) const; + float uncertaintyBin(unsigned fBin, float fY, bool fDirection) const; + float linearInterpolation (float fZ, const float fX[2], const float fY[2]) const; + JetCorrectorParameters* mParameters; +}; + +#endif + diff --git a/external/interface/Utilities.h b/external/interface/Utilities.h new file mode 100644 index 00000000..a676aa8b --- /dev/null +++ b/external/interface/Utilities.h @@ -0,0 +1,324 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/Utilities.h +// in 84ba5c3 commit +// Copied by Byeonghak Ko +// + + +#ifndef CondFormats_JetMETObjects_Utilities_h +#define CondFormats_JetMETObjects_Utilities_h + +#define STANDALONE // MODIFIED - For standalone of nanoAOD analysers + +#ifdef STANDALONE +#include +#else +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#endif + +#include +#include +#include +#include +#include +#include +#include + +namespace std +{ + //These functions print a tuple using a provided std::ostream + template + struct tuple_printer { + + static void print(std::ostream& out, const Type& value) { + out << std::get(value) << ", "; + tuple_printer::print(out, value); + } + }; + template + struct tuple_printer { + + static void print(std::ostream& out, const Type& value) { + out << std::get(value); + } + + }; + template + std::ostream& operator<<(std::ostream& out, const std::tuple& value) { + out << "("; + tuple_printer, 0, sizeof...(Types) - 1>::print(out, value); + out << ")"; + return out; + } + //---------------------------------------------------------------------- + //Returns a list of type indices + template + struct ct_integers_list { + template + struct push_back + { + typedef ct_integers_list type; + }; + }; + template + struct ct_iota_1 + { + typedef typename ct_iota_1::type::template push_back::type type; + }; + template <> + struct ct_iota_1<0> + { + typedef ct_integers_list<> type; + }; + //---------------------------------------------------------------------- + //Return a tuple which is a subset of the original tuple + //This function pops an entry off the font of the tuple + template + auto tuple_subset(const Tuple& tpl, ct_integers_list) + -> decltype(std::make_tuple(std::get(tpl)...)) + { + return std::make_tuple(std::get(tpl)...); + // this means: + // make_tuple(get(tpl), get(tpl), ...) + } + template + std::tuple tuple_tail(const std::tuple& tpl) + { + return tuple_subset(tpl, typename ct_iota_1::type()); + // this means: + // tuple_subset<1, 2, 3, ..., sizeof...(Tail)-1>(tpl, ..) + } + //---------------------------------------------------------------------- + //Recursive hashing function for tuples + template struct hash_specialization + { + typedef std::tuple argument_type; + typedef std::size_t result_type; + result_type operator()(const argument_type& t) const + { + const uint32_t& b = reinterpret_cast(std::get<0>(t)); + //const uint32_t& more = (*this)(tuple_tail(t)); + const uint32_t& more = hash_specialization()(tuple_tail(t)); + return b^more; + } + }; + //Base case + template<> struct hash_specialization + { + typedef std::tuple argument_type; + typedef std::size_t result_type; + result_type operator()(const argument_type& t) const + { + const uint32_t& b = reinterpret_cast(std::get<0>(t)); + const result_type& result = reinterpret_cast(b); + return result; + } + }; + //Overloaded verions of std::hash for tuples + template struct hash > + { + typedef std::tuple argument_type; + typedef std::size_t result_type; + result_type operator()(const argument_type& t) const + { + return hash_specialization()(t); + } + }; + template<> struct hash > + { + typedef std::tuple<> argument_type; + typedef std::size_t result_type; + result_type operator()(const argument_type& t) const + { + return -1; + } + }; +} + +namespace +{ + void handleError(const std::string& fClass, const std::string& fMessage) + { +#ifdef STANDALONE + std::stringstream sserr; + sserr< getTokens(const std::string& fLine) + { + std::vector tokens; + std::string currentToken; + for (unsigned ipos = 0; ipos < fLine.length (); ++ipos) + { + char c = fLine[ipos]; + if (c == '#') break; // ignore comments + else if (c == ' ') + { // flush current token if any + if (!currentToken.empty()) + { + tokens.push_back(currentToken); + currentToken.clear(); + } + } + else + currentToken += c; + } + if (!currentToken.empty()) tokens.push_back(currentToken); // flush end + return tokens; + } + //---------------------------------------------------------------------- + inline std::string getDefinitions(const std::string& token) + { + size_t iFirst = token.find ('{'); + size_t iLast = token.find ('}'); + if (iFirst != std::string::npos && iLast != std::string::npos && iFirst < iLast) + return std::string (token, iFirst+1, iLast-iFirst-1); + return ""; + } + //------------------------------------------------------------------------ + inline float quadraticInterpolation(float fZ, const float fX[3], const float fY[3]) + { + // Quadratic interpolation through the points (x[i],y[i]). First find the parabola that + // is defined by the points and then calculate the y(z). + float D[4],a[3]; + D[0] = fX[0]*fX[1]*(fX[0]-fX[1])+fX[1]*fX[2]*(fX[1]-fX[2])+fX[2]*fX[0]*(fX[2]-fX[0]); + D[3] = fY[0]*(fX[1]-fX[2])+fY[1]*(fX[2]-fX[0])+fY[2]*(fX[0]-fX[1]); + D[2] = fY[0]*(pow(fX[2],2)-pow(fX[1],2))+fY[1]*(pow(fX[0],2)-pow(fX[2],2))+fY[2]*(pow(fX[1],2)-pow(fX[0],2)); + D[1] = fY[0]*fX[1]*fX[2]*(fX[1]-fX[2])+fY[1]*fX[0]*fX[2]*(fX[2]-fX[0])+fY[2]*fX[0]*fX[1]*(fX[0]-fX[1]); + if (D[0] != 0) + { + a[0] = D[1]/D[0]; + a[1] = D[2]/D[0]; + a[2] = D[3]/D[0]; + } + else + { + a[0] = 0.0; + a[1] = 0.0; + a[2] = 0.0; + } + float r = a[0]+fZ*(a[1]+fZ*a[2]); + return r; + } + //------------------------------------------------------------------------ + //Generates a std::tuple type based on a stored type and the number of + // objects in the tuple. + //Note: All of the objects will be of the same type + template + struct join_tuples + { + }; + template + struct join_tuples, std::tuple> + { + typedef std::tuple type; + }; + template + struct generate_tuple_type + { + typedef typename generate_tuple_type::type left; + typedef typename generate_tuple_type::type right; + typedef typename join_tuples::type type; + }; + template + struct generate_tuple_type + { + typedef std::tuple type; + }; + template + struct generate_tuple_type + { + typedef std::tuple<> type; + }; + //------------------------------------------------------------------------ + //C++11 implementation of make_index_sequence, which is a C++14 function + // using aliases for cleaner syntax + template using Invoke = typename T::type; + + template struct seq{ using type = seq; }; + + template struct concat; + + template + struct concat, seq> + : seq{}; + + template + using Concat = Invoke>; + + template struct gen_seq; + template using GenSeq = Invoke>; + + template + struct gen_seq : Concat, GenSeq>{}; + + template<> struct gen_seq<0> : seq<>{}; + template<> struct gen_seq<1> : seq<0>{}; + //------------------------------------------------------------------------ + //Generates a tuple based on a given function (i.e. lambda expression) + template + auto gen_tuple_impl(F func, seq ) + -> decltype(std::make_tuple(func(Is)...)) + { + return std::make_tuple(func(Is)...); + } + template + auto gen_tuple(F func) + -> decltype(gen_tuple_impl(func, GenSeq() )) + { + return gen_tuple_impl(func, GenSeq() ); + } +} +#endif diff --git a/external/src/SimpleJetCorrectionUncertainty.cc b/external/src/SimpleJetCorrectionUncertainty.cc new file mode 100644 index 00000000..a617b0e8 --- /dev/null +++ b/external/src/SimpleJetCorrectionUncertainty.cc @@ -0,0 +1,139 @@ +// +// This code is from copying +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/SimpleJetCorrectionUncertainty.cc +// in d86f5ed commit +// Copied by Byeonghak Ko +// + + +//#include "CondFormats/JetMETObjects/interface/SimpleJetCorrectionUncertainty.h" // MODIFIED +//#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" // MODIFIED +#include "nano/external/interface/SimpleJetCorrectionUncertainty.h" // MODIFIED +#include "nano/external/interface/JetCorrectorParameters.h" // MODIFIED +//#include "FWCore/MessageLogger/interface/MessageLogger.h" // MODIFIED +#include +#include + + +// MODIFIED +namespace edm { + std::ostream &LogError(const char *szTitle) {return std::cerr << szTitle << ": ";} +}; + + +///////////////////////////////////////////////////////////////////////// +SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty () +{ + mParameters = new JetCorrectorParameters(); +} +///////////////////////////////////////////////////////////////////////// +SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const std::string& fDataFile) +{ + mParameters = new JetCorrectorParameters(fDataFile); +} +///////////////////////////////////////////////////////////////////////// +SimpleJetCorrectionUncertainty::SimpleJetCorrectionUncertainty(const JetCorrectorParameters& fParameters) +{ + mParameters = new JetCorrectorParameters(fParameters); +} +///////////////////////////////////////////////////////////////////////// +SimpleJetCorrectionUncertainty::~SimpleJetCorrectionUncertainty () +{ + delete mParameters; +} +///////////////////////////////////////////////////////////////////////// +float SimpleJetCorrectionUncertainty::uncertainty(const std::vector& fX, float fY, bool fDirection) const +{ + float result = 1.; + int bin = mParameters->binIndex(fX); + if (bin<0) { + edm::LogError("SimpleJetCorrectionUncertainty")<<" bin variables out of range"; + result = -999.0; + } else + result = uncertaintyBin((unsigned)bin,fY,fDirection); + return result; +} +///////////////////////////////////////////////////////////////////////// +float SimpleJetCorrectionUncertainty::uncertaintyBin(unsigned fBin, float fY, bool fDirection) const +{ + if (fBin >= mParameters->size()) { + edm::LogError("SimpleJetCorrectionUncertainty")<<" wrong bin: "<size()<<" are available"; + return -999.0; + } + const std::vector& p = mParameters->record(fBin).parameters(); + if ((p.size() % 3) != 0) + // MODIFIED + //throw cms::Exception ("SimpleJetCorrectionUncertainty")<<"wrong # of parameters: multiple of 3 expected, "< yGrid,value; + unsigned int N = p.size()/3; + float result = -1.0; + for(unsigned i=0;i= yGrid[N-1]) + result = value[N-1]; + else + { + int bin = findBin(yGrid,fY); + float vx[2],vy[2]; + for(int i=0;i<2;i++) + { + vx[i] = yGrid[bin+i]; + vy[i] = value[bin+i]; + } + result = linearInterpolation(fY,vx,vy); + } + return result; +} +///////////////////////////////////////////////////////////////////////// +float SimpleJetCorrectionUncertainty::linearInterpolation(float fZ, const float fX[2], const float fY[2]) const +{ + // Linear interpolation through the points (x[i],y[i]). First find the line that + // is defined by the points and then calculate the y(z). + float r = 0; + if (fX[0] == fX[1]) + { + if (fY[0] == fY[1]) + r = fY[0]; + else { + edm::LogError("SimpleJetCorrectionUncertainty")<<" interpolation error"; + return -999.0; + } + } + else + { + float a = (fY[1]-fY[0])/(fX[1]-fX[0]); + float b = (fY[0]*fX[1]-fY[1]*fX[0])/(fX[1]-fX[0]); + r = a*fZ+b; + } + return r; +} +///////////////////////////////////////////////////////////////////////// +int SimpleJetCorrectionUncertainty::findBin(const std::vector& v, float x) const +{ + int i; + int n = v.size()-1; + if (n<=0) return -1; + if (x=v[n]) + return -1; + for(i=0;i=v[i] && x Date: Tue, 18 Dec 2018 23:43:25 +0900 Subject: [PATCH 14/23] I guess Travis wants these files. This is just temporary and these text files must be added into CPLUOS/nanoData. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index a850b0c5..ee813730 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,4 +15,4 @@ script: - wget https://cern.ch/iawatson/mc_test_v5.root - wget https://cern.ch/iawatson/data_test_v5.root - docker exec -it cmssw bash -c 'yum --disablerepo=osg install -y unzip' - - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && cd nano && getFiles && testAnalyser' + - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt nano/analysis/data/jer && cd nano && getFiles && testAnalyser' From 4ddb662557c641006a8ea05bfcec8b9ce854a001 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Tue, 18 Dec 2018 23:56:39 +0900 Subject: [PATCH 15/23] Oops, a mistake... --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index ee813730..a2e3ef8f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,4 +15,4 @@ script: - wget https://cern.ch/iawatson/mc_test_v5.root - wget https://cern.ch/iawatson/data_test_v5.root - docker exec -it cmssw bash -c 'yum --disablerepo=osg install -y unzip' - - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt nano/analysis/data/jer && cd nano && getFiles && testAnalyser' + - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt -P nano/analysis/data/jer && cd nano && getFiles && testAnalyser' From 490f9b3118094c0da3327079d732fe9f0178b616 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 19 Dec 2018 00:16:32 +0900 Subject: [PATCH 16/23] Revert "Oops, a mistake..." This reverts commit 4ddb662557c641006a8ea05bfcec8b9ce854a001. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index a2e3ef8f..ee813730 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,4 +15,4 @@ script: - wget https://cern.ch/iawatson/mc_test_v5.root - wget https://cern.ch/iawatson/data_test_v5.root - docker exec -it cmssw bash -c 'yum --disablerepo=osg install -y unzip' - - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt -P nano/analysis/data/jer && cd nano && getFiles && testAnalyser' + - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt nano/analysis/data/jer && cd nano && getFiles && testAnalyser' From da5d011f159b15543e28bc8ac60e220176548562 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 19 Dec 2018 00:16:38 +0900 Subject: [PATCH 17/23] Revert "I guess Travis wants these files. This is just temporary and these text files must be added into CPLUOS/nanoData.": Okay, seems like Travis doesn't like me >:( This reverts commit acb1ac0ce4ec85a9f31f216c8e34bd6a170744b9. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index ee813730..a850b0c5 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,4 +15,4 @@ script: - wget https://cern.ch/iawatson/mc_test_v5.root - wget https://cern.ch/iawatson/data_test_v5.root - docker exec -it cmssw bash -c 'yum --disablerepo=osg install -y unzip' - - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P nano/analysis/data/jer && wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt nano/analysis/data/jer && cd nano && getFiles && testAnalyser' + - docker exec -it cmssw bash -c 'cd ~/CMSSW_9_4_6/src && ls && source /opt/cms/cmsset_default.sh; eval `scram runtime -sh` && cd nano && getFiles && testAnalyser' From 9877bc363d712fc55dcb363c17cf493de814c7ed Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 19 Dec 2018 00:18:40 +0900 Subject: [PATCH 18/23] Okay, guess Travis wants me to add these files as like this. But these files must be added into CPLUOS/nanoData, as mentioned --- ...mmer16_25nsV1_MC_PtResolution_AK4PFchs.txt | 183 ++++++++++++++++++ .../jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt | 27 +++ 2 files changed, 210 insertions(+) create mode 100644 analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt create mode 100644 analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt diff --git a/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt b/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt new file mode 100644 index 00000000..7633031e --- /dev/null +++ b/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt @@ -0,0 +1,183 @@ +{2 JetEta Rho 1 JetPt sqrt([0]*abs([0])/(x*x)+[1]*[1]*pow(x,[3])+[2]*[2]) Resolution} +-4.7 -3.2 0 6.69 6 15 3000 2.511 0.3167 0.09085 -0.7407 +-4.7 -3.2 6.69 12.39 6 15 3000 3.297 0.2091 6.258e-05 -0.2755 +-4.7 -3.2 12.39 18.09 6 15 3000 1.85 2.281 0.1042 -1.635 +-4.7 -3.2 18.09 23.79 6 15 3000 3.869 1.001 0.09955 -1.266 +-4.7 -3.2 23.79 29.49 6 15 3000 -23.98 24.11 0.1057 -1.988 +-4.7 -3.2 29.49 35.19 6 15 3000 5.403 0.2371 1.5e-05 -0.3177 +-4.7 -3.2 35.19 40.9 6 15 3000 5.753 0.2337 0.0002982 -0.3108 +-3.2 -3 0 6.69 6 15 3000 0.0002851 3.01 0.1382 -1.702 +-3.2 -3 6.69 12.39 6 15 3000 -33.01 33.04 0.1343 -1.991 +-3.2 -3 12.39 18.09 6 15 3000 -67.94 67.8 0.1342 -1.996 +-3.2 -3 18.09 23.79 6 15 3000 -47.81 48 0.1391 -1.996 +-3.2 -3 23.79 29.49 6 15 3000 7.162 0.9211 0.1395 -1.209 +-3.2 -3 29.49 35.19 6 15 3000 8.193 0.1995 2.822e-05 -0.132 +-3.2 -3 35.19 40.9 6 15 3000 8.133 0.9983 0.1349 -1.181 +-3 -2.8 0 6.69 6 15 3000 4.467 0.1997 -3.491e-06 -0.2623 +-3 -2.8 6.69 12.39 6 15 3000 4.17 0.928 0.07702 -1.063 +-3 -2.8 12.39 18.09 6 15 3000 -0.04491 3.67 0.08704 -1.641 +-3 -2.8 18.09 23.79 6 15 3000 5.528 1.286 0.07962 -1.187 +-3 -2.8 23.79 29.49 6 15 3000 -78.36 78.23 0.08448 -1.996 +-3 -2.8 29.49 35.19 6 15 3000 7.559 1.147 0.07023 -1.134 +-3 -2.8 35.19 40.9 6 15 3000 -59.03 59.03 -0.08184 -1.992 +-2.8 -2.5 0 6.69 6 15 3000 4.244 0.2766 -1.86e-08 -0.5068 +-2.8 -2.5 6.69 12.39 6 15 3000 4.919 0.3193 5.463e-06 -0.58 +-2.8 -2.5 12.39 18.09 6 15 3000 5.909 0.2752 4.144e-06 -0.5272 +-2.8 -2.5 18.09 23.79 6 15 3000 -47.31 47.18 0.05853 -1.991 +-2.8 -2.5 23.79 29.49 6 15 3000 -46.49 46.33 0.05698 -1.989 +-2.8 -2.5 29.49 35.19 6 15 3000 8.651 0.2522 6.592e-06 -0.4835 +-2.8 -2.5 35.19 40.9 6 15 3000 7.716 2.481 0.0531 -1.455 +-2.5 -2.3 0 6.69 6 15 3000 3.125 0.6026 0.02576 -0.8702 +-2.5 -2.3 6.69 12.39 6 15 3000 3.935 0.6533 0.02587 -0.889 +-2.5 -2.3 12.39 18.09 6 15 3000 4.198 1.024 0.03618 -1.069 +-2.5 -2.3 18.09 23.79 6 15 3000 2.948 2.386 0.04771 -1.382 +-2.5 -2.3 23.79 29.49 6 15 3000 4.415 2.086 0.04704 -1.294 +-2.5 -2.3 29.49 35.19 6 15 3000 -3.084 4.156 0.05366 -1.503 +-2.5 -2.3 35.19 40.9 6 15 3000 -6.144 5.969 0.05633 -1.602 +-2.3 -2.1 0 6.69 6 15 3000 0.3022 1.127 0.03826 -1.134 +-2.3 -2.1 6.69 12.39 6 15 3000 2.161 1.217 0.03826 -1.142 +-2.3 -2.1 12.39 18.09 6 15 3000 3.218 1.21 0.03662 -1.112 +-2.3 -2.1 18.09 23.79 6 15 3000 3.328 1.638 0.04398 -1.216 +-2.3 -2.1 23.79 29.49 6 15 3000 5.506 1.173 0.04403 -1.054 +-2.3 -2.1 29.49 35.19 6 15 3000 -2.444 3.613 0.05639 -1.437 +-2.3 -2.1 35.19 40.9 6 15 3000 2.217 3.133 0.05032 -1.338 +-2.1 -1.9 0 6.69 6 15 3000 1.184 0.8944 0.03233 -1.005 +-2.1 -1.9 6.69 12.39 6 15 3000 1.691 1.124 0.03736 -1.094 +-2.1 -1.9 12.39 18.09 6 15 3000 2.837 1.077 0.03437 -1.046 +-2.1 -1.9 18.09 23.79 6 15 3000 2.459 1.589 -0.04007 -1.18 +-2.1 -1.9 23.79 29.49 6 15 3000 4.058 1.369 -0.03922 -1.087 +-2.1 -1.9 29.49 35.19 6 15 3000 4.231 1.679 0.0432 -1.13 +-2.1 -1.9 35.19 40.9 6 15 3000 2.635 2.648 0.04929 -1.28 +-1.9 -1.7 0 6.69 6 15 3000 -0.8823 1.092 0.03599 -1.062 +-1.9 -1.7 6.69 12.39 6 15 3000 2.193 0.9891 0.03382 -1.012 +-1.9 -1.7 12.39 18.09 6 15 3000 2.9 1.043 0.03477 -1.019 +-1.9 -1.7 18.09 23.79 6 15 3000 2.371 1.488 -0.04053 -1.145 +-1.9 -1.7 23.79 29.49 6 15 3000 3.75 1.458 0.04346 -1.122 +-1.9 -1.7 29.49 35.19 6 15 3000 3.722 1.808 0.04668 -1.177 +-1.9 -1.7 35.19 40.9 6 15 3000 4.836 1.47 0.03875 -1.047 +-1.7 -1.3 0 6.69 6 15 3000 -1.692 1.192 0.05049 -1.06 +-1.7 -1.3 6.69 12.39 6 15 3000 -1.804 1.48 0.05315 -1.145 +-1.7 -1.3 12.39 18.09 6 15 3000 1.673 1.402 0.0536 -1.116 +-1.7 -1.3 18.09 23.79 6 15 3000 2.906 1.305 0.05377 -1.076 +-1.7 -1.3 23.79 29.49 6 15 3000 2.766 1.613 0.05511 -1.137 +-1.7 -1.3 29.49 35.19 6 15 3000 3.409 1.746 0.05585 -1.143 +-1.7 -1.3 35.19 40.9 6 15 3000 3.086 2.034 0.05795 -1.181 +-1.3 -1.1 0 6.69 6 15 3000 -0.7275 0.8099 0.04885 -0.9097 +-1.3 -1.1 6.69 12.39 6 15 3000 1.829 0.8156 0.04991 -0.9145 +-1.3 -1.1 12.39 18.09 6 15 3000 2.72 0.8454 0.05036 -0.9215 +-1.3 -1.1 18.09 23.79 6 15 3000 3.07 0.9201 0.05067 -0.9439 +-1.3 -1.1 23.79 29.49 6 15 3000 3.991 0.8715 0.05041 -0.9151 +-1.3 -1.1 29.49 35.19 6 15 3000 4.001 1.14 0.05214 -0.9987 +-1.3 -1.1 35.19 40.9 6 15 3000 4.522 1.22 0.05122 -1 +-1.1 -0.8 0 6.69 6 15 3000 1.423 0.4736 0.03233 -0.7093 +-1.1 -0.8 6.69 12.39 6 15 3000 2.249 0.5041 0.03355 -0.7316 +-1.1 -0.8 12.39 18.09 6 15 3000 2.961 0.4889 0.03129 -0.7091 +-1.1 -0.8 18.09 23.79 6 15 3000 3.4 0.5757 0.03541 -0.7742 +-1.1 -0.8 23.79 29.49 6 15 3000 3.884 0.6457 0.03731 -0.8146 +-1.1 -0.8 29.49 35.19 6 15 3000 4.433 0.7524 0.03962 -0.8672 +-1.1 -0.8 35.19 40.9 6 15 3000 4.681 0.9075 0.04182 -0.9304 +-0.8 -0.5 0 6.69 6 15 3000 1.003 0.4142 0.02486 -0.6698 +-0.8 -0.5 6.69 12.39 6 15 3000 2.134 0.3971 0.02264 -0.6469 +-0.8 -0.5 12.39 18.09 6 15 3000 2.66 0.4566 0.02755 -0.7058 +-0.8 -0.5 18.09 23.79 6 15 3000 3.264 0.4799 0.02702 -0.7156 +-0.8 -0.5 23.79 29.49 6 15 3000 3.877 0.5249 0.02923 -0.7479 +-0.8 -0.5 29.49 35.19 6 15 3000 4.441 0.581 0.03045 -0.7804 +-0.8 -0.5 35.19 40.9 6 15 3000 4.742 0.8003 0.03613 -0.9062 +-0.5 0 0 6.69 6 15 3000 0.6172 0.3908 0.02003 -0.6407 +-0.5 0 6.69 12.39 6 15 3000 1.775 0.4231 0.02199 -0.6701 +-0.5 0 12.39 18.09 6 15 3000 2.457 0.4626 0.02416 -0.7045 +-0.5 0 18.09 23.79 6 15 3000 2.996 0.5242 0.02689 -0.7508 +-0.5 0 23.79 29.49 6 15 3000 3.623 0.5591 0.0288 -0.7747 +-0.5 0 29.49 35.19 6 15 3000 4.167 0.6365 0.03045 -0.8179 +-0.5 0 35.19 40.9 6 15 3000 4.795 0.6819 0.03145 -0.8408 +0 0.5 0 6.69 6 15 3000 0.6172 0.3908 0.02003 -0.6407 +0 0.5 6.69 12.39 6 15 3000 1.775 0.4231 0.02199 -0.6701 +0 0.5 12.39 18.09 6 15 3000 2.457 0.4626 0.02416 -0.7045 +0 0.5 18.09 23.79 6 15 3000 2.996 0.5242 0.02689 -0.7508 +0 0.5 23.79 29.49 6 15 3000 3.623 0.5591 0.0288 -0.7747 +0 0.5 29.49 35.19 6 15 3000 4.167 0.6365 0.03045 -0.8179 +0 0.5 35.19 40.9 6 15 3000 4.795 0.6819 0.03145 -0.8408 +0.5 0.8 0 6.69 6 15 3000 1.003 0.4142 0.02486 -0.6698 +0.5 0.8 6.69 12.39 6 15 3000 2.134 0.3971 0.02264 -0.6469 +0.5 0.8 12.39 18.09 6 15 3000 2.66 0.4566 0.02755 -0.7058 +0.5 0.8 18.09 23.79 6 15 3000 3.264 0.4799 0.02702 -0.7156 +0.5 0.8 23.79 29.49 6 15 3000 3.877 0.5249 0.02923 -0.7479 +0.5 0.8 29.49 35.19 6 15 3000 4.441 0.581 0.03045 -0.7804 +0.5 0.8 35.19 40.9 6 15 3000 4.742 0.8003 0.03613 -0.9062 +0.8 1.1 0 6.69 6 15 3000 1.423 0.4736 0.03233 -0.7093 +0.8 1.1 6.69 12.39 6 15 3000 2.249 0.5041 0.03355 -0.7316 +0.8 1.1 12.39 18.09 6 15 3000 2.961 0.4889 0.03129 -0.7091 +0.8 1.1 18.09 23.79 6 15 3000 3.4 0.5757 0.03541 -0.7742 +0.8 1.1 23.79 29.49 6 15 3000 3.884 0.6457 0.03731 -0.8146 +0.8 1.1 29.49 35.19 6 15 3000 4.433 0.7524 0.03962 -0.8672 +0.8 1.1 35.19 40.9 6 15 3000 4.681 0.9075 0.04182 -0.9304 +1.1 1.3 0 6.69 6 15 3000 -0.7275 0.8099 0.04885 -0.9097 +1.1 1.3 6.69 12.39 6 15 3000 1.829 0.8156 0.04991 -0.9145 +1.1 1.3 12.39 18.09 6 15 3000 2.72 0.8454 0.05036 -0.9215 +1.1 1.3 18.09 23.79 6 15 3000 3.07 0.9201 0.05067 -0.9439 +1.1 1.3 23.79 29.49 6 15 3000 3.991 0.8715 0.05041 -0.9151 +1.1 1.3 29.49 35.19 6 15 3000 4.001 1.14 0.05214 -0.9987 +1.1 1.3 35.19 40.9 6 15 3000 4.522 1.22 0.05122 -1 +1.3 1.7 0 6.69 6 15 3000 -1.692 1.192 0.05049 -1.06 +1.3 1.7 6.69 12.39 6 15 3000 -1.804 1.48 0.05315 -1.145 +1.3 1.7 12.39 18.09 6 15 3000 1.673 1.402 0.0536 -1.116 +1.3 1.7 18.09 23.79 6 15 3000 2.906 1.305 0.05377 -1.076 +1.3 1.7 23.79 29.49 6 15 3000 2.766 1.613 0.05511 -1.137 +1.3 1.7 29.49 35.19 6 15 3000 3.409 1.746 0.05585 -1.143 +1.3 1.7 35.19 40.9 6 15 3000 3.086 2.034 0.05795 -1.181 +1.7 1.9 0 6.69 6 15 3000 -0.8823 1.092 0.03599 -1.062 +1.7 1.9 6.69 12.39 6 15 3000 2.193 0.9891 0.03382 -1.012 +1.7 1.9 12.39 18.09 6 15 3000 2.9 1.043 0.03477 -1.019 +1.7 1.9 18.09 23.79 6 15 3000 2.371 1.488 -0.04053 -1.145 +1.7 1.9 23.79 29.49 6 15 3000 3.75 1.458 0.04346 -1.122 +1.7 1.9 29.49 35.19 6 15 3000 3.722 1.808 0.04668 -1.177 +1.7 1.9 35.19 40.9 6 15 3000 4.836 1.47 0.03875 -1.047 +1.9 2.1 0 6.69 6 15 3000 1.184 0.8944 0.03233 -1.005 +1.9 2.1 6.69 12.39 6 15 3000 1.691 1.124 0.03736 -1.094 +1.9 2.1 12.39 18.09 6 15 3000 2.837 1.077 0.03437 -1.046 +1.9 2.1 18.09 23.79 6 15 3000 2.459 1.589 -0.04007 -1.18 +1.9 2.1 23.79 29.49 6 15 3000 4.058 1.369 -0.03922 -1.087 +1.9 2.1 29.49 35.19 6 15 3000 4.231 1.679 0.0432 -1.13 +1.9 2.1 35.19 40.9 6 15 3000 2.635 2.648 0.04929 -1.28 +2.1 2.3 0 6.69 6 15 3000 0.3022 1.127 0.03826 -1.134 +2.1 2.3 6.69 12.39 6 15 3000 2.161 1.217 0.03826 -1.142 +2.1 2.3 12.39 18.09 6 15 3000 3.218 1.21 0.03662 -1.112 +2.1 2.3 18.09 23.79 6 15 3000 3.328 1.638 0.04398 -1.216 +2.1 2.3 23.79 29.49 6 15 3000 5.506 1.173 0.04403 -1.054 +2.1 2.3 29.49 35.19 6 15 3000 -2.444 3.613 0.05639 -1.437 +2.1 2.3 35.19 40.9 6 15 3000 2.217 3.133 0.05032 -1.338 +2.3 2.5 0 6.69 6 15 3000 3.125 0.6026 0.02576 -0.8702 +2.3 2.5 6.69 12.39 6 15 3000 3.935 0.6533 0.02587 -0.889 +2.3 2.5 12.39 18.09 6 15 3000 4.198 1.024 0.03618 -1.069 +2.3 2.5 18.09 23.79 6 15 3000 2.948 2.386 0.04771 -1.382 +2.3 2.5 23.79 29.49 6 15 3000 4.415 2.086 0.04704 -1.294 +2.3 2.5 29.49 35.19 6 15 3000 -3.084 4.156 0.05366 -1.503 +2.3 2.5 35.19 40.9 6 15 3000 -6.144 5.969 0.05633 -1.602 +2.5 2.8 0 6.69 6 15 3000 4.244 0.2766 -1.86e-08 -0.5068 +2.5 2.8 6.69 12.39 6 15 3000 4.919 0.3193 5.463e-06 -0.58 +2.5 2.8 12.39 18.09 6 15 3000 5.909 0.2752 4.144e-06 -0.5272 +2.5 2.8 18.09 23.79 6 15 3000 -47.31 47.18 0.05853 -1.991 +2.5 2.8 23.79 29.49 6 15 3000 -46.49 46.33 0.05698 -1.989 +2.5 2.8 29.49 35.19 6 15 3000 8.651 0.2522 6.592e-06 -0.4835 +2.5 2.8 35.19 40.9 6 15 3000 7.716 2.481 0.0531 -1.455 +2.8 3 0 6.69 6 15 3000 4.467 0.1997 -3.491e-06 -0.2623 +2.8 3 6.69 12.39 6 15 3000 4.17 0.928 0.07702 -1.063 +2.8 3 12.39 18.09 6 15 3000 -0.04491 3.67 0.08704 -1.641 +2.8 3 18.09 23.79 6 15 3000 5.528 1.286 0.07962 -1.187 +2.8 3 23.79 29.49 6 15 3000 -78.36 78.23 0.08448 -1.996 +2.8 3 29.49 35.19 6 15 3000 7.559 1.147 0.07023 -1.134 +2.8 3 35.19 40.9 6 15 3000 -59.03 59.03 -0.08184 -1.992 +3 3.2 0 6.69 6 15 3000 0.0002851 3.01 0.1382 -1.702 +3 3.2 6.69 12.39 6 15 3000 -33.01 33.04 0.1343 -1.991 +3 3.2 12.39 18.09 6 15 3000 -67.94 67.8 0.1342 -1.996 +3 3.2 18.09 23.79 6 15 3000 -47.81 48 0.1391 -1.996 +3 3.2 23.79 29.49 6 15 3000 7.162 0.9211 0.1395 -1.209 +3 3.2 29.49 35.19 6 15 3000 8.193 0.1995 2.822e-05 -0.132 +3 3.2 35.19 40.9 6 15 3000 8.133 0.9983 0.1349 -1.181 +3.2 4.7 0 6.69 6 15 3000 2.511 0.3167 0.09085 -0.7407 +3.2 4.7 6.69 12.39 6 15 3000 3.297 0.2091 6.258e-05 -0.2755 +3.2 4.7 12.39 18.09 6 15 3000 1.85 2.281 0.1042 -1.635 +3.2 4.7 18.09 23.79 6 15 3000 3.869 1.001 0.09955 -1.266 +3.2 4.7 23.79 29.49 6 15 3000 -23.98 24.11 0.1057 -1.988 +3.2 4.7 29.49 35.19 6 15 3000 5.403 0.2371 1.5e-05 -0.3177 +3.2 4.7 35.19 40.9 6 15 3000 5.753 0.2337 0.0002982 -0.3108 diff --git a/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt b/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt new file mode 100644 index 00000000..4c080015 --- /dev/null +++ b/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt @@ -0,0 +1,27 @@ +{1 JetEta 0 None ScaleFactor} +-5.191 -3.139 3 1.1922 1.0434 1.341 +-3.139 -2.964 3 1.1869 1.0626 1.3112 +-2.964 -2.853 3 1.7788 1.578 1.9796 +-2.853 -2.500 3 1.3418 1.1327 1.5509 +-2.500 -2.322 3 1.2963 1.0592 1.5334 +-2.322 -2.043 3 1.1512 1.0372 1.2652 +-2.043 -1.930 3 1.1426 1.0212 1.264 +-1.930 -1.740 3 1.1000 0.9921 1.2079 +-1.740 -1.305 3 1.1278 1.0292 1.2264 +-1.305 -1.131 3 1.1609 1.0584 1.2634 +-1.131 -0.783 3 1.1464 1.0832 1.2096 +-0.783 -0.522 3 1.1948 1.1296 1.26 +-0.522 0.000 3 1.1595 1.095 1.224 +0.000 0.522 3 1.1595 1.095 1.224 +0.522 0.783 3 1.1948 1.1296 1.26 +0.783 1.131 3 1.1464 1.0832 1.2096 +1.131 1.305 3 1.1609 1.0584 1.2634 +1.305 1.740 3 1.1278 1.0292 1.2264 +1.740 1.930 3 1.1000 0.9921 1.2079 +1.930 2.043 3 1.1426 1.0212 1.264 +2.043 2.322 3 1.1512 1.0372 1.2652 +2.322 2.500 3 1.2963 1.0592 1.5334 +2.500 2.853 3 1.3418 1.1327 1.5509 +2.853 2.964 3 1.7788 1.578 1.9796 +2.964 3.139 3 1.1869 1.0626 1.3112 +3.139 5.191 3 1.1922 1.0434 1.341 From 535dd8b471f4b17934e2cdd46a7275cd2ef07326 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 19 Dec 2018 03:20:37 +0900 Subject: [PATCH 19/23] Fixed the source paths --- external/src/JetResolution.cc | 2 +- external/src/JetResolutionObject.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/src/JetResolution.cc b/external/src/JetResolution.cc index af375238..499202e4 100644 --- a/external/src/JetResolution.cc +++ b/external/src/JetResolution.cc @@ -1,6 +1,6 @@ // // This code is from copying -// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/JetResolution.cc +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/JetMETCorrections/Modules/src/JetResolution.cc // in 4d7b6ef commit // Copied by Byeonghak Ko // diff --git a/external/src/JetResolutionObject.cc b/external/src/JetResolutionObject.cc index 373251f7..9daeab4f 100644 --- a/external/src/JetResolutionObject.cc +++ b/external/src/JetResolutionObject.cc @@ -1,6 +1,6 @@ // // This code is from copying -// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/interface/JetResolutionObject.h +// https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/CondFormats/JetMETObjects/src/JetResolutionObject.cc // in a2f9cbb commit // Copied by Byeonghak Ko // From 04afb129c070842762156e48f517c35c61e09db4 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Wed, 19 Dec 2018 13:40:18 +0900 Subject: [PATCH 20/23] Changed some name for avoiding conflicts --- analysis/interface/topObjectSelection.h | 4 ++-- analysis/src/topObjectSelection.cc | 10 +++++----- external/interface/JetCorrectorParametersHelper.h | 4 ++-- external/interface/JetResolution.h | 2 +- external/interface/JetResolutionObject.h | 8 ++++---- external/src/JetResolution.cc | 2 +- external/src/JetResolutionObject.cc | 4 ++-- 7 files changed, 17 insertions(+), 17 deletions(-) diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 9063aed2..36b94712 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -23,8 +23,8 @@ class topObjectSelection : public nanoBase JetCorrectionUncertainty *jecUnc; - JME::JetResolution jetResObj; - JME::JetResolutionScaleFactor jetResSFObj; + JMENano::JetResolution jetResObj; + JMENano::JetResolutionScaleFactor jetResSFObj; TRandom3 *rndEngine; public: diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index ebd8e7b6..78fc93cd 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -11,11 +11,11 @@ topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jer/" "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; - jetResSFObj = JME::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); + jetResSFObj = JMENano::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); std::string strPathJetResObj = env + "/src/nano/analysis/data/jer/" "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; - jetResObj = JME::JetResolution(strPathJetResObj.c_str()); + jetResObj = JMENano::JetResolution(strPathJetResObj.c_str()); rndEngine = new TRandom3(12345); @@ -261,9 +261,9 @@ void topObjectSelection::GetJetMassPt(UInt_t nIdx, //if ( m_isMC && ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) if ( m_isMC ) { // Evaluating the central part cJER of the factor - JME::JetParameters jetPars = {{JME::Binning::JetPt, fJetPt}, - {JME::Binning::JetEta, fJetEta}, - {JME::Binning::Rho, fixedGridRhoFastjetAll}}; + JMENano::JetParameters jetPars = {{JMENano::Binning::JetPt, fJetPt}, + {JMENano::Binning::JetEta, fJetEta}, + {JMENano::Binning::Rho, fixedGridRhoFastjetAll}}; const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. const double cJER = jetResSFObj.getScaleFactor(jetPars, diff --git a/external/interface/JetCorrectorParametersHelper.h b/external/interface/JetCorrectorParametersHelper.h index 136e3ff3..8edf1008 100644 --- a/external/interface/JetCorrectorParametersHelper.h +++ b/external/interface/JetCorrectorParametersHelper.h @@ -27,8 +27,8 @@ #include #include #include // MODIFIED -//#include "FWCore/Utilities/interface/Exception.h" -//#include "FWCore/MessageLogger/interface/MessageLogger.h" +//#include "FWCore/Utilities/interface/Exception.h" // MODIFIED +//#include "FWCore/MessageLogger/interface/MessageLogger.h" // MODIFIED //---------- JetCorrectorParametersHelper class ---------------- diff --git a/external/interface/JetResolution.h b/external/interface/JetResolution.h index 7ad735b8..69ee2703 100644 --- a/external/interface/JetResolution.h +++ b/external/interface/JetResolution.h @@ -27,7 +27,7 @@ namespace edm { #endif -namespace JME { +namespace JMENano { class JetResolution { public: JetResolution(const std::string& filename); diff --git a/external/interface/JetResolutionObject.h b/external/interface/JetResolutionObject.h index aad1f3c3..29d35275 100644 --- a/external/interface/JetResolutionObject.h +++ b/external/interface/JetResolutionObject.h @@ -47,7 +47,7 @@ T clip(const T& n, const T& lower, const T& upper) { return std::max(lower, std::min(n, upper)); } -namespace JME { +namespace JMENano { template struct bimap { typedef std::unordered_map left_type; @@ -89,8 +89,8 @@ namespace JME { // Hash function for Binning enum class namespace std { template<> - struct hash { - typedef JME::Binning argument_type; + struct hash { + typedef JMENano::Binning argument_type; typedef std::size_t result_type; hash int_hash; @@ -101,7 +101,7 @@ namespace std { }; }; -namespace JME { +namespace JMENano { class JetParameters { public: diff --git a/external/src/JetResolution.cc b/external/src/JetResolution.cc index 499202e4..0dfbc6c4 100644 --- a/external/src/JetResolution.cc +++ b/external/src/JetResolution.cc @@ -21,7 +21,7 @@ #include "nano/external/interface/JetResolution.h" // MODIFIED #endif -namespace JME { +namespace JMENano { JetResolution::JetResolution(const std::string& filename) { m_object = std::make_shared(filename); diff --git a/external/src/JetResolutionObject.cc b/external/src/JetResolutionObject.cc index 9daeab4f..b65804a4 100644 --- a/external/src/JetResolutionObject.cc +++ b/external/src/JetResolutionObject.cc @@ -38,7 +38,7 @@ namespace edm { #include #include -namespace JME { +namespace JMENano { std::string getDefinitionLine(const std::string& line) { size_t first = line.find ('{'); @@ -137,7 +137,7 @@ namespace JME { if (it == m_values.cend()) { throwException(edm::errors::NotFound, "JER parametrisation depends on '" + JetParameters::binning_to_string.left.at(bin) + - "' but no value for this parameter has been specified. Please call the appropriate 'set' function of the JME::JetParameters object"); + "' but no value for this parameter has been specified. Please call the appropriate 'set' function of the JMENano::JetParameters object"); } values.push_back(it->second); From e010f7e34c48a307224d483a15419d8c26b0e346 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Thu, 20 Dec 2018 19:41:01 +0900 Subject: [PATCH 21/23] Modifying getFiles instead of including the JER/JES text files into the repositories (neither nano nor nanoData) --- ...mmer16_25nsV1_MC_PtResolution_AK4PFchs.txt | 183 ------------------ .../jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt | 27 --- nanoAOD/scripts/getFiles | 5 + 3 files changed, 5 insertions(+), 210 deletions(-) delete mode 100644 analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt delete mode 100644 analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt diff --git a/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt b/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt deleted file mode 100644 index 7633031e..00000000 --- a/analysis/data/jer/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt +++ /dev/null @@ -1,183 +0,0 @@ -{2 JetEta Rho 1 JetPt sqrt([0]*abs([0])/(x*x)+[1]*[1]*pow(x,[3])+[2]*[2]) Resolution} --4.7 -3.2 0 6.69 6 15 3000 2.511 0.3167 0.09085 -0.7407 --4.7 -3.2 6.69 12.39 6 15 3000 3.297 0.2091 6.258e-05 -0.2755 --4.7 -3.2 12.39 18.09 6 15 3000 1.85 2.281 0.1042 -1.635 --4.7 -3.2 18.09 23.79 6 15 3000 3.869 1.001 0.09955 -1.266 --4.7 -3.2 23.79 29.49 6 15 3000 -23.98 24.11 0.1057 -1.988 --4.7 -3.2 29.49 35.19 6 15 3000 5.403 0.2371 1.5e-05 -0.3177 --4.7 -3.2 35.19 40.9 6 15 3000 5.753 0.2337 0.0002982 -0.3108 --3.2 -3 0 6.69 6 15 3000 0.0002851 3.01 0.1382 -1.702 --3.2 -3 6.69 12.39 6 15 3000 -33.01 33.04 0.1343 -1.991 --3.2 -3 12.39 18.09 6 15 3000 -67.94 67.8 0.1342 -1.996 --3.2 -3 18.09 23.79 6 15 3000 -47.81 48 0.1391 -1.996 --3.2 -3 23.79 29.49 6 15 3000 7.162 0.9211 0.1395 -1.209 --3.2 -3 29.49 35.19 6 15 3000 8.193 0.1995 2.822e-05 -0.132 --3.2 -3 35.19 40.9 6 15 3000 8.133 0.9983 0.1349 -1.181 --3 -2.8 0 6.69 6 15 3000 4.467 0.1997 -3.491e-06 -0.2623 --3 -2.8 6.69 12.39 6 15 3000 4.17 0.928 0.07702 -1.063 --3 -2.8 12.39 18.09 6 15 3000 -0.04491 3.67 0.08704 -1.641 --3 -2.8 18.09 23.79 6 15 3000 5.528 1.286 0.07962 -1.187 --3 -2.8 23.79 29.49 6 15 3000 -78.36 78.23 0.08448 -1.996 --3 -2.8 29.49 35.19 6 15 3000 7.559 1.147 0.07023 -1.134 --3 -2.8 35.19 40.9 6 15 3000 -59.03 59.03 -0.08184 -1.992 --2.8 -2.5 0 6.69 6 15 3000 4.244 0.2766 -1.86e-08 -0.5068 --2.8 -2.5 6.69 12.39 6 15 3000 4.919 0.3193 5.463e-06 -0.58 --2.8 -2.5 12.39 18.09 6 15 3000 5.909 0.2752 4.144e-06 -0.5272 --2.8 -2.5 18.09 23.79 6 15 3000 -47.31 47.18 0.05853 -1.991 --2.8 -2.5 23.79 29.49 6 15 3000 -46.49 46.33 0.05698 -1.989 --2.8 -2.5 29.49 35.19 6 15 3000 8.651 0.2522 6.592e-06 -0.4835 --2.8 -2.5 35.19 40.9 6 15 3000 7.716 2.481 0.0531 -1.455 --2.5 -2.3 0 6.69 6 15 3000 3.125 0.6026 0.02576 -0.8702 --2.5 -2.3 6.69 12.39 6 15 3000 3.935 0.6533 0.02587 -0.889 --2.5 -2.3 12.39 18.09 6 15 3000 4.198 1.024 0.03618 -1.069 --2.5 -2.3 18.09 23.79 6 15 3000 2.948 2.386 0.04771 -1.382 --2.5 -2.3 23.79 29.49 6 15 3000 4.415 2.086 0.04704 -1.294 --2.5 -2.3 29.49 35.19 6 15 3000 -3.084 4.156 0.05366 -1.503 --2.5 -2.3 35.19 40.9 6 15 3000 -6.144 5.969 0.05633 -1.602 --2.3 -2.1 0 6.69 6 15 3000 0.3022 1.127 0.03826 -1.134 --2.3 -2.1 6.69 12.39 6 15 3000 2.161 1.217 0.03826 -1.142 --2.3 -2.1 12.39 18.09 6 15 3000 3.218 1.21 0.03662 -1.112 --2.3 -2.1 18.09 23.79 6 15 3000 3.328 1.638 0.04398 -1.216 --2.3 -2.1 23.79 29.49 6 15 3000 5.506 1.173 0.04403 -1.054 --2.3 -2.1 29.49 35.19 6 15 3000 -2.444 3.613 0.05639 -1.437 --2.3 -2.1 35.19 40.9 6 15 3000 2.217 3.133 0.05032 -1.338 --2.1 -1.9 0 6.69 6 15 3000 1.184 0.8944 0.03233 -1.005 --2.1 -1.9 6.69 12.39 6 15 3000 1.691 1.124 0.03736 -1.094 --2.1 -1.9 12.39 18.09 6 15 3000 2.837 1.077 0.03437 -1.046 --2.1 -1.9 18.09 23.79 6 15 3000 2.459 1.589 -0.04007 -1.18 --2.1 -1.9 23.79 29.49 6 15 3000 4.058 1.369 -0.03922 -1.087 --2.1 -1.9 29.49 35.19 6 15 3000 4.231 1.679 0.0432 -1.13 --2.1 -1.9 35.19 40.9 6 15 3000 2.635 2.648 0.04929 -1.28 --1.9 -1.7 0 6.69 6 15 3000 -0.8823 1.092 0.03599 -1.062 --1.9 -1.7 6.69 12.39 6 15 3000 2.193 0.9891 0.03382 -1.012 --1.9 -1.7 12.39 18.09 6 15 3000 2.9 1.043 0.03477 -1.019 --1.9 -1.7 18.09 23.79 6 15 3000 2.371 1.488 -0.04053 -1.145 --1.9 -1.7 23.79 29.49 6 15 3000 3.75 1.458 0.04346 -1.122 --1.9 -1.7 29.49 35.19 6 15 3000 3.722 1.808 0.04668 -1.177 --1.9 -1.7 35.19 40.9 6 15 3000 4.836 1.47 0.03875 -1.047 --1.7 -1.3 0 6.69 6 15 3000 -1.692 1.192 0.05049 -1.06 --1.7 -1.3 6.69 12.39 6 15 3000 -1.804 1.48 0.05315 -1.145 --1.7 -1.3 12.39 18.09 6 15 3000 1.673 1.402 0.0536 -1.116 --1.7 -1.3 18.09 23.79 6 15 3000 2.906 1.305 0.05377 -1.076 --1.7 -1.3 23.79 29.49 6 15 3000 2.766 1.613 0.05511 -1.137 --1.7 -1.3 29.49 35.19 6 15 3000 3.409 1.746 0.05585 -1.143 --1.7 -1.3 35.19 40.9 6 15 3000 3.086 2.034 0.05795 -1.181 --1.3 -1.1 0 6.69 6 15 3000 -0.7275 0.8099 0.04885 -0.9097 --1.3 -1.1 6.69 12.39 6 15 3000 1.829 0.8156 0.04991 -0.9145 --1.3 -1.1 12.39 18.09 6 15 3000 2.72 0.8454 0.05036 -0.9215 --1.3 -1.1 18.09 23.79 6 15 3000 3.07 0.9201 0.05067 -0.9439 --1.3 -1.1 23.79 29.49 6 15 3000 3.991 0.8715 0.05041 -0.9151 --1.3 -1.1 29.49 35.19 6 15 3000 4.001 1.14 0.05214 -0.9987 --1.3 -1.1 35.19 40.9 6 15 3000 4.522 1.22 0.05122 -1 --1.1 -0.8 0 6.69 6 15 3000 1.423 0.4736 0.03233 -0.7093 --1.1 -0.8 6.69 12.39 6 15 3000 2.249 0.5041 0.03355 -0.7316 --1.1 -0.8 12.39 18.09 6 15 3000 2.961 0.4889 0.03129 -0.7091 --1.1 -0.8 18.09 23.79 6 15 3000 3.4 0.5757 0.03541 -0.7742 --1.1 -0.8 23.79 29.49 6 15 3000 3.884 0.6457 0.03731 -0.8146 --1.1 -0.8 29.49 35.19 6 15 3000 4.433 0.7524 0.03962 -0.8672 --1.1 -0.8 35.19 40.9 6 15 3000 4.681 0.9075 0.04182 -0.9304 --0.8 -0.5 0 6.69 6 15 3000 1.003 0.4142 0.02486 -0.6698 --0.8 -0.5 6.69 12.39 6 15 3000 2.134 0.3971 0.02264 -0.6469 --0.8 -0.5 12.39 18.09 6 15 3000 2.66 0.4566 0.02755 -0.7058 --0.8 -0.5 18.09 23.79 6 15 3000 3.264 0.4799 0.02702 -0.7156 --0.8 -0.5 23.79 29.49 6 15 3000 3.877 0.5249 0.02923 -0.7479 --0.8 -0.5 29.49 35.19 6 15 3000 4.441 0.581 0.03045 -0.7804 --0.8 -0.5 35.19 40.9 6 15 3000 4.742 0.8003 0.03613 -0.9062 --0.5 0 0 6.69 6 15 3000 0.6172 0.3908 0.02003 -0.6407 --0.5 0 6.69 12.39 6 15 3000 1.775 0.4231 0.02199 -0.6701 --0.5 0 12.39 18.09 6 15 3000 2.457 0.4626 0.02416 -0.7045 --0.5 0 18.09 23.79 6 15 3000 2.996 0.5242 0.02689 -0.7508 --0.5 0 23.79 29.49 6 15 3000 3.623 0.5591 0.0288 -0.7747 --0.5 0 29.49 35.19 6 15 3000 4.167 0.6365 0.03045 -0.8179 --0.5 0 35.19 40.9 6 15 3000 4.795 0.6819 0.03145 -0.8408 -0 0.5 0 6.69 6 15 3000 0.6172 0.3908 0.02003 -0.6407 -0 0.5 6.69 12.39 6 15 3000 1.775 0.4231 0.02199 -0.6701 -0 0.5 12.39 18.09 6 15 3000 2.457 0.4626 0.02416 -0.7045 -0 0.5 18.09 23.79 6 15 3000 2.996 0.5242 0.02689 -0.7508 -0 0.5 23.79 29.49 6 15 3000 3.623 0.5591 0.0288 -0.7747 -0 0.5 29.49 35.19 6 15 3000 4.167 0.6365 0.03045 -0.8179 -0 0.5 35.19 40.9 6 15 3000 4.795 0.6819 0.03145 -0.8408 -0.5 0.8 0 6.69 6 15 3000 1.003 0.4142 0.02486 -0.6698 -0.5 0.8 6.69 12.39 6 15 3000 2.134 0.3971 0.02264 -0.6469 -0.5 0.8 12.39 18.09 6 15 3000 2.66 0.4566 0.02755 -0.7058 -0.5 0.8 18.09 23.79 6 15 3000 3.264 0.4799 0.02702 -0.7156 -0.5 0.8 23.79 29.49 6 15 3000 3.877 0.5249 0.02923 -0.7479 -0.5 0.8 29.49 35.19 6 15 3000 4.441 0.581 0.03045 -0.7804 -0.5 0.8 35.19 40.9 6 15 3000 4.742 0.8003 0.03613 -0.9062 -0.8 1.1 0 6.69 6 15 3000 1.423 0.4736 0.03233 -0.7093 -0.8 1.1 6.69 12.39 6 15 3000 2.249 0.5041 0.03355 -0.7316 -0.8 1.1 12.39 18.09 6 15 3000 2.961 0.4889 0.03129 -0.7091 -0.8 1.1 18.09 23.79 6 15 3000 3.4 0.5757 0.03541 -0.7742 -0.8 1.1 23.79 29.49 6 15 3000 3.884 0.6457 0.03731 -0.8146 -0.8 1.1 29.49 35.19 6 15 3000 4.433 0.7524 0.03962 -0.8672 -0.8 1.1 35.19 40.9 6 15 3000 4.681 0.9075 0.04182 -0.9304 -1.1 1.3 0 6.69 6 15 3000 -0.7275 0.8099 0.04885 -0.9097 -1.1 1.3 6.69 12.39 6 15 3000 1.829 0.8156 0.04991 -0.9145 -1.1 1.3 12.39 18.09 6 15 3000 2.72 0.8454 0.05036 -0.9215 -1.1 1.3 18.09 23.79 6 15 3000 3.07 0.9201 0.05067 -0.9439 -1.1 1.3 23.79 29.49 6 15 3000 3.991 0.8715 0.05041 -0.9151 -1.1 1.3 29.49 35.19 6 15 3000 4.001 1.14 0.05214 -0.9987 -1.1 1.3 35.19 40.9 6 15 3000 4.522 1.22 0.05122 -1 -1.3 1.7 0 6.69 6 15 3000 -1.692 1.192 0.05049 -1.06 -1.3 1.7 6.69 12.39 6 15 3000 -1.804 1.48 0.05315 -1.145 -1.3 1.7 12.39 18.09 6 15 3000 1.673 1.402 0.0536 -1.116 -1.3 1.7 18.09 23.79 6 15 3000 2.906 1.305 0.05377 -1.076 -1.3 1.7 23.79 29.49 6 15 3000 2.766 1.613 0.05511 -1.137 -1.3 1.7 29.49 35.19 6 15 3000 3.409 1.746 0.05585 -1.143 -1.3 1.7 35.19 40.9 6 15 3000 3.086 2.034 0.05795 -1.181 -1.7 1.9 0 6.69 6 15 3000 -0.8823 1.092 0.03599 -1.062 -1.7 1.9 6.69 12.39 6 15 3000 2.193 0.9891 0.03382 -1.012 -1.7 1.9 12.39 18.09 6 15 3000 2.9 1.043 0.03477 -1.019 -1.7 1.9 18.09 23.79 6 15 3000 2.371 1.488 -0.04053 -1.145 -1.7 1.9 23.79 29.49 6 15 3000 3.75 1.458 0.04346 -1.122 -1.7 1.9 29.49 35.19 6 15 3000 3.722 1.808 0.04668 -1.177 -1.7 1.9 35.19 40.9 6 15 3000 4.836 1.47 0.03875 -1.047 -1.9 2.1 0 6.69 6 15 3000 1.184 0.8944 0.03233 -1.005 -1.9 2.1 6.69 12.39 6 15 3000 1.691 1.124 0.03736 -1.094 -1.9 2.1 12.39 18.09 6 15 3000 2.837 1.077 0.03437 -1.046 -1.9 2.1 18.09 23.79 6 15 3000 2.459 1.589 -0.04007 -1.18 -1.9 2.1 23.79 29.49 6 15 3000 4.058 1.369 -0.03922 -1.087 -1.9 2.1 29.49 35.19 6 15 3000 4.231 1.679 0.0432 -1.13 -1.9 2.1 35.19 40.9 6 15 3000 2.635 2.648 0.04929 -1.28 -2.1 2.3 0 6.69 6 15 3000 0.3022 1.127 0.03826 -1.134 -2.1 2.3 6.69 12.39 6 15 3000 2.161 1.217 0.03826 -1.142 -2.1 2.3 12.39 18.09 6 15 3000 3.218 1.21 0.03662 -1.112 -2.1 2.3 18.09 23.79 6 15 3000 3.328 1.638 0.04398 -1.216 -2.1 2.3 23.79 29.49 6 15 3000 5.506 1.173 0.04403 -1.054 -2.1 2.3 29.49 35.19 6 15 3000 -2.444 3.613 0.05639 -1.437 -2.1 2.3 35.19 40.9 6 15 3000 2.217 3.133 0.05032 -1.338 -2.3 2.5 0 6.69 6 15 3000 3.125 0.6026 0.02576 -0.8702 -2.3 2.5 6.69 12.39 6 15 3000 3.935 0.6533 0.02587 -0.889 -2.3 2.5 12.39 18.09 6 15 3000 4.198 1.024 0.03618 -1.069 -2.3 2.5 18.09 23.79 6 15 3000 2.948 2.386 0.04771 -1.382 -2.3 2.5 23.79 29.49 6 15 3000 4.415 2.086 0.04704 -1.294 -2.3 2.5 29.49 35.19 6 15 3000 -3.084 4.156 0.05366 -1.503 -2.3 2.5 35.19 40.9 6 15 3000 -6.144 5.969 0.05633 -1.602 -2.5 2.8 0 6.69 6 15 3000 4.244 0.2766 -1.86e-08 -0.5068 -2.5 2.8 6.69 12.39 6 15 3000 4.919 0.3193 5.463e-06 -0.58 -2.5 2.8 12.39 18.09 6 15 3000 5.909 0.2752 4.144e-06 -0.5272 -2.5 2.8 18.09 23.79 6 15 3000 -47.31 47.18 0.05853 -1.991 -2.5 2.8 23.79 29.49 6 15 3000 -46.49 46.33 0.05698 -1.989 -2.5 2.8 29.49 35.19 6 15 3000 8.651 0.2522 6.592e-06 -0.4835 -2.5 2.8 35.19 40.9 6 15 3000 7.716 2.481 0.0531 -1.455 -2.8 3 0 6.69 6 15 3000 4.467 0.1997 -3.491e-06 -0.2623 -2.8 3 6.69 12.39 6 15 3000 4.17 0.928 0.07702 -1.063 -2.8 3 12.39 18.09 6 15 3000 -0.04491 3.67 0.08704 -1.641 -2.8 3 18.09 23.79 6 15 3000 5.528 1.286 0.07962 -1.187 -2.8 3 23.79 29.49 6 15 3000 -78.36 78.23 0.08448 -1.996 -2.8 3 29.49 35.19 6 15 3000 7.559 1.147 0.07023 -1.134 -2.8 3 35.19 40.9 6 15 3000 -59.03 59.03 -0.08184 -1.992 -3 3.2 0 6.69 6 15 3000 0.0002851 3.01 0.1382 -1.702 -3 3.2 6.69 12.39 6 15 3000 -33.01 33.04 0.1343 -1.991 -3 3.2 12.39 18.09 6 15 3000 -67.94 67.8 0.1342 -1.996 -3 3.2 18.09 23.79 6 15 3000 -47.81 48 0.1391 -1.996 -3 3.2 23.79 29.49 6 15 3000 7.162 0.9211 0.1395 -1.209 -3 3.2 29.49 35.19 6 15 3000 8.193 0.1995 2.822e-05 -0.132 -3 3.2 35.19 40.9 6 15 3000 8.133 0.9983 0.1349 -1.181 -3.2 4.7 0 6.69 6 15 3000 2.511 0.3167 0.09085 -0.7407 -3.2 4.7 6.69 12.39 6 15 3000 3.297 0.2091 6.258e-05 -0.2755 -3.2 4.7 12.39 18.09 6 15 3000 1.85 2.281 0.1042 -1.635 -3.2 4.7 18.09 23.79 6 15 3000 3.869 1.001 0.09955 -1.266 -3.2 4.7 23.79 29.49 6 15 3000 -23.98 24.11 0.1057 -1.988 -3.2 4.7 29.49 35.19 6 15 3000 5.403 0.2371 1.5e-05 -0.3177 -3.2 4.7 35.19 40.9 6 15 3000 5.753 0.2337 0.0002982 -0.3108 diff --git a/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt b/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt deleted file mode 100644 index 4c080015..00000000 --- a/analysis/data/jer/Summer16_25nsV1_MC_SF_AK4PFchs.txt +++ /dev/null @@ -1,27 +0,0 @@ -{1 JetEta 0 None ScaleFactor} --5.191 -3.139 3 1.1922 1.0434 1.341 --3.139 -2.964 3 1.1869 1.0626 1.3112 --2.964 -2.853 3 1.7788 1.578 1.9796 --2.853 -2.500 3 1.3418 1.1327 1.5509 --2.500 -2.322 3 1.2963 1.0592 1.5334 --2.322 -2.043 3 1.1512 1.0372 1.2652 --2.043 -1.930 3 1.1426 1.0212 1.264 --1.930 -1.740 3 1.1000 0.9921 1.2079 --1.740 -1.305 3 1.1278 1.0292 1.2264 --1.305 -1.131 3 1.1609 1.0584 1.2634 --1.131 -0.783 3 1.1464 1.0832 1.2096 --0.783 -0.522 3 1.1948 1.1296 1.26 --0.522 0.000 3 1.1595 1.095 1.224 -0.000 0.522 3 1.1595 1.095 1.224 -0.522 0.783 3 1.1948 1.1296 1.26 -0.783 1.131 3 1.1464 1.0832 1.2096 -1.131 1.305 3 1.1609 1.0584 1.2634 -1.305 1.740 3 1.1278 1.0292 1.2264 -1.740 1.930 3 1.1000 0.9921 1.2079 -1.930 2.043 3 1.1426 1.0212 1.264 -2.043 2.322 3 1.1512 1.0372 1.2652 -2.322 2.500 3 1.2963 1.0592 1.5334 -2.500 2.853 3 1.3418 1.1327 1.5509 -2.853 2.964 3 1.7788 1.578 1.9796 -2.964 3.139 3 1.1869 1.0626 1.3112 -3.139 5.191 3 1.1922 1.0434 1.341 diff --git a/nanoAOD/scripts/getFiles b/nanoAOD/scripts/getFiles index e26b2621..53d98a58 100755 --- a/nanoAOD/scripts/getFiles +++ b/nanoAOD/scripts/getFiles @@ -9,3 +9,8 @@ mv nanoData-master/* $CMSSW_BASE/src/nano/analysis/data/ rm -rf nanoData-master rm -rf master.zip +wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt -P $CMSSW_BASE/src/nano/analysis/data/jer +wget https://github.com/cms-jet/JRDatabase/raw/master/textFiles/Summer16_25nsV1_MC/Summer16_25nsV1_MC_SF_AK4PFchs.txt -P $CMSSW_BASE/src/nano/analysis/data/jer +wget https://github.com/cms-jet/JECDatabase/raw/master/textFiles/Summer16_23Sep2016V4_MC/Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt -P $CMSSW_BASE/src/nano/analysis/data/jer + + From 01b7ae7f5a1be31c35f86bdf77cfb9654501f2cf Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Mon, 4 Feb 2019 08:04:07 +0900 Subject: [PATCH 22/23] Moved the codes into nanoBase class --- analysis/interface/nanoBase.h | 37 +++++++ analysis/interface/topObjectSelection.h | 44 +-------- analysis/src/nanoBase.cc | 116 ++++++++++++++++++++++ analysis/src/topObjectSelection.cc | 122 +----------------------- 4 files changed, 157 insertions(+), 162 deletions(-) diff --git a/analysis/interface/nanoBase.h b/analysis/interface/nanoBase.h index a251c758..3ae3a8b6 100644 --- a/analysis/interface/nanoBase.h +++ b/analysis/interface/nanoBase.h @@ -16,6 +16,11 @@ #include "nano/external/interface/MuonScaleFactorEvaluator.h" #include "nano/external/interface/ElecScaleFactorEvaluator.h" +#include "nano/external/interface/JetCorrectionUncertainty.h" +#include "nano/external/interface/JetCorrectorParameters.h" +#include "nano/external/interface/JetResolution.h" +#include + #include "nano/external/interface/BTagCalibrationStandalone.h" //#include "nano/external/interface/TopTriggerSF.h" //#include "nano/external/interface/TTbarModeDefs.h" @@ -32,6 +37,17 @@ class nanoBase : public Events nanoBase(TTree *tree, Bool_t isMC=false) : nanoBase(tree,0,0,isMC) {} virtual ~nanoBase(); virtual void Loop() = 0; + + // In uncertainty study we need to switch the kinematic variables of jets + // The following variables are for this switch + // In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on. + // In default, these are same as the original ones, but when a user wants to study systematic uncertainty + // so that he/she needs to switch them to the evaluated ones, + // just touching them in anlalyser class will be okay, and this is for it. + virtual void GetJetMassPt(UInt_t nIdx, + Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, UInt_t unFlag); + + Int_t GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution); // Output Variables TFile* m_output; @@ -44,8 +60,29 @@ class nanoBase : public Events MuonScaleFactorEvaluator m_muonSF; ElecScaleFactorEvaluator m_elecSF; BTagCalibrationReader m_btagSF, m_btagSF_up, m_btagSF_dn; + JetCorrectionUncertainty *m_jecUnc; + JMENano::JetResolution m_jetResObj; + JMENano::JetResolutionScaleFactor m_jetResSFObj; + TRandom3 *m_rndEngine; + + Float_t m_fDRcone_JER; + Float_t m_fResFactorMathcer; Bool_t m_isMC; + + enum { + OptBit_JER_Up = 0, + OptBit_JER_Dn, + OptBit_JES_Up, + OptBit_JES_Dn + }; + + enum { + OptFlag_JER_Up = ( 1 << OptBit_JER_Up ), + OptFlag_JER_Dn = ( 1 << OptBit_JER_Dn ), + OptFlag_JES_Up = ( 1 << OptBit_JES_Up ), + OptFlag_JES_Dn = ( 1 << OptBit_JES_Dn ) + }; }; #endif diff --git a/analysis/interface/topObjectSelection.h b/analysis/interface/topObjectSelection.h index 36b94712..d20d7ae7 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -5,11 +5,6 @@ #include "nano/external/interface/TopTriggerSF.h" //#include "nano/external/interface/TTbarModeDefs.h" -#include "nano/external/interface/JetCorrectionUncertainty.h" -#include "nano/external/interface/JetCorrectorParameters.h" -#include "nano/external/interface/JetResolution.h" -#include - class topObjectSelection : public nanoBase { protected: @@ -21,12 +16,6 @@ class topObjectSelection : public nanoBase UInt_t m_unFlag; - JetCorrectionUncertainty *jecUnc; - - JMENano::JetResolution jetResObj; - JMENano::JetResolutionScaleFactor jetResSFObj; - TRandom3 *rndEngine; - public: // YOU MUST SET UP ALL IN THE BELOW!!! // (SetCutValues() will force you to do it) @@ -71,9 +60,6 @@ class topObjectSelection : public nanoBase Float_t cut_BJetConeSizeOverlap; Float_t *cut_BJetTypeBTag; // For example, set it as cut_BJetTypeBTag = Jet_btagCSVV2; Float_t cut_BJetBTagCut; - - Float_t m_fDRcone_JER; - Float_t m_fResFactorMathcer; public: // Tip: If you want to use your own additional cut with the existing cut, @@ -95,17 +81,6 @@ class topObjectSelection : public nanoBase virtual bool additionalConditionForBJet(UInt_t nIdx, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, Float_t &fJetMass) {return true;}; - // In uncertainty study we need to switch the kinematic variables of jets - // The following variables are for this switch - // In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on. - // In default, these are same as the original ones, but when a user wants to study systematic uncertainty - // so that he/she needs to switch them to the evaluated ones, - // just touching them in anlalyser class will be okay, and this is for it. - virtual void GetJetMassPt(UInt_t nIdx, - Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi); - - Int_t GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution); - public: std::vector muonSelection(); std::vector elecSelection(); @@ -120,30 +95,13 @@ class topObjectSelection : public nanoBase topObjectSelection(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, UInt_t unFlag = 0); topObjectSelection(TTree *tree=0, Bool_t isMC=false, UInt_t unFlag = 0) : topObjectSelection(tree, 0, 0, isMC, unFlag) {} - ~topObjectSelection() { - if ( jecUnc != NULL ) delete jecUnc; - if ( rndEngine != NULL ) delete rndEngine; - } + ~topObjectSelection() {} // In this function you need to set all the cut conditions in the above // If you do not set this function up (so that you didn't set the cuts), the compiler will deny your code, // so you can be noticed that you forgot the setting up. // And you don't need to run this function indivisually; it will be run in the creator of this class. virtual int SetCutValues() = 0; - - enum { - OptBit_JER_Up = 0, - OptBit_JER_Dn, - OptBit_JES_Up, - OptBit_JES_Dn - }; - - enum { - OptFlag_JER_Up = ( 1 << OptBit_JER_Up ), - OptFlag_JER_Dn = ( 1 << OptBit_JER_Dn ), - OptFlag_JES_Up = ( 1 << OptBit_JES_Up ), - OptFlag_JES_Dn = ( 1 << OptBit_JES_Dn ) - }; }; #endif diff --git a/analysis/src/nanoBase.cc b/analysis/src/nanoBase.cc index 5e680e2e..73a04e53 100644 --- a/analysis/src/nanoBase.cc +++ b/analysis/src/nanoBase.cc @@ -33,6 +33,38 @@ nanoBase::nanoBase(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) : m_btagSF.load(calib, BTagEntry::FLAV_B, "comb"); m_btagSF_up.load(calib, BTagEntry::FLAV_B, "comb"); m_btagSF_dn.load(calib, BTagEntry::FLAV_B, "comb"); + + m_jecUnc = NULL; + m_rndEngine = NULL; + + if (isMC) { + std::string env = getenv("CMSSW_BASE"); + + std::string strPathJetResSFObj = env+"/src/nano/analysis/data/jer/" + "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; + std::string strPathJetResObj = env+"/src/nano/analysis/data/jer/" + "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; + + if (!exists_test(strPathJetResSFObj) || !exists_test(strPathJetResObj)) { + std::cout << "Missing data file, run getFiles and try again" << std::endl; + exit(50); + } + + m_jetResSFObj = JMENano::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); + m_jetResObj = JMENano::JetResolution(strPathJetResObj.c_str()); + + m_rndEngine = new TRandom3(12345); + + std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" + //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt"; + "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; + + JetCorrectorParameters JetCorPar(strPathJecUnc, "Total"); + m_jecUnc = new JetCorrectionUncertainty(JetCorPar); + } + + m_fDRcone_JER = 0.4; // For AK4 jets + m_fResFactorMathcer = 3; // According to https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetResolution } nanoBase::~nanoBase() @@ -41,6 +73,90 @@ nanoBase::~nanoBase() m_output->Write(); m_output->Close(); } + + if (m_jecUnc != NULL) delete m_jecUnc; + if (m_rndEngine != NULL) delete m_rndEngine; +} + +// In uncertainty study we need to switch the kinematic variables of jets +// In default, these are same as the original ones, but when a user wants to study systematic uncertainty +// so that he/she needs to switch them to the evaluated ones, +// just touching them in anlalyser class will be okay, and this is for it. +void nanoBase::GetJetMassPt(UInt_t nIdx, + Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, UInt_t unFlag) +{ + Float_t fCorrFactor = 1.0; + + fJetMass = Jet_mass[nIdx]; + fJetPt = Jet_pt[nIdx]; + fJetEta = Jet_eta[nIdx]; + fJetPhi = Jet_phi[nIdx]; + + //if ( m_isMC && ( unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) + if (m_isMC) { + // Evaluating the central part cJER of the factor + JMENano::JetParameters jetPars = {{JMENano::Binning::JetPt, fJetPt}, + {JMENano::Binning::JetEta, fJetEta}, + {JMENano::Binning::Rho, fixedGridRhoFastjetAll}}; + + const double jetRes = m_jetResObj.getResolution(jetPars); // Note: this is relative resolution. + const double cJER = m_jetResSFObj.getScaleFactor(jetPars, + ((unFlag & (OptFlag_JER_Up | OptFlag_JER_Dn)) == 0 ? Variation::NOMINAL : + ((unFlag & OptFlag_JER_Up) != 0 ? Variation::UP : Variation::DOWN))); + + // We need corresponding genJet + Int_t nIdxGen = GetMatchGenJet(nIdx, fJetPt * jetRes); + const double genJetPt = GenJet_pt[nIdxGen]; + + if ( nIdxGen >= 0 ) { + fCorrFactor = (genJetPt+(fJetPt-genJetPt)*cJER)/fJetPt; + } else { + const double smear = m_rndEngine->Gaus(0, 1); + fCorrFactor = (cJER <= 1 ? 1 : 1+smear*jetRes*sqrt(cJER*cJER-1)); + } + + if ((unFlag & (OptFlag_JES_Up | OptFlag_JES_Dn)) != 0) { // JES + // The evaluator needs pT and eta of the current jet + m_jecUnc->setJetPt(fCorrFactor*fJetPt); + m_jecUnc->setJetEta(fJetEta); + + if ((unFlag & OptFlag_JES_Up) != 0) { + fCorrFactor *= 1+m_jecUnc->getUncertainty(true); + } else { + fCorrFactor *= 1-m_jecUnc->getUncertainty(true); + } + } + } + + fJetMass *= fCorrFactor; + fJetPt *= fCorrFactor; +} + + +Int_t nanoBase::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { + UInt_t i; + + double dEta, dPhi, dR; + double dRFound = m_fDRcone_JER; + UInt_t nIdxFound = 999; + + for (i = 0; i < nGenJet; i++) { + dEta = Jet_eta[nIdxJet]-GenJet_eta[i]; + dPhi = std::abs(Jet_phi[nIdxJet]-GenJet_phi[i]); + if (dPhi > (double)M_PI) dPhi -= (double)(2*M_PI); + + dR = std::sqrt(dEta*dEta+dPhi*dPhi); + + if (dR >= (double)m_fDRcone_JER*0.5) continue; + if (dRFound > dR) { + if (std::abs(Jet_pt[nIdxJet] - GenJet_pt[i]) >= m_fResFactorMathcer*fResolution) continue; + + dRFound = dR; + nIdxFound = i; + } + } + + return (dRFound < 0.75*(double)m_fDRcone_JER ? (Int_t)nIdxFound : -1); } void nanoBase::Loop() diff --git a/analysis/src/topObjectSelection.cc b/analysis/src/topObjectSelection.cc index 78fc93cd..9ead3072 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -3,34 +3,8 @@ using std::vector; topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, UInt_t unFlag) : - nanoBase(tree, had, hadTruth, isMC), m_unFlag(unFlag), jecUnc(NULL), rndEngine(NULL) + nanoBase(tree, had, hadTruth, isMC), m_unFlag(unFlag) { - //if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn | OptFlag_JER_Up | OptFlag_JER_Dn ) ) != 0 ) - if ( m_isMC ) { - std::string env = getenv("CMSSW_BASE"); - - std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jer/" - "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; - jetResSFObj = JMENano::JetResolutionScaleFactor(strPathJetResSFObj.c_str()); - - std::string strPathJetResObj = env + "/src/nano/analysis/data/jer/" - "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; - jetResObj = JMENano::JetResolution(strPathJetResObj.c_str()); - - rndEngine = new TRandom3(12345); - - if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { - std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" - //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt"; - "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; - - JetCorrectorParameters JetCorPar(strPathJecUnc, "Total"); - jecUnc = new JetCorrectionUncertainty(JetCorPar); - } - } - - m_fDRcone_JER = 0.4; // For AK4 jets - m_fResFactorMathcer = 3; // According to https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetResolution } vector topObjectSelection::elecSelection() { @@ -154,7 +128,7 @@ vector topObjectSelection::jetSelection() { vector jets; for (UInt_t i = 0; i < nJet; ++i){ Float_t fJetMass, fJetPt, fJetEta, fJetPhi; - GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi); + GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi, m_unFlag); if (fJetPt < cut_JetPt) continue; if (std::abs(fJetEta) > cut_JetEta) continue; @@ -218,7 +192,7 @@ vector topObjectSelection::bjetSelection() { vector bjets; for (UInt_t i = 0; i < nJet; ++i ) { Float_t fJetMass, fJetPt, fJetEta, fJetPhi; - GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi); + GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi, m_unFlag); if (fJetPt < cut_BJetPt) continue; if (std::abs(fJetEta) > cut_BJetEta) continue; @@ -242,93 +216,3 @@ vector topObjectSelection::bjetSelection() { } -// In uncertainty study we need to switch the kinematic variables of jets -// The following variables are for this switch -// In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on. -// In default, these are same as the original ones, but when a user wants to study systematic uncertainty -// so that he/she needs to switch them to the evaluated ones, -// just touching them in anlalyser class will be okay, and this is for it. -void topObjectSelection::GetJetMassPt(UInt_t nIdx, - Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi) -{ - Float_t fCorrFactor = 1.0; - - fJetMass = Jet_mass[ nIdx ]; - fJetPt = Jet_pt[ nIdx ]; - fJetEta = Jet_eta[ nIdx ]; - fJetPhi = Jet_phi[ nIdx ]; - - //if ( m_isMC && ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) - if ( m_isMC ) { - // Evaluating the central part cJER of the factor - JMENano::JetParameters jetPars = {{JMENano::Binning::JetPt, fJetPt}, - {JMENano::Binning::JetEta, fJetEta}, - {JMENano::Binning::Rho, fixedGridRhoFastjetAll}}; - - const double jetRes = jetResObj.getResolution(jetPars); // Note: this is relative resolution. - const double cJER = jetResSFObj.getScaleFactor(jetPars, - ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn ) ) == 0 ? Variation::NOMINAL : - ( ( m_unFlag & OptFlag_JER_Up ) != 0 ? Variation::UP : Variation::DOWN ) )); - - // We need corresponding genJet - //Int_t nIdxGen = Jet_genJetIdx[ nIdx ]; - Int_t nIdxGen = GetMatchGenJet(nIdx, fJetPt * jetRes); - const double genJetPt = GenJet_pt[ nIdxGen ]; - - // JER (nominal and up and down) - apply scaling method if matched genJet is found, - // apply gaussian smearing method if unmatched - /*if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool - std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) - fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt);*/ - if ( nIdxGen >= 0 ) { - fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt; - } else { - const double smear = rndEngine->Gaus(0, 1); - fCorrFactor = ( cJER <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJER * cJER - 1) ); - } - - if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { // JES - // The evaluator needs pT and eta of the current jet - jecUnc->setJetPt(fCorrFactor * fJetPt); - jecUnc->setJetEta(fJetEta); - - if ( ( m_unFlag & OptFlag_JES_Up ) != 0 ) { - fCorrFactor *= 1 + jecUnc->getUncertainty(true); - } else { - fCorrFactor *= 1 - jecUnc->getUncertainty(true); - } - } - } - - fJetMass *= fCorrFactor; - fJetPt *= fCorrFactor; -} - - -Int_t topObjectSelection::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) { - UInt_t i; - - double dEta, dPhi, dR; - double dRFound = m_fDRcone_JER; - UInt_t nIdxFound = 999; - - for ( i = 0 ; i < nGenJet ; i++ ) { - dEta = Jet_eta[ nIdxJet ] - GenJet_eta[ i ]; - dPhi = std::abs(Jet_phi[ nIdxJet ] - GenJet_phi[ i ]); - if ( dPhi > (double)M_PI ) dPhi -= (double)( 2 * M_PI ); - - dR = std::sqrt(dEta * dEta + dPhi * dPhi); - - if ( dR >= (double)m_fDRcone_JER * 0.5 ) continue; - if ( dRFound > dR ) { - if ( std::abs(Jet_pt[ nIdxJet ] - GenJet_pt[ i ]) >= m_fResFactorMathcer * fResolution ) continue; - - dRFound = dR; - nIdxFound = i; - } - } - - return ( dRFound < 0.75 * (double)m_fDRcone_JER ? (Int_t)nIdxFound : -1 ); -} - - From f41580468ea2c7b7b5a711f85dc72fb48e824253 Mon Sep 17 00:00:00 2001 From: Byeonghak Ko Date: Mon, 4 Feb 2019 14:19:45 +0900 Subject: [PATCH 23/23] Checking the existence of JEC data files --- analysis/src/nanoBase.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/analysis/src/nanoBase.cc b/analysis/src/nanoBase.cc index 73a04e53..f40e20f6 100644 --- a/analysis/src/nanoBase.cc +++ b/analysis/src/nanoBase.cc @@ -44,8 +44,12 @@ nanoBase::nanoBase(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) : "Summer16_25nsV1_MC_SF_AK4PFchs.txt"; std::string strPathJetResObj = env+"/src/nano/analysis/data/jer/" "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt"; + std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" + "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; - if (!exists_test(strPathJetResSFObj) || !exists_test(strPathJetResObj)) { + if (!exists_test(strPathJetResSFObj) || !exists_test(strPathJetResObj) || + !exists_test(strPathJecUnc)) + { std::cout << "Missing data file, run getFiles and try again" << std::endl; exit(50); } @@ -55,9 +59,6 @@ nanoBase::nanoBase(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) : m_rndEngine = new TRandom3(12345); - std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" - //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt"; - "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; JetCorrectorParameters JetCorPar(strPathJecUnc, "Total"); m_jecUnc = new JetCorrectionUncertainty(JetCorPar);