diff --git a/CommonTools/doc/example.info b/CommonTools/doc/example.info
index 5fbf80d521e5bfdbd9b775d0affc28e39eb5ebef..c10f87d6c82e8c56058c4880de5feac0736d8c31 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/CMakeLists.txt b/JEC/CMakeLists.txt
index 798c500daa899a84e64ac134c0ad8806be7b4742..e20c2a790374849a9edb21d33ca9aaf435884f48 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/bin/applyJERsmearing.cc b/JEC/bin/applyJERsmearing.cc
index 95c93348096bbe0af464bba6b9ae25c615249b5e..75d7f7793463493ae24a20980148589dfc9fb75b 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);
@@ -110,26 +111,43 @@ 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 = [&SF,&vars](const RecJet& recjet) {
+ auto getSFs = [year,&SF,&vars](const RecJet& recJet, ostream& cout) {
vector<float> sfs;
- for (const string var: vars)
- sfs.push_back(SF->evaluate({recjet.p4.Eta(), var}));
+ 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 JER->evaluate({recjet.p4.Eta(), recjet.CorrPt(), pu->rho});
- };
for (DT::Looper looper(tIn); looper(); ++looper) {
[[ maybe_unused ]]
@@ -143,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;
@@ -151,59 +186,65 @@ 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);;
- vector<float> sfs = getSFs(*rec_it);
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, cout);
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;
- applySFs(*rec_it, sfs);
+ applySFs(*rec_it, sfs, cout);
}
+ 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) {
+ cout << "---------------" << endl;
- float resolution = getJER(*fake_it);
- normal_distribution<float> d{0,resolution};
- 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;
+ float resolution = Evaluate(JER, cout, *fake_it, pu->rho);;
+ vector<float> sfs = getSFs(*fake_it, cout);
- vector<float> sfs = getSFs(*fake_it);
for (float& sf: sfs) {
- cout << "|- " << sf;
- sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
- cout << " -> " << sf << endl;
+ 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);
+ 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 << "|- applying SFs" << endl;
-
- applySFs(*fake_it, sfs);
+ applySFs(*fake_it, sfs, cout);
}
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)
@@ -244,10 +285,13 @@ 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`)")
.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/JEC/bin/applyJEScorrections.cc b/JEC/bin/applyJEScorrections.cc
index 6978db869af2f519921eee09c0821c1f3b8fb6dc..30b8c7932a393ae45ae2dd26e7ecf4e4d594cd5d 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");
@@ -125,8 +126,7 @@ void applyJEScorrections
BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected JES corrections", tIn) );
// nominal value
- auto corr = nomSF->evaluate({recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho});
- cout << recJet.p4.Pt() << ' ' << corr << '\n';
+ float corr = Evaluate(nomSF, cout, recJet, pu->rho, {}, *recEvt);
recJet.scales.front() = corr;
if (!applySyst) continue;
@@ -134,12 +134,11 @@ void applyJEScorrections
// load uncertainties in JECs
recJet.scales.reserve(JES_variations.size()*2+1);
for (const correction::Correction::Ref& varSF: varSFs) {
- auto var = varSF->evaluate({recJet.p4.Eta(), recJet.p4.Pt()}); /// \todo CorrPt?
+ 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 478b82db3e1c78e048f486a7be9b465b7be7d061..fbbe77ffeef2cdabc9e1850c2bfdee4d0bfef80c 100644
--- a/JEC/bin/common.h
+++ b/JEC/bin/common.h
@@ -2,14 +2,19 @@
#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>
+#include <correction.h>
#include <darwin.h>
namespace DAS::JetEnergy {
@@ -26,20 +31,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)
@@ -153,12 +155,74 @@ 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 with DAS objects from
+/// correctionlib with exceptions.
+///
+/// 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
+ 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 {
+ 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();
+ /// \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);
+}
+
////////////////////////////////////////////////////////////////////////////////
/// Determine the R values for which JetMET provides corrections.
int GetR (const Darwin::Tools::UserInfo& metainfo)
diff --git a/JEC/bin/getJMEtables.cc b/JEC/bin/getJMEtables.cc
index ee41b3f419c4f3ec8280550d6e31c05f5af39b12..4c8b1f432d9ec1f0015d191af0d1c820bd9aea16 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,12 @@ void getJMEtable
inputs.push_back(named_inputs.at(input.name()));
// get and fill the value
+
double value = sf->evaluate(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 +191,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 5b759e1d55d7f33259b658a18d57bdc646eae30f..175a6c45f11d6a06e04818b6bcd34b9b5a2292b3 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
diff --git a/JEC/interface/JMEmatching.h b/JEC/interface/JMEmatching.h
index ba9789dabb8e573c389f6e91b1a33484432e6e3d..bab02c670696b3167ad333a9cd53c4e0d084d375 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);
};
diff --git a/JEC/test/CMakeLists.txt b/JEC/test/CMakeLists.txt
index 8331da218a6a6d79f96b1b6f91de5645342fc651..1db38186f7a89a8823bfb44b94b5d880d845fc88 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()
diff --git a/Ntupliser/test/runAnalysis.sh b/Ntupliser/test/runAnalysis.sh
index bffcba6923cc04d1d0ef713055bf05ee698c5f86..ed89c9243acb8ecf6db2b25bef8851f5dfbcc518 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