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/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/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 42d49f01..d20d7ae7 100644 --- a/analysis/interface/topObjectSelection.h +++ b/analysis/interface/topObjectSelection.h @@ -14,6 +14,8 @@ class topObjectSelection : public nanoBase Float_t b_maxBDiscr_nonb; + UInt_t m_unFlag; + public: // YOU MUST SET UP ALL IN THE BELOW!!! // (SetCutValues() will force you to do it) @@ -73,8 +75,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;}; public: std::vector muonSelection(); @@ -87,8 +92,9 @@ 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(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() {} // In this function you need to set all the cut conditions in the above diff --git a/analysis/src/nanoBase.cc b/analysis/src/nanoBase.cc index 5e680e2e..f40e20f6 100644 --- a/analysis/src/nanoBase.cc +++ b/analysis/src/nanoBase.cc @@ -33,6 +33,39 @@ 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"; + std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/" + "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt"; + + if (!exists_test(strPathJetResSFObj) || !exists_test(strPathJetResObj) || + !exists_test(strPathJecUnc)) + { + 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); + + + 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 +74,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/topEventSelectionDL.cc b/analysis/src/topEventSelectionDL.cc index 88a6166f..cec06b93 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 8ded0c48..9ead3072 100644 --- a/analysis/src/topObjectSelection.cc +++ b/analysis/src/topObjectSelection.cc @@ -2,9 +2,10 @@ 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) +{ +} vector topObjectSelection::elecSelection() { vector elecs; @@ -126,17 +127,20 @@ vector topObjectSelection::jetSelection() { double csvWgtHF_dn = 1.0, csvWgtLF_dn = 1.0, csvWgtC_dn = 1.0, csvWgtTotal_dn = 1.0; vector jets; 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, m_unFlag); + + 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; } if (hasOverLap) continue; - if ( !additionalConditionForJet(i) ) continue; + if ( !additionalConditionForJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue; auto jet = TParticle(); jet.SetMomentum(mom); @@ -187,18 +191,21 @@ 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, m_unFlag); + + 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; } if (hasOverLap) continue; - if ( !additionalConditionForBJet(i) ) continue; + if ( !additionalConditionForBJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue; auto bjet = TParticle(); bjet.SetMomentum(mom); @@ -207,3 +214,5 @@ vector topObjectSelection::bjetSelection() { } return bjets; } + + 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..8edf1008 --- /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" // MODIFIED +//#include "FWCore/MessageLogger/interface/MessageLogger.h" // MODIFIED + + +//---------- 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..69ee2703 --- /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 JMENano { + 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..29d35275 --- /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 JMENano { + 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 JMENano::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 JMENano { + + 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/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/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 JMENano { + + 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..b65804a4 --- /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/src/JetResolutionObject.cc +// 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 JMENano { + + 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 JMENano::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 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