From b07a6df4551d5012808b1fad5d08285e7832d4a0 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Mon, 10 Feb 2025 22:13:11 +0100
Subject: [PATCH 1/8] extending the application of the JECs to Run 3 (more
input variables are needed)
---
JEC/bin/applyJERsmearing.cc | 22 +++++++++++++---------
JEC/bin/applyJEScorrections.cc | 19 ++++++++++++++++---
JEC/bin/common.h | 20 ++++++++++++++++++++
3 files changed, 49 insertions(+), 12 deletions(-)
diff --git a/JEC/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index 95c93348..b645c5ef 100644
--- a/JEC/bin/applyJERsmearing.cc
+++ b/JEC/bin/applyJERsmearing.cc
@@ -63,6 +63,7 @@ void applyJERsmearing
auto R = GetR(metainfo);
if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
make_unique<TFile>(inputs.front().c_str() )) );
+ auto year = metainfo.Get<int>("flags", "year");
JMEmatching<vector<RecJet>, vector<GenJet>>::maxDR = (R/10.)/2;
auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
@@ -74,7 +75,7 @@ void applyJERsmearing
const auto& tables = config.get<fs::path>("corrections.JER.tables");
metainfo.Set<fs::path>("corrections", "JER", "tables", tables);
if (!fs::exists(tables))
- BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exists.",
+ BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exist.",
tables, make_error_code(errc::no_such_file_or_directory)) );
const auto campaign = config.get<string>("corrections.JER.campaign"); // correction campaign
metainfo.Set<string>("corrections", "JER", "campaign", campaign);
@@ -111,10 +112,13 @@ void applyJERsmearing
*cset, campaign, "MC", "PtResolution", suffix);
// a few lambda functions
- auto getSFs = [&SF,&vars](const RecJet& recjet) {
+ auto getSFs = [year,&SF,&vars](const RecJet& recJet) {
vector<float> sfs;
- for (const string var: vars)
- sfs.push_back(SF->evaluate({recjet.p4.Eta(), var}));
+ for (const string& var: vars) {
+ vector<correction::Variable::Type> inputs {recJet.p4.Eta(), var};
+ if (year > 2020) inputs.insert(inputs.begin()+1, recJet.CorrPt());
+ sfs.push_back(Evaluate(SF, inputs));
+ }
return sfs;
};
auto applySFs = [applySyst](RecJet& recjet, const vector<float>& sfs) {
@@ -158,7 +162,6 @@ void applyJERsmearing
<< "|- response = " << response << '\n'
<< "|- resolution = " << resolution << endl;
- vector<float> sfs = getSFs(*rec_it);
if (abs(response - 1) > core_width * resolution) {
cout << "|- match in the tail" << endl;
fake_its.push_back(rec_it); // left for later, with stochastic smearing
@@ -166,10 +169,11 @@ void applyJERsmearing
}
cout << "|- match in the core" << endl;
+ vector<float> sfs = getSFs(*rec_it);
for (float& sf: sfs) {
- cout << "|- " << sf;
+ cout << "|- res sf = " << sf;
sf = max(0.f, 1+(sf-1)*(recpt-genpt)/recpt);
- cout << " -> " << sf << endl;
+ cout << " -> pt sf = " << sf << endl;
}
cout << "|- applying SFs" << endl;
@@ -191,9 +195,9 @@ void applyJERsmearing
vector<float> sfs = getSFs(*fake_it);
for (float& sf: sfs) {
- cout << "|- " << sf;
+ cout << "|- res sf = " << sf;
sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
- cout << " -> " << sf << endl;
+ cout << " -> pt sf = " << sf << endl;
}
cout << "|- applying SFs" << endl;
diff --git a/JEC/bin/applyJEScorrections.cc b/JEC/bin/applyJEScorrections.cc
index 6978db86..b59111b8 100644
--- a/JEC/bin/applyJEScorrections.cc
+++ b/JEC/bin/applyJEScorrections.cc
@@ -51,6 +51,7 @@ void applyJEScorrections
metainfo.Check(config);
auto isMC = metainfo.Get<bool>("flags", "isMC");
auto R = GetR(metainfo);
+ auto year = metainfo.Get<int>("flags", "year");
auto genEvt = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
@@ -124,17 +125,29 @@ void applyJEScorrections
if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected JES corrections", tIn) );
+ cout << "--------\n"
+ << "|- area = " << recJet.area << '\n'
+ << "|- pt = " << recJet.CorrPt() << '\n'
+ << "|- eta = " << recJet.p4.Eta() << '\n'
+ << "|- phi = " << recJet.p4.Phi() << '\n'
+ << "|- rho = " << pu->rho << endl;
+
// nominal value
- auto corr = nomSF->evaluate({recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho});
- cout << recJet.p4.Pt() << ' ' << corr << '\n';
+ vector<correction::Variable::Type> inputs {recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho};
+ if (year > 2020) inputs.insert(inputs.begin()+2, recJet.p4.Phi()); // \todo check 2022 and 2023
+ float corr = Evaluate(nomSF, inputs);
+ cout << "|- corr = " << corr << endl;
recJet.scales.front() = corr;
if (!applySyst) continue;
// load uncertainties in JECs
recJet.scales.reserve(JES_variations.size()*2+1);
+ inputs = {recJet.p4.Eta(), recJet.CorrPt()};
+ if (year > 2020) inputs.insert(inputs.begin()+1, recJet.p4.Phi()); // \todo check 2022 and 2023
for (const correction::Correction::Ref& varSF: varSFs) {
- auto var = varSF->evaluate({recJet.p4.Eta(), recJet.p4.Pt()}); /// \todo CorrPt?
+ float var = Evaluate(varSF, inputs);
+ cout << "|- var (" << varSF->name() << ") = " << var << endl;
recJet.scales.push_back(corr*(1-var));
recJet.scales.push_back(corr*(1+var));
}
diff --git a/JEC/bin/common.h b/JEC/bin/common.h
index 478b82db..a7d0e728 100644
--- a/JEC/bin/common.h
+++ b/JEC/bin/common.h
@@ -9,7 +9,9 @@
#include "Core/CommonTools/interface/binnings.h"
#include <boost/algorithm/string.hpp>
+#include <boost/assert/source_location.hpp>
+#include <correction.h>
#include <darwin.h>
namespace DAS::JetEnergy {
@@ -159,6 +161,24 @@ CorrectionType GetCorrection (const auto& corrections, //!< from `correctionlib`
);
}
+////////////////////////////////////////////////////////////////////////////////
+/// \brief Wrapper to evaluate scale factors from correctionlib with exceptions.
+float Evaluate (const auto& correction, //!< simple or compound correction
+ const std::vector<correction::Variable::Type>& inputs //!< input variables
+ )
+try {
+ return correction->evaluate(inputs);
+}
+catch (std::runtime_error& e) {
+ const char * name = typeid(*correction).name();
+ auto type = new std::string;
+ *type = boost::core::demangle(name);
+ auto location = boost::source_location(__FILE__, __LINE__-6, type->c_str());
+ /// \note Unlike in most cases where we can use BOOST_THROW_EXCEPTION,
+ /// here, we provide the location by hand.
+ boost::throw_exception(e, location);
+}
+
////////////////////////////////////////////////////////////////////////////////
/// Determine the R values for which JetMET provides corrections.
int GetR (const Darwin::Tools::UserInfo& metainfo)
--
GitLab
From be4c7df081e1e06f7ea1648b766dedb29f100e95 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Thu, 20 Feb 2025 22:16:54 +0100
Subject: [PATCH 2/8] fix memory leak in exception
---
JEC/bin/common.h | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/JEC/bin/common.h b/JEC/bin/common.h
index a7d0e728..c847e463 100644
--- a/JEC/bin/common.h
+++ b/JEC/bin/common.h
@@ -163,6 +163,8 @@ CorrectionType GetCorrection (const auto& corrections, //!< from `correctionlib`
////////////////////////////////////////////////////////////////////////////////
/// \brief Wrapper to evaluate scale factors from correctionlib with exceptions.
+///
+/// \todo Use `boost::core::demangle()` to provide a more readable function name?
float Evaluate (const auto& correction, //!< simple or compound correction
const std::vector<correction::Variable::Type>& inputs //!< input variables
)
@@ -171,9 +173,7 @@ try {
}
catch (std::runtime_error& e) {
const char * name = typeid(*correction).name();
- auto type = new std::string;
- *type = boost::core::demangle(name);
- auto location = boost::source_location(__FILE__, __LINE__-6, type->c_str());
+ auto location = boost::source_location(__FILE__, __LINE__-4, name);
/// \note Unlike in most cases where we can use BOOST_THROW_EXCEPTION,
/// here, we provide the location by hand.
boost::throw_exception(e, location);
--
GitLab
From 122896130fac9afc63ef5a521425ce37d38c0419 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Thu, 20 Feb 2025 22:17:27 +0100
Subject: [PATCH 3/8] minor cpp changes in JME matching
---
JEC/interface/JMEmatching.h | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/JEC/interface/JMEmatching.h b/JEC/interface/JMEmatching.h
index ba9789da..bab02c67 100644
--- a/JEC/interface/JMEmatching.h
+++ b/JEC/interface/JMEmatching.h
@@ -49,11 +49,11 @@ template<typename RecContainer = std::vector<RecJet>,
// 2) declare a few utilities
- typedef pair<size_t,size_t> Indices;
+ using Indices = pair<size_t,size_t>;
auto DeltaR2 = [&recJets,&genJets](const Indices& indices) {
const FourVector& recJet = recJets[indices.first ].p4,
- genJet = genJets[indices.second].p4;
+ & genJet = genJets[indices.second].p4;
return ROOT::Math::VectorUtil::DeltaR2(recJet, genJet);
};
--
GitLab
From f4d2fb9387da140dab4a87d56b854e4a4b9c7590 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Thu, 20 Feb 2025 22:19:25 +0100
Subject: [PATCH 4/8] implementing JME smearing with correctionlib as an
alternative to builtin smearing
---
CommonTools/doc/example.info | 1 +
JEC/bin/applyJERsmearing.cc | 51 ++++++++++++++++++++++++++++-------
Ntupliser/test/runAnalysis.sh | 2 +-
3 files changed, 43 insertions(+), 11 deletions(-)
diff --git a/CommonTools/doc/example.info b/CommonTools/doc/example.info
index 5fbf80d5..c10f87d6 100644
--- a/CommonTools/doc/example.info
+++ b/CommonTools/doc/example.info
@@ -52,6 +52,7 @@ corrections ; for event-by-event corrections
tables ${DARWIN_TABLES}/JER/tables.json
campaign Summer20UL16
core_width 2
+ smear /dev/null
}
normalisation
{
diff --git a/JEC/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index b645c5ef..a1748bcb 100644
--- a/JEC/bin/applyJERsmearing.cc
+++ b/JEC/bin/applyJERsmearing.cc
@@ -111,6 +111,20 @@ void applyJERsmearing
JER = GetCorrection<correction::Correction::Ref>(
*cset, campaign, "MC", "PtResolution", suffix);
+ const auto& smear_tables = config.get<fs::path>("corrections.JER.smear");
+ metainfo.Set<fs::path>("corrections", "JER", "smear", smear_tables);
+
+ // preparing the JER smearing from JME
+ if (!fs::exists(smear_tables))
+ BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exist.",
+ smear_tables, make_error_code(errc::no_such_file_or_directory)) );
+ correction::Correction::Ref smear;
+ bool JME_smearing = smear_tables != "/dev/null";
+ if (JME_smearing) {
+ auto cset_smear = correction::CorrectionSet::from_file(smear_tables.string());
+ smear = cset_smear->at("JERSmear");
+ }
+
// a few lambda functions
auto getSFs = [year,&SF,&vars](const RecJet& recJet) {
vector<float> sfs;
@@ -132,7 +146,7 @@ void applyJERsmearing
}
};
auto getJER = [pu,&JER](const RecJet& recjet) {
- return JER->evaluate({recjet.p4.Eta(), recjet.CorrPt(), pu->rho});
+ return Evaluate(JER, {recjet.p4.Eta(), recjet.CorrPt(), pu->rho});
};
for (DT::Looper looper(tIn); looper(); ++looper) {
@@ -176,32 +190,46 @@ void applyJERsmearing
cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs" << endl;
+ cout << "|- applying SFs";
applySFs(*rec_it, sfs);
}
cout << "Applying stochastic smearing for all other jets "
"(either unmatched, or matched in the tails of the response)" << endl;
for (auto& fake_it: fake_its) {
+ cout << "---------------" << endl;
float resolution = getJER(*fake_it);
- normal_distribution<float> d{0,resolution};
+ vector<float> sfs = getSFs(*fake_it);
+
static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
- float response = d(m_random_generator);
cout << "|- recpt = " << fake_it->CorrPt() << '\n'
- << "|- response = " << response << '\n'
<< "|- resolution = " << resolution << endl;
- vector<float> sfs = getSFs(*fake_it);
for (float& sf: sfs) {
cout << "|- res sf = " << sf;
- sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
+ if (JME_smearing) {
+ int evtID = recEvt->evtNo;
+ float recpt = fake_it->CorrPt(),
+ receta = fake_it->p4.Eta(),
+ genpt = -1, // this means no match for correctionlib
+ rho = pu->rho;
+
+ // JetPt, JetEta, GenPt, Rho, EventID, JER, JERSF
+ vector<correction::Variable::Type> inputs { recpt, receta, genpt, rho, evtID, resolution, sf};
+ sf = Evaluate(smear, inputs);
+ }
+ else {
+ normal_distribution<float> d{0,resolution};
+ float response = d(m_random_generator);
+ cout << "|- response = " << response << '\n';
+ sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
+ }
cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs" << endl;
-
+ cout << "|- applying SFs";
applySFs(*fake_it, sfs);
}
@@ -251,7 +279,10 @@ int main (int argc, char * argv[])
"campaign (e.g. `Summer19UL18_RunA`)")
.arg<float>("core_width", "corrections.JER.core_width",
"width of the Gaussian core in units of sigma "
- "(for the matching algorithm)");
+ "(for the matching algorithm)")
+ .arg<fs::path>("smear", "corrections.JER.smear",
+ "path to JER smearing correctionlib table "
+ "(use /dev/null to rely on a built-in smearing)");
const auto& config = options(argc, argv);
const auto& slice = options.slice();
diff --git a/Ntupliser/test/runAnalysis.sh b/Ntupliser/test/runAnalysis.sh
index bffcba69..ed89c924 100755
--- a/Ntupliser/test/runAnalysis.sh
+++ b/Ntupliser/test/runAnalysis.sh
@@ -43,7 +43,7 @@ if [ "$isMC" = "yes" ]; then
getJetResponse applyJEScorrections.root getJetResponse.root false
fitJetResponse getJetResponse.root fitJetResponse.root
fitJetResolution fitJetResponse.root fitJetResolution.root
- applyJERsmearing applyJEScorrections.root applyJERsmearing.root $JEC_TABLES Summer19UL18 1 -f
+ applyJERsmearing applyJEScorrections.root applyJERsmearing.root $JEC_TABLES Summer19UL18 1 /dev/null -f
getBTagBinnedDiscriminant applyJERsmearing.root getBTagBinnedDiscriminant.root
getBTagPerformance getBTagBinnedDiscriminant.root getBTagPerformance.root
--
GitLab
From 184eb30cfe2a5b8e54083b03f28f94b2d136af34 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Thu, 20 Feb 2025 22:21:01 +0100
Subject: [PATCH 5/8] increasing verbosity in JER smearing
---
JEC/bin/applyJERsmearing.cc | 33 +++++++++++++++++++++++++++++++--
1 file changed, 31 insertions(+), 2 deletions(-)
diff --git a/JEC/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index a1748bcb..54375edd 100644
--- a/JEC/bin/applyJERsmearing.cc
+++ b/JEC/bin/applyJERsmearing.cc
@@ -161,6 +161,23 @@ void applyJERsmearing
JMEmatching matching(*recJets, *genJets);
auto fake_its = matching.fake_its; // we extract it because we will modify it
+ cout << "GenJets:";
+ for (const auto& genJet: *genJets)
+ cout << '\t' << genJet.p4.Pt();
+ cout << "\nRecJets:";
+ for (const auto& recJet: *recJets)
+ cout << '\t' << recJet.CorrPt();
+ cout << "\nCouples:";
+ for (const auto& [recJet,genJet]: matching.match_its)
+ cout << '\t' << genJet->p4.Pt() << ',' << recJet->CorrPt();
+ cout << "\nMisses:";
+ for (const auto& genJet: matching.miss_its)
+ cout << '\t' << genJet->p4.Pt();
+ cout << "\nFakes:";
+ for (const auto& recJet: matching.fake_its)
+ cout << '\t' << recJet->CorrPt();
+ cout << endl;
+
cout << "Applying scaling smearing for jets matched in the core of the response" << endl;
for (const auto& [rec_it, gen_it]: matching.match_its) {
cout << "---------------" << endl;
@@ -190,10 +207,16 @@ void applyJERsmearing
cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs";
+ cout << "|- applying SFs: " << recpt;
applySFs(*rec_it, sfs);
+ cout << " -> " << rec_it->CorrPt() << endl;
}
+ cout << "Fakes after identification of matches in tails: ";
+ for (const auto& recJet: fake_its)
+ cout << '\t' << recJet->CorrPt();
+ cout << endl;
+
cout << "Applying stochastic smearing for all other jets "
"(either unmatched, or matched in the tails of the response)" << endl;
for (auto& fake_it: fake_its) {
@@ -229,13 +252,19 @@ void applyJERsmearing
cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs";
+ cout << "|- applying SFs: " << fake_it->CorrPt();
applySFs(*fake_it, sfs);
+ cout << " -> " << fake_it->CorrPt() << endl;
}
cout << "Sorting by descending pt" << endl;
sort(recJets->begin(), recJets->end(), pt_sort);
+ cout << "RecJets after smearing:";
+ for (const auto& recJet: *recJets)
+ cout << '\t' << recJet.CorrPt();
+ cout << endl;
+
for (size_t i = 0; i < calib.size(); ++i) {
calib.at(i)(*genJets, genEvt->weights.front());
calib.at(i)(*recJets, genEvt->weights.front() * recEvt->weights.front(), i); /// \todo fix in case of offset (e.g. JES applied on MC)
--
GitLab
From cd1502c4d8ac5bd1469324f0a8f1bb4b5490b75e Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Mon, 17 Mar 2025 11:17:20 +0100
Subject: [PATCH 6/8] fix compilation of binaries dependent of correctionlib on
macOS
---
JEC/CMakeLists.txt | 12 ++++++------
JEC/test/CMakeLists.txt | 2 +-
2 files changed, 7 insertions(+), 7 deletions(-)
diff --git a/JEC/CMakeLists.txt b/JEC/CMakeLists.txt
index 798c500d..e20c2a79 100644
--- a/JEC/CMakeLists.txt
+++ b/JEC/CMakeLists.txt
@@ -9,13 +9,13 @@ core_add_library(
JMEmatching)
core_add_executable(applyJERsmearing LIBRARIES correctionlib)
core_add_executable(applyJEScorrections LIBRARIES correctionlib)
-core_add_executable(fitJetResolution NO_EXAMPLE)
-core_add_executable(fitJetResponse NO_EXAMPLE)
+core_add_executable(fitJetResolution NO_EXAMPLE LIBRARIES correctionlib)
+core_add_executable(fitJetResponse NO_EXAMPLE LIBRARIES correctionlib)
core_add_executable(getJMEtables NO_EXAMPLE LIBRARIES correctionlib)
-core_add_executable(getJetResponse NO_EXAMPLE LIBRARIES Objects)
-core_add_executable(getHLTJetResponse NO_EXAMPLE LIBRARIES Objects)
-core_add_executable(getPtBalance NO_EXAMPLE)
-core_add_executable(getDijetBalance NO_EXAMPLE)
+core_add_executable(getJetResponse NO_EXAMPLE LIBRARIES correctionlib Objects)
+core_add_executable(getHLTJetResponse NO_EXAMPLE LIBRARIES correctionlib Objects)
+core_add_executable(getPtBalance NO_EXAMPLE LIBRARIES correctionlib)
+core_add_executable(getDijetBalance NO_EXAMPLE LIBRARIES correctionlib)
#core_add_executable(getZjetBalance NO_EXAMPLE)
add_subdirectory(test)
diff --git a/JEC/test/CMakeLists.txt b/JEC/test/CMakeLists.txt
index 8331da21..1db38186 100644
--- a/JEC/test/CMakeLists.txt
+++ b/JEC/test/CMakeLists.txt
@@ -10,6 +10,6 @@ foreach(X
)
add_executable(${X} ${X}.cc)
target_include_directories(${X} PRIVATE "${CMAKE_SOURCE_DIR}/interface")
- target_link_libraries(${X} ROOT::Hist ROOT::Graf ROOT::HistPainter Darwin::Darwin JEC)
+ target_link_libraries(${X} ROOT::Hist ROOT::Graf ROOT::HistPainter Darwin::Darwin JEC correctionlib)
add_test(NAME ${X} COMMAND ${X})
endforeach()
--
GitLab
From c140ac8ceddaf6d8a6873f90790bcf6f3b6f0ed5 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Wed, 12 Mar 2025 11:42:05 +0100
Subject: [PATCH 7/8] adapting getJMEtables for Run 3 tables (only tested for
pt resolution)
---
JEC/bin/applyJERsmearing.cc | 2 +-
JEC/bin/common.h | 19 +++---
JEC/bin/getJMEtables.cc | 111 +++++++++++++++++++++++++++++-------
JEC/bin/getJetResponse.cc | 20 ++++---
4 files changed, 109 insertions(+), 43 deletions(-)
diff --git a/JEC/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index 54375edd..89709102 100644
--- a/JEC/bin/applyJERsmearing.cc
+++ b/JEC/bin/applyJERsmearing.cc
@@ -305,7 +305,7 @@ int main (int argc, char * argv[])
.arg<fs::path>("tables", "corrections.JER.tables",
"JSON file containing the corrections")
.arg<string>("campaign", "corrections.JER.campaign",
- "campaign (e.g. `Summer19UL18_RunA`)")
+ "campaign (e.g. `Summer19UL18_MC`)")
.arg<float>("core_width", "corrections.JER.core_width",
"width of the Gaussian core in units of sigma "
"(for the matching algorithm)")
diff --git a/JEC/bin/common.h b/JEC/bin/common.h
index c847e463..fce856df 100644
--- a/JEC/bin/common.h
+++ b/JEC/bin/common.h
@@ -28,20 +28,17 @@ inline static const std::vector<TString> absetaBins = MakeTitle(abseta_edges, "|
inline static const std::map<int,std::vector<double>> rho_edges = {
{ 2016, {0, 6.69, 12.39, 18.09, 23.79, 29.49, 35.19, 40.9, 999} }, // EOY16
{ 2017, {0, 7.47, 13.49, 19.52, 25.54, 31.57, 37.59, 999} }, // UL17
- { 2018, {0, 7.35, 13.26, 19.17, 25.08, 30.99, 36.9, 999} } // UL18
-};
+ { 2018, {0, 7.35, 13.26, 19.17, 25.08, 30.99, 36.9, 999} }, // UL18
-inline static const std::map<int,int> nRhoBins = {
- {2016, rho_edges.at(2016).size()-1 },
- {2017, rho_edges.at(2017).size()-1 },
- {2018, rho_edges.at(2018).size()-1 }
+ { 2024, {0.0, 8.23, 14.49, 20.75, 27.01, 33.27, 39.53, 45.79, 52.05, 999} }, // Winter24 (copy of Summer23)
};
-inline static const std::map<int,std::vector<TString>> rhoBins = {
- { 2016, MakeTitle(rho_edges.at(2016), "#rho", false, false, [](double v) { return Form("%.2f", v);}) },
- { 2017, MakeTitle(rho_edges.at(2017), "#rho", false, false, [](double v) { return Form("%.2f", v);}) },
- { 2018, MakeTitle(rho_edges.at(2018), "#rho", false, false, [](double v) { return Form("%.2f", v);}) }
-};
+inline int GetNRhoBins (int year) { return rho_edges.at(year).size()-1; }
+inline std::vector<TString> rhoBins (int year)
+{
+ return MakeTitle(rho_edges.at(year), "#rho", false, false,
+ [](double v) { return Form("%.2f", v); });
+}
////////////////////////////////////////////////////////////////////////////////
/// binning for JES uncertainties (taken from JES files)
diff --git a/JEC/bin/getJMEtables.cc b/JEC/bin/getJMEtables.cc
index ee41b3f4..5789f7b5 100644
--- a/JEC/bin/getJMEtables.cc
+++ b/JEC/bin/getJMEtables.cc
@@ -4,8 +4,7 @@
#include <iostream>
#include <string>
-#include <TH3D.h>
-#include <TF1.h>
+#include <TH3.h>
#include <TFile.h>
#include "Core/CommonTools/interface/binnings.h"
@@ -41,24 +40,21 @@ void getJMEtable
auto fOut = DT::GetOutputFile(output);
+ auto isMC = config.get<bool>("flags.isMC");
auto R = config.get<int>("flags.R");
auto year = config.get<int>("flags.year");
- string algo = "chs"; /// \todo identify algo
+ DT::UserInfo userinfo(config);
+ string algo = GetAlgo(userinfo);
cout << R << ' ' << year << ' ' << algo << endl;
- const auto tag = config.get<string>("corrections.JES.tag"), // correction tag
- level = config.get<string>("corrections.JES.level"); // correction level
- string key = tag + "_" + level + "_AK" + to_string(R) + "PF" + algo;
- cout << key << endl;
+ auto campaign = config.get<string>("corrections.JME.campaign"); // calibration campaign
+ auto level = config.get<string>("corrections.JME.level"); // calibration level
+ string suffix = "AK" + to_string(R) + "PF" + algo;
auto cset = correction::CorrectionSet::from_file(input.string());
- correction::Correction::Ref sf;
- try {
- sf = cset->at(key);
- }
- catch (std::out_of_range& e) {
- BOOST_THROW_EXCEPTION( std::invalid_argument("Nothing corresponding to that key") );
- }
+ auto sf = GetCorrection<correction::Correction::Ref>(
+ *cset, campaign, isMC ? "MC" : "DATA", level, suffix);
+ /// \todo Compound corrections?
cout << "Declaring histogram" << endl;
unique_ptr<TH1> h;
@@ -70,17 +66,79 @@ void getJMEtable
rho_edges.at(year).size()-1,
rho_edges.at(year).data());
named_inputs["JetA"] = M_PI * pow(R*0.1, 2); // Jet area estimation
+ named_inputs["JetPhi"] = 0; /// \todo Use `THn` to get one more dimension
}
else if (level == "L2Relative" || level == "L2L3Residual")
h = make_unique<TH2D>(level.c_str(), level.c_str(),
nAbsEtaBins, abseta_edges .data(),
nPtJERCbins, pt_JERC_edges.data());
+ else if (level == "PtResolution") {
+ // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_PtResolution_AK4PFPuppi (v1)
+ // │ PtResolution for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
+ // │ https://gitlab.cern.ch/cms-jetmet/jerc2json
+ // │ Node counts: Binning: 39, Formula: 305
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ JetEta (real) │
+ // │ │ pseudorapidity of the jet │
+ // │ │ Range: [-5.191, 5.191), overflow ok │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ JetPt (real) │
+ // │ │ pT of the jet before specific correction (for JER and uncertainties: after all corrections │
+ // │ │ applied) │
+ // │ │ Range: [-inf, inf), overflow ok │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ Rho (real) │
+ // │ │ energy density rho (as measure of PU) │
+ // │ │ Range: [0.0, 52.05), overflow ok │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€ â—€ output ─────╮
+ // │ │ correction (real) │
+ // │ │ No description │
+ // │ ╰───────────────────╯
+ h = make_unique<TH3D>(level.c_str(), level.c_str(),
+ nAbsEtaBins, abseta_edges .data(), /// \todo check binning
+ nPtJERCbins, pt_JERC_edges.data(),
+ rho_edges.at(year).size()-1,
+ rho_edges.at(year).data());
+ }
else if (level == "ScaleFactor") {
- h = make_unique<TH1D>(level.c_str(), level.c_str(),
- nAbsEtaBins, abseta_edges .data()); /// \todo check binning
+ // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_ScaleFactor_AK4PFPuppi (v1)
+ // │ ScaleFactor for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
+ // │ https://gitlab.cern.ch/cms-jetmet/jerc2json
+ // │ Node counts: Binning: 43, Category: 1470
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ JetEta (real) │
+ // │ │ pseudorapidity of the jet │
+ // │ │ Range: [-5.191, 5.191), overflow ok │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ JetPt (real) │
+ // │ │ pT of the jet before specific correction (for JER and uncertainties: after all corrections │
+ // │ │ applied) │
+ // │ │ Range: [10.0, 5373.0), overflow ok │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€â”€ â–¶ input ───────────────────────────────────────────╮
+ // │ │ systematic (string) │
+ // │ │ systematics: nom, up, down │
+ // │ │ Values: down, nom, up │
+ // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
+ // │ â•â”€â”€â”€â”€ â—€ output ─────╮
+ // │ │ correction (real) │
+ // │ │ No description │
+ // │ ╰───────────────────╯
+ //
+ if (year < 2020) // run 2
+ h = make_unique<TH1D>(level.c_str(), level.c_str(),
+ nAbsEtaBins, abseta_edges.data()); /// \todo check binning
+ else // run 3
+ h = make_unique<TH2D>(level.c_str(), level.c_str(),
+ nAbsEtaBins, abseta_edges.data(), /// \todo check binning
+ nPtJERCbins, pt_JERC_edges.data());
named_inputs["systematic"] = "nom";
}
- else // assuming uncertainties
+ else // assuming JES uncertainties
h = make_unique<TH2D>(level.c_str(), level.c_str(),
nEtaUncBins, eta_unc_edges.data(),
nPtJERCbins, pt_JERC_edges.data());
@@ -89,6 +147,9 @@ void getJMEtable
for (int etabin = 1; etabin <= h->GetNbinsX(); ++etabin)
for (int ptbin = 1; ptbin <= h->GetNbinsY(); ++ ptbin)
for (int rhobin = 1; rhobin <= h->GetNbinsZ(); ++rhobin) {
+ [[ maybe_unused ]]
+ static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
+
// set values by name
named_inputs["JetEta"] = h->GetXaxis()->GetBinCenter(etabin);
named_inputs["JetPt" ] = h->GetYaxis()->GetBinCenter( ptbin);
@@ -100,11 +161,11 @@ void getJMEtable
inputs.push_back(named_inputs.at(input.name()));
// get and fill the value
- double value = sf->evaluate(inputs);
+ double value = Evaluate(sf, inputs);
cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
- if (level == "L1FastJet")
+ if (level == "L1FastJet" || level == "PtResolution")
h->SetBinContent(etabin, ptbin, rhobin, value);
else if (level == "ScaleFactor")
h->SetBinContent(etabin, value);
@@ -129,10 +190,16 @@ int main (int argc, char * argv[])
auto options = DAS::Options("Obtain JES corrections JetMET tables in ROOT histograms");
options.input ("input" , &input , "JSON file containing the corrections")
.output("output", &output, "output ROOT file")
- .arg<string>("tag", "corrections.JES.tag", "tag (e.g. Summer19UL18_V5_MC)")
- .arg<string>("level", "corrections.JES.level", "level (e.g. L2Relative)")
+ .arg<string>("campaign", "corrections.JME.campaign",
+ "campaign (e.g. Summer19UL18_V5)")
+ .arg<string>("level", "corrections.JME.level", "level (e.g. L2Relative)")
+ /// \todo JES and JER campaigns are often different, but it would be nice to
+ /// 1/ extract the config from the metainfo of a n-tuple
+ /// 2/ use this config as argument to the present executable to get a plot of the tables
+ .arg<bool>("isMC", "flags.isMC", "boolean flag")
.arg<int>("R", "flags.R", "R (x10) parameter in jet clustering algorithm")
- .arg<int>("year", "flags.year", "year (4 digits)");
+ .arg<int>("year", "flags.year", "year (4 digits)")
+ .args ("labels", "flags.labels", "e.g. CHS or PUPPI");
const auto& config = options(argc, argv);
const int steering = options.steering();
diff --git a/JEC/bin/getJetResponse.cc b/JEC/bin/getJetResponse.cc
index 5b759e1d..175a6c45 100644
--- a/JEC/bin/getJetResponse.cc
+++ b/JEC/bin/getJetResponse.cc
@@ -76,10 +76,11 @@ void getJetResponse
mubinsResp; // Response in bins of mu (true pileup)
const auto year = metainfo.Get<int>("flags", "year");
+ const size_t nRhoBins = GetNRhoBins(year);
if (hiPU) {
- rhobinsResp.reserve(nRhoBins.at(year));
- mubinsResp.reserve(nRhoBins.at(year));
- for (int rhobin = 1; rhobin <= nRhoBins.at(year); ++rhobin) {
+ rhobinsResp.reserve(nRhoBins);
+ mubinsResp.reserve(nRhoBins);
+ for (int rhobin = 1; rhobin <= nRhoBins; ++rhobin) {
rhobinsResp.push_back( makeRespHist( Form("rhobin%d", rhobin) , "JER") );
mubinsResp.push_back( makeRespHist( Form( "mubin%d", rhobin) , "JER") );
}
@@ -92,11 +93,12 @@ void getJetResponse
auto recWgt = rEv->weights.front();
optional<size_t> irho, imu;
if (hiPU) {
- static const size_t nrho = static_cast<size_t>(nRhoBins.at(year));
- for (irho = 0; pileUp->rho > rho_edges.at(year).at(*irho+1) && *irho < nrho; ++*irho);
- if (*irho >= nrho)
- BOOST_THROW_EXCEPTION( DE::AnomalousEvent(Form("irho = %lu > nRhoBins.at(%d) = %d",
- *irho, year, nRhoBins.at(year)), tIn) );
+ for (irho = 0; pileUp->rho > rho_edges.at(year).at(*irho+1)
+ && *irho < nRhoBins; ++*irho);
+ if (*irho >= nRhoBins)
+ BOOST_THROW_EXCEPTION(
+ DE::AnomalousEvent(Form("irho = %lu > nRhoBins = %zu",
+ *irho, nRhoBins), tIn) );
static const size_t imuMax = mubinsResp.size()-1;
imu = static_cast<size_t>(pileUp->GetTrPU())/10;
*imu = min(*imu,imuMax);
@@ -129,7 +131,7 @@ void getJetResponse
inclResp->SetTitle("Response (inclusive)");
inclResp->Write("Response");
if (hiPU)
- for (int irho = 0; irho < nRhoBins.at(year); ++irho) {
+ for (int irho = 0; irho < nRhoBins; ++irho) {
/// \todo remove upper boundaries from titles of last bins
--
GitLab
From 73bd609c5d003d6a2e932435f3d5db5bf20ce7e2 Mon Sep 17 00:00:00 2001
From: "Patrick L.S. Connor" <patrick.connor@cern.ch>
Date: Tue, 15 Apr 2025 15:20:27 +0200
Subject: [PATCH 8/8] making DAS::JetEnergy::Evaluate more flexbile (possibly
at the cost of performance) to face the different number of input variables
for the various (versions of the) tables
---
JEC/bin/applyJERsmearing.cc | 84 +++++++++++++---------------------
JEC/bin/applyJEScorrections.cc | 18 +-------
JEC/bin/common.h | 63 +++++++++++++++++++++----
JEC/bin/getJMEtables.cc | 3 +-
4 files changed, 91 insertions(+), 77 deletions(-)
diff --git a/JEC/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index 89709102..75d7f779 100644
--- a/JEC/bin/applyJERsmearing.cc
+++ b/JEC/bin/applyJERsmearing.cc
@@ -126,28 +126,28 @@ void applyJERsmearing
}
// a few lambda functions
- auto getSFs = [year,&SF,&vars](const RecJet& recJet) {
+ auto getSFs = [year,&SF,&vars](const RecJet& recJet, ostream& cout) {
vector<float> sfs;
- for (const string& var: vars) {
- vector<correction::Variable::Type> inputs {recJet.p4.Eta(), var};
- if (year > 2020) inputs.insert(inputs.begin()+1, recJet.CorrPt());
- sfs.push_back(Evaluate(SF, inputs));
- }
+ sfs.reserve(vars.size());
+ for (const string& var: vars)
+ sfs.push_back(Evaluate(SF, cout, recJet, {}, {}, {}, var));
return sfs;
};
- auto applySFs = [applySyst](RecJet& recjet, const vector<float>& sfs) {
- float nom_scale = recjet.scales.front();
- for (float& scale: recjet.scales)
+ auto applySFs = [applySyst](RecJet& recJet, const vector<float>& sfs, ostream& cout) {
+ cout << "applying SFs: " << recJet.CorrPt();
+ float nom_scale = recJet.scales.front();
+ for (float& scale: recJet.scales)
scale *= sfs.front();
+ float recpt = recJet.CorrPt();
+ cout << " -> " << recpt;
if (applySyst) {
- recjet.scales.push_back(nom_scale*sfs[1]);
- recjet.scales.push_back(nom_scale*sfs[2]);
+ recJet.scales.push_back(nom_scale*sfs[1]);
+ recJet.scales.push_back(nom_scale*sfs[2]);
+ cout << ' ' << recpt * sfs[1]
+ << ' ' << recpt * sfs[2] << endl;
}
};
- auto getJER = [pu,&JER](const RecJet& recjet) {
- return Evaluate(JER, {recjet.p4.Eta(), recjet.CorrPt(), pu->rho});
- };
for (DT::Looper looper(tIn); looper(); ++looper) {
[[ maybe_unused ]]
@@ -186,30 +186,23 @@ void applyJERsmearing
recpt = rec_it->CorrPt();
float response = recpt / genpt,
- resolution = getJER(*rec_it);
-
- cout << "|- genpt = " << genpt << '\n'
- << "|- recpt = " << recpt << '\n'
- << "|- response = " << response << '\n'
- << "|- resolution = " << resolution << endl;
+ resolution = Evaluate(JER, cout, *rec_it, pu->rho);;
if (abs(response - 1) > core_width * resolution) {
- cout << "|- match in the tail" << endl;
+ cout << "match in the tail" << endl;
fake_its.push_back(rec_it); // left for later, with stochastic smearing
continue;
}
- cout << "|- match in the core" << endl;
+ cout << "match in the core" << endl;
- vector<float> sfs = getSFs(*rec_it);
+ vector<float> sfs = getSFs(*rec_it, cout);
for (float& sf: sfs) {
- cout << "|- res sf = " << sf;
+ cout << "res sf = " << sf;
sf = max(0.f, 1+(sf-1)*(recpt-genpt)/recpt);
cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs: " << recpt;
- applySFs(*rec_it, sfs);
- cout << " -> " << rec_it->CorrPt() << endl;
+ applySFs(*rec_it, sfs, cout);
}
cout << "Fakes after identification of matches in tails: ";
@@ -222,39 +215,26 @@ void applyJERsmearing
for (auto& fake_it: fake_its) {
cout << "---------------" << endl;
- float resolution = getJER(*fake_it);
- vector<float> sfs = getSFs(*fake_it);
-
- static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
-
- cout << "|- recpt = " << fake_it->CorrPt() << '\n'
- << "|- resolution = " << resolution << endl;
+ float resolution = Evaluate(JER, cout, *fake_it, pu->rho);;
+ vector<float> sfs = getSFs(*fake_it, cout);
for (float& sf: sfs) {
- cout << "|- res sf = " << sf;
- if (JME_smearing) {
- int evtID = recEvt->evtNo;
- float recpt = fake_it->CorrPt(),
- receta = fake_it->p4.Eta(),
- genpt = -1, // this means no match for correctionlib
- rho = pu->rho;
-
- // JetPt, JetEta, GenPt, Rho, EventID, JER, JERSF
- vector<correction::Variable::Type> inputs { recpt, receta, genpt, rho, evtID, resolution, sf};
- sf = Evaluate(smear, inputs);
- }
+ if (JME_smearing) // recJet, rho, genJet, recEvt, syst, JER, JERSF
+ sf = Evaluate(smear, cout, *fake_it, pu->rho, {}, *recEvt, {}, resolution, sf);
else {
normal_distribution<float> d{0,resolution};
+ static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
float response = d(m_random_generator);
- cout << "|- response = " << response << '\n';
sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
+ cout << "smearing "
+ << "RecPt=" << fake_it->CorrPt() << ' '
+ << "JER=" << resolution << ' '
+ << "response = " << response << ' '
+ << "-> " << sf << endl;
}
- cout << " -> pt sf = " << sf << endl;
}
- cout << "|- applying SFs: " << fake_it->CorrPt();
- applySFs(*fake_it, sfs);
- cout << " -> " << fake_it->CorrPt() << endl;
+ applySFs(*fake_it, sfs, cout);
}
cout << "Sorting by descending pt" << endl;
@@ -305,7 +285,7 @@ int main (int argc, char * argv[])
.arg<fs::path>("tables", "corrections.JER.tables",
"JSON file containing the corrections")
.arg<string>("campaign", "corrections.JER.campaign",
- "campaign (e.g. `Summer19UL18_MC`)")
+ "campaign (e.g. `Summer19UL18`)")
.arg<float>("core_width", "corrections.JER.core_width",
"width of the Gaussian core in units of sigma "
"(for the matching algorithm)")
diff --git a/JEC/bin/applyJEScorrections.cc b/JEC/bin/applyJEScorrections.cc
index b59111b8..30b8c793 100644
--- a/JEC/bin/applyJEScorrections.cc
+++ b/JEC/bin/applyJEScorrections.cc
@@ -125,34 +125,20 @@ void applyJEScorrections
if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected JES corrections", tIn) );
- cout << "--------\n"
- << "|- area = " << recJet.area << '\n'
- << "|- pt = " << recJet.CorrPt() << '\n'
- << "|- eta = " << recJet.p4.Eta() << '\n'
- << "|- phi = " << recJet.p4.Phi() << '\n'
- << "|- rho = " << pu->rho << endl;
-
// nominal value
- vector<correction::Variable::Type> inputs {recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho};
- if (year > 2020) inputs.insert(inputs.begin()+2, recJet.p4.Phi()); // \todo check 2022 and 2023
- float corr = Evaluate(nomSF, inputs);
- cout << "|- corr = " << corr << endl;
+ float corr = Evaluate(nomSF, cout, recJet, pu->rho, {}, *recEvt);
recJet.scales.front() = corr;
if (!applySyst) continue;
// load uncertainties in JECs
recJet.scales.reserve(JES_variations.size()*2+1);
- inputs = {recJet.p4.Eta(), recJet.CorrPt()};
- if (year > 2020) inputs.insert(inputs.begin()+1, recJet.p4.Phi()); // \todo check 2022 and 2023
for (const correction::Correction::Ref& varSF: varSFs) {
- float var = Evaluate(varSF, inputs);
- cout << "|- var (" << varSF->name() << ") = " << var << endl;
+ float var = Evaluate(varSF, cout, recJet);
recJet.scales.push_back(corr*(1-var));
recJet.scales.push_back(corr*(1+var));
}
}
- cout << flush;
// CP after corrections
for (size_t i = 0; i < calib.size(); ++i)
diff --git a/JEC/bin/common.h b/JEC/bin/common.h
index fce856df..fbbe77ff 100644
--- a/JEC/bin/common.h
+++ b/JEC/bin/common.h
@@ -2,11 +2,14 @@
#include <cassert>
#include <iostream>
+#include <optional>
#include <regex>
#include <string>
#include <sstream>
#include "Core/CommonTools/interface/binnings.h"
+#include "Core/Objects/interface/Jet.h"
+#include "Core/Objects/interface/Event.h"
#include <boost/algorithm/string.hpp>
#include <boost/assert/source_location.hpp>
@@ -152,27 +155,71 @@ CorrectionType GetCorrection (const auto& corrections, //!< from `correctionlib`
BOOST_THROW_EXCEPTION(
invalid_argument(
"No `"s + campaign + "*"s + type + "*" + level + "*"s + suffix
- + "` correction can be found in the given tables."s
+ + "` correction can be found in the given tables. "s
+ ScanCorrections(corrections)
)
);
}
////////////////////////////////////////////////////////////////////////////////
-/// \brief Wrapper to evaluate scale factors from correctionlib with exceptions.
+/// \brief Wrapper to evaluate scale factors with DAS objects from
+/// correctionlib with exceptions.
///
-/// \todo Use `boost::core::demangle()` to provide a more readable function name?
+/// The implementation is generic, namely for all years, corrections, and
+/// versions of the corrections, which may depend on various variables. The
+/// variables names (e.g. JetPt) are only applicable to JERC. This brute-force
+/// approach is necessary because of the unpredictable changes of conventions.
float Evaluate (const auto& correction, //!< simple or compound correction
- const std::vector<correction::Variable::Type>& inputs //!< input variables
+ std::ostream& cout, //!< output stream, possible redirected
+ const RecJet& recJet, //!< for all corrections
+ const std::optional<float>& rho = {}, //!< e.g. for smearing
+ const std::optional<GenJet>& genJet = {}, //!< for stoch smearing
+ const std::optional<RecEvent>& recEvt = {}, //!< for the run or event number
+ const std::optional<std::string>& systematic = {}, //!< e.g. for JER
+ const std::optional<float>& JER = {}, //!< e.g. for stoch smearing
+ const std::optional<float>& JERSF = {} //!< e.g. for stoch smearing
)
try {
- return correction->evaluate(inputs);
+ using namespace std;
+
+ cout << correction->name();
+ vector<correction::Variable::Type> inputs;
+ inputs.reserve(10);
+
+ auto push_back = [&cout,&inputs](const auto& v) {
+ inputs.push_back(v);
+ cout << '=' << v;
+ };
+
+ for (const correction::Variable& input: correction->inputs()) {
+ const string& n = input.name();
+ cout << ' ' << n;
+ if (n == "JetPt" ) push_back(recJet.CorrPt());
+ else if (n == "JetEta" ) push_back(recJet.p4.Eta());
+ else if (n == "JetPhi" ) push_back(recJet.p4.Phi());
+ else if (n == "JetEta" ) push_back(recJet.p4.Phi());
+ else if (n == "JetA" ) push_back(recJet.area);
+ else if (n == "Rho" ) push_back(*rho);
+ else if (n == "systematic") push_back(*systematic);
+ else if (n == "GenPt" ) push_back(genJet ? genJet->CorrPt() : -1);
+ else if (n == "EventID" ) push_back((int) recEvt->evtNo);
+ else if (n == "run" ) push_back((float) recEvt->runNo);
+ else if (n == "JER" ) push_back(*JER);
+ else if (n == "JERSF" ) push_back(*JERSF);
+ else BOOST_THROW_EXCEPTION( invalid_argument("`"s + n + "` is needed by "s
+ + correction->name() + " but not recognized."s) );
+ /// \exception A `std::invalid_argument` is thrown when an unrecognised variable
+ /// is requested by the tables. It should then be sufficient to edit
+ /// function.
+ }
+ auto corr = correction->evaluate(inputs);
+ cout << " -> " << corr << endl;
+ return corr;
}
catch (std::runtime_error& e) {
const char * name = typeid(*correction).name();
- auto location = boost::source_location(__FILE__, __LINE__-4, name);
- /// \note Unlike in most cases where we can use BOOST_THROW_EXCEPTION,
- /// here, we provide the location by hand.
+ /// \todo Use `boost::core::demangle()` to provide a more readable function name?
+ auto location = boost::source_location(__FILE__, __LINE__, name);
boost::throw_exception(e, location);
}
diff --git a/JEC/bin/getJMEtables.cc b/JEC/bin/getJMEtables.cc
index 5789f7b5..4c8b1f43 100644
--- a/JEC/bin/getJMEtables.cc
+++ b/JEC/bin/getJMEtables.cc
@@ -161,7 +161,8 @@ void getJMEtable
inputs.push_back(named_inputs.at(input.name()));
// get and fill the value
- double value = Evaluate(sf, inputs);
+
+ double value = sf->evaluate(inputs);
cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
--
GitLab