 |
DAS
3.0
Das Analysis System
|
|
void | applyJERsmearing (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0}) |
|
void | applyJEScorrections (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0}) |
|
vector< DAS::Uncertainties::Variation > | GetScaleVariations (const DT::MetaInfo &metainfo, bool applySyst) |
|
int | GetNRhoBins (int year) |
|
std::vector< TString > | rhoBins (int year) |
|
std::vector< double > | getBinning (int nBins, float first, float last) |
|
void | assertValidBinning (const std::vector< double > &v) |
|
std::vector< float > | GetLogBinning (float minpt, float maxpt, int nbins) |
|
std::string | GetAlgo (const Darwin::Tools::UserInfo &metainfo) |
|
std::string | GetShortCampaign (const std::string &campaign) |
|
std::string | ScanCorrections (const auto &corrections, const std::string &key="") |
|
template<typename CorrectionType > |
CorrectionType | GetCorrection (const auto &corrections, const std::string &campaign, const std::string &type, const std::string &level, const std::string &suffix) |
|
float | Evaluate (const auto &correction, std::ostream &cout, const RecJet &recJet, const std::optional< float > &rho={}, const std::optional< GenJet > &genJet={}, const std::optional< RecEvent > &recEvt={}, const std::optional< std::string > &systematic={}, const std::optional< float > &JER={}, const std::optional< float > &JERSF={}) |
|
int | GetR (const Darwin::Tools::UserInfo &metainfo) |
|
void | FitResolution (TDirectory *dir, unique_ptr< TH1 > h, const int steering) |
|
void | loopDirsFromFitResponse (TDirectory *dIn, TDirectory *dOut, const int steering) |
|
void | fitJetResolution (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering) |
|
vector< double > | GetMergedBinning (TH1 *h, const int threshold=30) |
|
pair< float, float > | findDLogExtrema (const unique_ptr< TH1 > &DLog, const double mu, const pair< float, float > Range) |
|
double | round_to (double value, double precision=1.0) |
|
TH1 * | CumulativeResponse (TH1 *res1D) |
|
void | FitResponseByBin (TDirectory *dir, unique_ptr< TH2 > res2D, const int steering) |
|
void | loopDirsFromGetResponse (TDirectory *dIn, TDirectory *dOut, const int steering, const DT::Slice slice) |
|
void | fitJetResponse (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0}) |
|
bool | passDijetSelection (const auto &dijet, const auto &jets, const Variation &v, const double alpha=0.3) |
|
void | getDijetBalance (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0}) |
|
void | getJetResponse (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0}) |
|
void | getJMEtable (const fs::path input, const fs::path output, const pt::ptree &config, const int steering) |
|
vector< double > | GetBinning (TAxis *axis) |
|
|
static const std::vector< double > | pt_JERC_edges = {15, 17, 20, 23, 27, 30, 35, 40, 45, 57, 72, 90, 120, 150, 200, 300, 400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500} |
|
static const int | nPtJERCbins = pt_JERC_edges.size()-1 |
|
static const std::vector< double > | abseta_edges = { 0.000, 0.261, 0.522, 0.783, 1.044, 1.305, 1.566, 1.740, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 2.964, 3.139, 3.489, 3.839, 5.191 } |
|
static const int | nAbsEtaBins = abseta_edges.size()-1 |
|
static const std::vector< TString > | absetaBins = MakeTitle(abseta_edges, "|#eta|", false, true, [](double v) { return Form("%.3f", v);}) |
|
static const std::map< int, std::vector< double > > | rho_edges |
|
static const std::vector< double > | eta_unc_edges = {-5.4, -5.0, -4.4, -4.0, -3.5, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.2 , 2.4 , 2.6 , 2.8 , 3.0, 3.5 , 4.0 , 4.4 , 5.0 , 5.4} |
|
static const int | nEtaUncBins = eta_unc_edges.size()-1 |
|
static const std::vector< std::string > | JES_variations {"AbsoluteStat", "AbsoluteScale", "AbsoluteMPFBias", "Fragmentation", "SinglePionECAL", "SinglePionHCAL", "FlavorQCD", "RelativeJEREC1", "RelativeJEREC2", "RelativeJERHF", "RelativePtBB", "RelativePtEC1", "RelativePtEC2", "RelativePtHF", "RelativeBal", "RelativeSample", "RelativeFSR", "RelativeStatFSR", "RelativeStatEC", "RelativeStatHF", "PileUpDataMC", "PileUpPtRef", "PileUpPtBB", "PileUpPtEC1", "PileUpPtEC2", "PileUpPtHF"} |
|
static const float | w = 0.8 |
|
static const auto | deps = numeric_limits<double>::epsilon() |
|
◆ applyJERsmearing()
void DAS::JetEnergy::applyJERsmearing |
( |
const vector< fs::path > & |
inputs, |
|
|
const fs::path & |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice = {1,0} |
|
) |
| |
Get JER smearing corrections in the form of factors to store in the scales
member of RecJet
.
- Todo:
- investigate JER uncertainty breakdown (see Core#35)
- Todo:
- fix in case of offset (e.g. JES applied on MC)
- Parameters
-
inputs | input ROOT files (n-tuples) |
output | name of output ROOT files (n-tuples) |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
slice | slices for running |
54 cout << __func__ <<
' ' << slice <<
" start" << endl;
57 auto tIn = flow.GetInputTree(slice);
58 auto [fOut, tOut] = flow.GetOutput(output);
62 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
63 auto R =
GetR(metainfo);
64 if (!
isMC) BOOST_THROW_EXCEPTION(
DE::BadInput(
"Only MC may be used as input.",
65 make_unique<TFile>(
inputs.front().c_str() )) );
66 auto year = metainfo.Get<
int>(
"flags",
"year");
67 JMEmatching<vector<RecJet>, vector<GenJet>>::maxDR = (R/10.)/2;
69 auto genEvt = flow.GetBranchReadOnly<GenEvent>(
"genEvent");
70 auto recEvt = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
71 auto pu = flow.GetBranchReadOnly<PileUp >(
"pileup" );
72 auto genJets = flow.GetBranchReadOnly <vector<GenJet>>(
"genJets");
73 auto recJets = flow.GetBranchReadWrite<vector<RecJet>>(
"recJets");
75 const auto& tables =
config.get<fs::path>(
"corrections.JER.tables");
76 metainfo.Set<fs::path>(
"corrections",
"JER",
"tables", tables);
77 if (!fs::exists(tables))
78 BOOST_THROW_EXCEPTION( fs::filesystem_error(
"File does not exist.",
79 tables, make_error_code(errc::no_such_file_or_directory)) );
80 const auto campaign =
config.get<
string>(
"corrections.JER.campaign");
81 metainfo.Set<
string>(
"corrections",
"JER",
"campaign",
campaign);
82 const auto& core_width =
config.get<
float>(
"corrections.JER.core_width");
84 metainfo.Set<
float>(
"corrections",
"JER",
"core_width", core_width);
88 ControlPlots raw(
"raw");
89 vector<ControlPlots> calib { ControlPlots(
"nominal") };
91 const bool applySyst = steering &
DT::syst;
92 vector<string> vars {
"nom"};
95 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
"JER" +
SysDown);
96 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
"JER" +
SysUp);
98 calib.push_back(ControlPlots(
"JER" +
SysDown));
99 calib.push_back(ControlPlots(
"JER" +
SysUp));
101 vars.push_back(
"up");
102 vars.push_back(
"down");
106 auto cset = correction::CorrectionSet::from_file(tables.string());
107 string suffix =
"AK" + to_string(R) +
"PF" +
algo;
109 auto SF = GetCorrection<correction::Correction::Ref>(
111 JER = GetCorrection<correction::Correction::Ref>(
114 const auto& smear_tables =
config.get<fs::path>(
"corrections.JER.smear");
115 metainfo.Set<fs::path>(
"corrections",
"JER",
"smear", smear_tables);
118 if (!fs::exists(smear_tables))
119 BOOST_THROW_EXCEPTION( fs::filesystem_error(
"File does not exist.",
120 smear_tables, make_error_code(errc::no_such_file_or_directory)) );
121 correction::Correction::Ref smear;
122 bool JME_smearing = smear_tables !=
"/dev/null";
124 auto cset_smear = correction::CorrectionSet::from_file(smear_tables.string());
125 smear = cset_smear->at(
"JERSmear");
129 auto getSFs = [
year,&SF,&vars](
const RecJet& recJet, ostream& cout) {
131 sfs.reserve(vars.size());
132 for (
const string& var: vars)
133 sfs.push_back(
Evaluate(SF, cout, recJet, {}, {}, {}, var));
136 auto applySFs = [applySyst](RecJet& recJet,
const vector<float>& sfs, ostream& cout) {
137 cout <<
"applying SFs: " << recJet.CorrPt();
138 float nom_scale = recJet.scales.front();
139 for (
float& scale: recJet.scales)
140 scale *= sfs.front();
141 float recpt = recJet.CorrPt();
142 cout <<
" -> " << recpt;
145 recJet.scales.push_back(nom_scale*sfs[1]);
146 recJet.scales.push_back(nom_scale*sfs[2]);
147 cout <<
' ' << recpt * sfs[1]
148 <<
' ' << recpt * sfs[2] << endl;
152 for (
DT::Looper looper(tIn); looper(); ++looper) {
155 cout <<
"===============" << endl;
157 raw(*genJets, genEvt->weights.front());
158 raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
160 cout <<
"Running matching" << endl;
161 JMEmatching matching(*recJets, *genJets);
162 auto fake_its = matching.fake_its;
165 for (
const auto& genJet: *genJets)
166 cout <<
'\t' << genJet.p4.Pt();
167 cout <<
"\nRecJets:";
168 for (
const auto& recJet: *recJets)
169 cout <<
'\t' << recJet.CorrPt();
170 cout <<
"\nCouples:";
171 for (
const auto& [recJet,genJet]: matching.match_its)
172 cout <<
'\t' << genJet->p4.Pt() <<
',' << recJet->CorrPt();
174 for (
const auto& genJet: matching.miss_its)
175 cout <<
'\t' << genJet->p4.Pt();
177 for (
const auto& recJet: matching.fake_its)
178 cout <<
'\t' << recJet->CorrPt();
181 cout <<
"Applying scaling smearing for jets matched in the core of the response" << endl;
182 for (
const auto& [rec_it, gen_it]: matching.match_its) {
183 cout <<
"---------------" << endl;
185 float genpt = gen_it->p4.Pt(),
186 recpt = rec_it->CorrPt();
189 resolution =
Evaluate(JER, cout, *rec_it, pu->rho);;
191 if (abs(
response - 1) > core_width * resolution) {
192 cout <<
"match in the tail" << endl;
193 fake_its.push_back(rec_it);
196 cout <<
"match in the core" << endl;
198 vector<float> sfs = getSFs(*rec_it, cout);
199 for (
float&
sf: sfs) {
200 cout <<
"res sf = " <<
sf;
201 sf = max(0.
f, 1+(
sf-1)*(recpt-genpt)/recpt);
202 cout <<
" -> pt sf = " <<
sf << endl;
205 applySFs(*rec_it, sfs, cout);
208 cout <<
"Fakes after identification of matches in tails: ";
209 for (
const auto& recJet: fake_its)
210 cout <<
'\t' << recJet->CorrPt();
213 cout <<
"Applying stochastic smearing for all other jets "
214 "(either unmatched, or matched in the tails of the response)" << endl;
215 for (
auto& fake_it: fake_its) {
216 cout <<
"---------------" << endl;
218 float resolution =
Evaluate(JER, cout, *fake_it, pu->rho);;
219 vector<float> sfs = getSFs(*fake_it, cout);
221 for (
float&
sf: sfs) {
223 sf =
Evaluate(smear, cout, *fake_it, pu->rho, {}, *recEvt, {}, resolution,
sf);
225 normal_distribution<float> d{0,resolution};
226 static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
227 float response = d(m_random_generator);
230 <<
"RecPt=" << fake_it->CorrPt() <<
' '
231 <<
"JER=" << resolution <<
' '
233 <<
"-> " <<
sf << endl;
237 applySFs(*fake_it, sfs, cout);
240 cout <<
"Sorting by descending pt" << endl;
241 sort(recJets->begin(), recJets->end(),
pt_sort);
243 cout <<
"RecJets after smearing:";
244 for (
const auto& recJet: *recJets)
245 cout <<
'\t' << recJet.CorrPt();
248 for (
size_t i = 0; i < calib.size(); ++i) {
249 calib.at(i)(*genJets, genEvt->weights.front());
250 calib.at(i)(*recJets, genEvt->weights.front() * recEvt->weights.front(), i);
253 cout <<
"Filling" << endl;
254 if (steering &
DT::fill) tOut->Fill();
257 metainfo.Set<
bool>(
"git",
"complete",
true);
262 cout << __func__ <<
' ' << slice <<
" end" << endl;
◆ applyJEScorrections()
void DAS::JetEnergy::applyJEScorrections |
( |
const vector< fs::path > |
inputs, |
|
|
const fs::path |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice = {1,0} |
|
) |
| |
Apply JES corrections to data and MC *n*tuples, as well as uncertainties in data (all sources are considered separately)
- Parameters
-
inputs | input ROOT files (n-tuples) |
output | name of output root file containing the histograms |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
slice | slices for running |
44 cout << __func__ <<
' ' << slice <<
" start" << endl;
47 auto tIn = flow.GetInputTree(slice);
48 auto [fOut, tOut] = flow.GetOutput(output);
52 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
53 auto R =
GetR(metainfo);
54 auto year = metainfo.Get<
int>(
"flags",
"year");
56 auto genEvt =
isMC ? flow.GetBranchReadOnly<GenEvent>(
"genEvent") :
nullptr;
57 auto recEvt = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
58 auto pu = flow.GetBranchReadOnly<PileUp >(
"pileup" );
59 auto recJets = flow.GetBranchReadWrite<vector<RecJet>>(
"recJets");
61 const auto tables =
config.get<fs::path>(
"corrections.JES.tables");
62 metainfo.Set<fs::path>(
"corrections",
"JES",
"tables", tables);
63 if (!fs::exists(tables))
64 BOOST_THROW_EXCEPTION( fs::filesystem_error(
"File does not exists.",
65 tables, make_error_code(errc::no_such_file_or_directory)) );
67 metainfo.Set<
string>(
"corrections",
"JES",
"campaign",
campaign);
72 ControlPlots raw(
"raw");
73 vector<ControlPlots> calib { ControlPlots(
"nominal") };
75 const bool applySyst = steering &
DT::syst;
79 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
source +
SysDown);
80 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
source +
SysUp);
87 auto cset = correction::CorrectionSet::from_file(tables.string());
89 string suffix =
"AK" + to_string(R) +
"PF" +
algo;
90 auto nomSF = GetCorrection<correction::CompoundCorrection::Ref>(
103 vector<correction::Correction::Ref> varSFs;
109 varSFs.push_back(GetCorrection<correction::Correction::Ref>(
110 *
cset, short_campaign,
"MC",
source, suffix));
112 for (
DT::Looper looper(tIn); looper(); ++looper) {
116 double weight = recEvt->weights.front();
122 for (
auto& recJet: *recJets) {
125 if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
129 float corr =
Evaluate(nomSF, cout, recJet, pu->rho, {}, *recEvt);
130 recJet.scales.front() = corr;
132 if (!applySyst)
continue;
136 for (
const correction::Correction::Ref& varSF: varSFs) {
137 float var =
Evaluate(varSF, cout, recJet);
138 recJet.scales.push_back(corr*(1-var));
139 recJet.scales.push_back(corr*(1+var));
144 for (
size_t i = 0; i < calib.size(); ++i)
145 calib.at(i)(*recJets,
weight, i);
147 sort(recJets->begin(), recJets->end(),
pt_sort);
148 if (steering &
DT::fill) tOut->Fill();
151 metainfo.Set<
bool>(
"git",
"complete",
true);
156 cout << __func__ <<
' ' << slice <<
" end" << endl;
◆ assertValidBinning()
void DAS::JetEnergy::assertValidBinning |
( |
const std::vector< double > & |
v | ) |
|
|
inline |
66 for (
size_t i = 1; i<v.size(); ++i) {
67 if (v[i] > v[i-1])
continue;
68 std::cerr << i <<
' ' << v[i] <<
' ' << v[i-1] <<
'\n';
◆ CumulativeResponse()
TH1* DAS::JetEnergy::CumulativeResponse |
( |
TH1 * |
res1D | ) |
|
Creating cumulative response histogram. Given a response 1D histogram, it adds to a bin the content of all previous bins.
309 const int NptBins = res1D->GetNbinsX();
310 for (
int n = 1; n <= NptBins; ++n) {
311 cont+=res1D->GetBinContent(n);
312 res1D->SetBinContent(n, cont);
◆ Evaluate()
float DAS::JetEnergy::Evaluate |
( |
const auto & |
correction, |
|
|
std::ostream & |
cout, |
|
|
const RecJet & |
recJet, |
|
|
const std::optional< float > & |
rho = {} , |
|
|
const std::optional< GenJet > & |
genJet = {} , |
|
|
const std::optional< RecEvent > & |
recEvt = {} , |
|
|
const std::optional< std::string > & |
systematic = {} , |
|
|
const std::optional< float > & |
JER = {} , |
|
|
const std::optional< float > & |
JERSF = {} |
|
) |
| |
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.
- Exceptions
-
A | std::invalid_argument is thrown when an unrecognised variable is requested by the tables. It should then be sufficient to edit function. |
- Todo:
- Use
boost::core::demangle()
to provide a more readable function name?
- Parameters
-
correction | simple or compound correction |
cout | output stream, possible redirected |
recJet | for all corrections |
rho | e.g. for smearing |
genJet | for stoch smearing |
recEvt | for the run or event number |
systematic | e.g. for JER |
JER | e.g. for stoch smearing |
JERSF | e.g. for stoch smearing |
176 const std::optional<GenJet>& genJet = {},
177 const std::optional<RecEvent>& recEvt = {},
178 const std::optional<std::string>& systematic = {},
179 const std::optional<float>& JER = {},
180 const std::optional<float>& JERSF = {}
185 cout << correction->name();
186 vector<correction::Variable::Type>
inputs;
189 auto push_back = [&cout,&
inputs](
const auto& v) {
194 for (
const correction::Variable&
input: correction->inputs()) {
195 const string& n =
input.name();
197 if (n ==
"JetPt" ) push_back(recJet.CorrPt());
198 else if (n ==
"JetEta" ) push_back(recJet.p4.Eta());
199 else if (n ==
"JetPhi" ) push_back(recJet.p4.Phi());
200 else if (n ==
"JetEta" ) push_back(recJet.p4.Phi());
201 else if (n ==
"JetA" ) push_back(recJet.area);
202 else if (n ==
"Rho" ) push_back(*
rho);
203 else if (n ==
"systematic") push_back(*systematic);
204 else if (n ==
"GenPt" ) push_back(genJet ? genJet->CorrPt() : -1);
205 else if (n ==
"EventID" ) push_back((
int) recEvt->evtNo);
206 else if (n ==
"run" ) push_back((
float) recEvt->runNo);
207 else if (n ==
"JER" ) push_back(*JER);
208 else if (n ==
"JERSF" ) push_back(*JERSF);
209 else BOOST_THROW_EXCEPTION( invalid_argument(
"`"s + n +
"` is needed by "s
210 + correction->name() +
" but not recognized."s) );
215 auto corr = correction->evaluate(
inputs);
216 cout <<
" -> " << corr << endl;
219 catch (std::runtime_error& e) {
220 const char *
name =
typeid(*correction).name();
222 auto location = boost::source_location(__FILE__, __LINE__,
name);
223 boost::throw_exception(e, location);
◆ findDLogExtrema()
pair<float,float> DAS::JetEnergy::findDLogExtrema |
( |
const unique_ptr< TH1 > & |
DLog, |
|
|
const double |
mu, |
|
|
const pair< float, float > |
Range |
|
) |
| |
Find the positions of the extrema, starting from the mean and within a given range (to avoid regions with too large fluctuations).
80 const int m = DLog->FindBin(mu);
84 for (
int i = m; i >= DLog->FindBin(Range.first); --i)
85 if (DLog->GetBinContent(i) > DLog->GetBinContent(imax)) imax = i;
86 float x1 = DLog->GetBinLowEdge(imax+1);
90 for (
int i = m; i <= DLog->FindBin(Range.second); ++i)
91 if (DLog->GetBinContent(i) < DLog->GetBinContent(imin)) imin = i;
92 float x2 = DLog->GetBinLowEdge(imin);
◆ fitJetResolution()
void DAS::JetEnergy::fitJetResolution |
( |
const fs::path & |
input, |
|
|
const fs::path & |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering |
|
) |
| |
Fit resolution curves from fitJetResponse
- Parameters
-
input | input ROOT file (histograms) |
output | name of output ROOT (histograms) |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
235 cout << __func__ <<
" start" << endl;
237 auto fIn = make_unique<TFile>(
input .c_str(),
"READ"),
238 fOut = make_unique<TFile>(output.c_str(),
"RECREATE");
241 cout << __func__ <<
" stop" << endl;
◆ fitJetResponse()
void DAS::JetEnergy::fitJetResponse |
( |
const fs::path & |
input, |
|
|
const fs::path & |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice = {1,0} |
|
) |
| |
Fit response distributions from getJetResponse
- Parameters
-
input | input ROOT file (histograms) |
output | name of output ROOT (histograms) |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
slice | number and index of slice |
520 cout << __func__ <<
' ' << slice <<
" start" << endl;
522 auto fIn = make_unique<TFile>(
input .c_str(),
"READ"),
523 fOut = make_unique<TFile>(output.c_str(),
"RECREATE");
526 cout << __func__ <<
' ' << slice <<
" stop" << endl;
◆ FitResolution()
void DAS::JetEnergy::FitResolution |
( |
TDirectory * |
dir, |
|
|
unique_ptr< TH1 > |
h, |
|
|
const int |
steering |
|
) |
| |
Fit the resolution obtained from fitJetResponse
in a given eta/rho bin.
- Todo:
- parameters in bins of eta (and rho/mu?)
--> one should probably give the directory with all eta bins as input
172 TString stdPrint = h->GetTitle();
173 stdPrint.ReplaceAll(TString( h->GetName() ),
"");
174 cout << h->GetName() <<
'\t' << stdPrint << endl;
177 ResolutionFit res(h, cout);
◆ FitResponseByBin()
void DAS::JetEnergy::FitResponseByBin |
( |
TDirectory * |
dir, |
|
|
unique_ptr< TH2 > |
res2D, |
|
|
const int |
steering |
|
) |
| |
Projection 2D histograms onto 1D histograms and fit response
Fit is currently performed in 3 stages:
- Fit the Gaussian core to get a better estimation of the mean and sigma.
- Extend the fit right range, transposing to a Single-Sided Crystal Ball fit, using the mean and sigma better define the fit parameters.
- Extend the fit left range, transposing to a Double-Sided Crystal Ball fit.
Fit iteratively on each stage by adjusting the fit range so that you minimize the amount of failed fits.
334 auto h_ptBins = unique_ptr<TH1>(
dynamic_cast<TH1*
>(res2D->ProjectionX(
"h_ptBins", 1, -1)) );
335 h_ptBins->GetXaxis()->SetTitle(
"p_{T}^{gen}");
338 map<TString, unique_ptr<TH1>> hists;
339 for (TString
name: {
"meanCore",
"sigmaCore",
343 "chi2ndf",
"Median"}) {
345 hists.insert( {
name, unique_ptr<TH1>(
dynamic_cast<TH1*
>( h_ptBins->Clone(
name) )) } );
346 hists[
name]->Reset();
348 hists[
name]->SetDirectory(
nullptr);
352 hists.insert( {
name, unique_ptr<TH1>(
dynamic_cast<TH1*
>( h_ptBins->Clone(
name) )) } );
353 hists[
name]->Reset();
355 hists[
name]->SetDirectory(
nullptr);
358 int Npt = res2D->GetXaxis()->GetNbins();
359 for (
int ptbin = 1; ptbin <= Npt; ++ptbin) {
360 double lowEdge = res2D->GetXaxis()->GetBinLowEdge(ptbin),
361 upEdge = res2D->GetXaxis()->GetBinUpEdge (ptbin);
362 const char * level = TString( res2D->GetXaxis()->GetTitle() ).Contains(
"gen") ?
"gen" :
"rec";
363 const char * hlt_level = TString( res2D->GetXaxis()->GetTitle() ).Contains(
"HLT") ?
"HLT" : level;
364 TString title = Form(
"%.0f < p_{T}^{%s} < %.0f", lowEdge, hlt_level, upEdge),
365 stdPrint = TString( res2D->GetTitle() );
366 stdPrint.ReplaceAll(
"zx projection",
"");
367 stdPrint.ReplaceAll(
")",
"");
368 stdPrint +=
", " + title +
")";
369 cout << stdPrint <<
def << endl;
372 auto original = unique_ptr<TH1>(
dynamic_cast<TH1*
>(res2D->ProjectionY(
"original", ptbin, ptbin)) );
373 if (original->GetEntries() < 100) {
374 cout <<
orange <<
"Skipped due to low number of events!" <<
def << endl;
379 auto h = unique_ptr<TH1>( original->Rebin(edges.size()-1,
"h", edges.data()) );
380 const char * h_title = TString( res2D->GetXaxis()->GetTitle() ).Contains(
"HLT") ?
"HLT bin low edge" :
"gen";
381 h->GetXaxis()->SetTitle(Form(
"#frac{p_{T}^{rec}}{p_{T}^{%s}}", h_title));
384 TDirectory * d_pt = dir->mkdir(Form(
"ptbin%d", ptbin), title);
387 h->SetDirectory(d_pt);
388 float normalisation = abs( h->Integral() );
389 h->Scale(1.
f/normalisation,
"width");
397 vector<double> probabilities{0.25, 0.5, 0.75};
398 const int nQuantiles = probabilities.size();
399 vector<double> quantiles(nQuantiles);
400 h->GetQuantiles(nQuantiles, quantiles.data(), probabilities.data());
401 double& median = quantiles[1];
402 double IQRError = M_PI / 2.0 * (quantiles[2] - quantiles[0]) / sqrt(h->GetEntries());
404 TString hlt_title = Form(
"%.0f < p_{T}^{%s} < %.0f", lowEdge, hlt_level, upEdge);
407 auto hc_tmp =
dynamic_cast<TH1*
>( h->Clone(
"hc_tmp") );
409 hc->SetNameTitle(
"h_cumulative", hlt_title);
410 hc->SetDirectory(d_pt);
414 ResponseFit resp(h, cout);
416 hists.at(
"Lcore")->SetBinContent(ptbin, resp.p[ResponseFit::KL]);
417 hists.at(
"Rcore")->SetBinContent(ptbin, resp.p[ResponseFit::KR]);
423 TString suffix = resp.good() ?
"" :
"Fail";
426 hists.at(
"meanCore" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::MU]);
427 hists.at(
"meanCore" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::MU]);
428 hists.at(
"sigmaCore" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::SIGMA]);
429 hists.at(
"sigmaCore" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::SIGMA]);
431 hists.at(
"kLtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::KL]);
432 hists.at(
"kLtail" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::KL]);
433 hists.at(
"nLtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::NL]);
434 hists.at(
"nLtail" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::NL]);
435 hists.at(
"kRtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::KR]);
436 hists.at(
"kRtail" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::KR]);
437 hists.at(
"nRtail" + suffix)->SetBinContent(ptbin, resp.p[ResponseFit::NR]);
438 hists.at(
"nRtail" + suffix)->SetBinError (ptbin, resp.e[ResponseFit::NR]);
440 hists.at(
"chi2ndf" + suffix)->SetBinContent(ptbin, resp.chi2ndf .value_or(0));
441 hists.at(
"chi2ndf" + suffix)->SetBinError (ptbin, resp.chi2ndfErr.value_or(0));
443 hists.at(
"Median" + suffix)->SetBinContent(ptbin, median);
444 hists.at(
"Median" + suffix)->SetBinError(ptbin, IQRError);
448 TString title = res2D->GetTitle();
449 title.ReplaceAll(
"Response distribution",
"");
450 title.ReplaceAll(
"zx projection",
"");
451 for (
auto&& h: hists) {
452 h.second->SetDirectory(dir);
453 h.second->SetTitle(h.second->GetTitle() + title);
◆ GetAlgo()
Determine from metainfo if CHS or PUPPI has been used to reconstruct the jets.
94 if (metainfo.
Find(
"flags",
"labels",
"CHS" ))
return "chs";
95 if (metainfo.
Find(
"flags",
"labels",
"PUPPI"))
return "Puppi";
97 std::cerr <<
orange <<
"Couldn't identify CHS or PUPPI. Running default "
98 "CHS for Run 2 (PUPPI foir Run 3)\n" <<
def;
100 auto year = metainfo.
Get<
int>(
"flags",
"year");
101 if (
year > 2019)
return "Puppi";
102 if (
year > 2014)
return "chs";
◆ getBinning()
std::vector<double> DAS::JetEnergy::getBinning |
( |
int |
nBins, |
|
|
float |
first, |
|
|
float |
last |
|
) |
| |
|
inline |
57 std::vector<double>
bins(nBins+1);
58 for(
int i = 0; i <= nBins; ++i)
59 bins[i] = first + ((last-first)/nBins) * i;
◆ GetBinning()
vector<double> DAS::JetEnergy::GetBinning |
( |
TAxis * |
axis | ) |
|
42 int n = axis->GetNbins();
43 vector<double>
bins(n+1);
44 for (
int i = 0; i <= n; ++i)
45 bins[i] = axis->GetBinLowEdge(i+1);
46 bins[n] = axis->GetBinUpEdge(n);
◆ GetCorrection()
CorrectionType DAS::JetEnergy::GetCorrection |
( |
const auto & |
corrections, |
|
|
const std::string & |
campaign, |
|
|
const std::string & |
type, |
|
|
const std::string & |
level, |
|
|
const std::string & |
suffix |
|
) |
| |
Returns the correction corresponding to a key.
< expecting correction::[Compound]Correction::Ref
- Parameters
-
corrections | from `correctionlib` |
campaign | e.g. Summer19UL18_RunA |
type | MC or DATA |
level | e.g. L1L2L3Res, TimePtEta |
suffix | e.g. AK4chs |
144 namespace al = boost::algorithm;
145 for (
const auto& correction: corrections) {
146 if (!al::starts_with(correction.first,
campaign) ||
147 !al::contains (correction.first, type) ||
148 !al::contains (correction.first, level) ||
149 !al::ends_with (correction.first, suffix))
152 return correction.second;
155 BOOST_THROW_EXCEPTION(
157 "No `"s +
campaign +
"*"s + type +
"*" + level +
"*"s + suffix
158 +
"` correction can be found in the given tables. "s
◆ getDijetBalance()
void DAS::JetEnergy::getDijetBalance |
( |
const vector< fs::path > |
inputs, |
|
|
const fs::path |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice = {1,0} |
|
) |
| |
Get Pt balance from dijet events.
- Parameters
-
inputs | input ROOT files (n-tuples) |
output | name of output root file containing the histograms |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
slice | slices for running |
82 cout << __func__ <<
' ' << slice <<
" start" << endl;
85 auto tIn = flow.GetInputTree(slice);
86 auto [fOut, tOut] = flow.GetOutput(output);
90 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
91 const auto alpha =
config.get<
float>(
"skim.dijet.alpha");
92 metainfo.Set<
float>(
"skim",
"dijet",
"alpha", alpha);
94 auto gEv =
isMC ? flow.GetBranchReadOnly<GenEvent>(
"genEvent"):
nullptr;
95 auto rEv = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
97 auto genJets =
isMC ? flow.GetBranchReadOnly<vector<GenJet>>(
"genJets") :
nullptr;
98 auto recJets = flow.GetBranchReadOnly<vector<RecJet>>(
"recJets");
100 auto genDijet =
isMC ? flow.GetBranchReadOnly<
GenDijet>(
"genDijet") :
nullptr;
101 auto recDijet = flow.GetBranchReadOnly<
RecDijet>(
"recDijet");
107 vector<Balance> recBalances;
108 TDirectory * rec = fOut->mkdir(
"rec");
113 vector<Balance> genBalances;
114 TDirectory * gen =
nullptr;
116 gen = fOut->mkdir(
"gen");
120 TRandom3 r3(metainfo.Seed<39856>(slice));
121 for (
DT::Looper looper(tIn); looper(); ++looper) {
126 auto r = r3.Uniform();
127 int itag = (r < 0.5),
130 for (Balance& bal: recBalances) {
131 auto gw =
isMC ? gEv->Weight(bal.v) : Weight{1.,0},
132 rw = rEv->Weight(bal.v);
135 (abs((itag ? recDijet->first : recDijet->second)->p4.Eta()) < 1.3))
136 bal((itag ? *recDijet->first : *recDijet->second),
137 (iprobe ? *recDijet->first : *recDijet->second), gw*rw);
142 for (Balance& bal: genBalances) {
143 auto gw = gEv->Weight(bal.v);
145 (abs((itag ? genDijet->first : genDijet->second)->p4.Eta()) < 1.3))
146 bal((itag ? *genDijet->first : *genDijet->second),
147 (iprobe ? *genDijet->first : *genDijet->second), gw);
152 for (
auto&
p: recBalances)
p.Write();
153 for (
auto&
p: genBalances)
p.Write();
155 metainfo.Set<
bool>(
"git",
"complete",
true);
157 cout << __func__ <<
' ' << slice <<
" stop" << endl;
◆ getJetResponse()
void DAS::JetEnergy::getJetResponse |
( |
const vector< fs::path > & |
inputs, |
|
|
const fs::path & |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice = {1,0} |
|
) |
| |
Get JER differential resolution and balance.
- Todo:
- with gen level weight too?
- Todo:
- remove upper boundaries from titles of last bins
- Parameters
-
inputs | input ROOT files (n-tuples) |
output | name of output ROOT file (histograms) |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
slice | slices for running |
50 cout << __func__ <<
' ' << slice <<
" start" << endl;
53 auto tIn = flow.GetInputTree(slice);
54 auto [fOut, tOut] = flow.GetOutput(output);
58 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
60 BOOST_THROW_EXCEPTION(
DE::BadInput(
"Only MC may be used as input.", metainfo) );
62 const auto R = metainfo.Get<
int>(
"flags",
"R");
64 JMEmatching<>::maxDR = R/10./2.;
65 cout <<
"Radius for matching: " << JMEmatching<>::maxDR << endl;
67 auto rEv = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
68 const auto hiPU =
config.get<
bool>(
"skims.pileup");
69 auto pileUp = hiPU ? flow.GetBranchReadOnly<PileUp>(
"pileup") :
nullptr;
71 auto recjets = flow.GetBranchReadOnly<vector<RecJet>>(
"recJets");
72 auto genjets = flow.GetBranchReadOnly<vector<GenJet>>(
"genJets");
74 unique_ptr<TH3> inclResp =
makeRespHist(
"inclusive",
"JER");
75 vector<unique_ptr<TH3>> rhobinsResp,
78 const auto year = metainfo.Get<
int>(
"flags",
"year");
81 rhobinsResp.reserve(nRhoBins);
82 mubinsResp.reserve(nRhoBins);
83 for (
int rhobin = 1; rhobin <= nRhoBins; ++rhobin) {
84 rhobinsResp.push_back(
makeRespHist( Form(
"rhobin%d", rhobin) ,
"JER") );
85 mubinsResp.push_back(
makeRespHist( Form(
"mubin%d", rhobin) ,
"JER") );
89 for (
DT::Looper looper(tIn); looper(); ++looper) {
93 auto recWgt = rEv->weights.front();
94 optional<size_t> irho, imu;
97 && *irho < nRhoBins; ++*irho);
98 if (*irho >= nRhoBins)
99 BOOST_THROW_EXCEPTION(
101 *irho, nRhoBins), tIn) );
102 static const size_t imuMax = mubinsResp.size()-1;
103 imu =
static_cast<size_t>(pileUp->GetTrPU())/10;
104 *imu = min(*imu,imuMax);
110 matching.match_its.resize(min(matching.match_its.size(),3lu));
112 for (
auto& [rec_it,gen_it]: matching.match_its) {
113 auto ptrec = rec_it->CorrPt();
114 auto ptgen = gen_it->p4.Pt();
118 auto etarec = abs( rec_it->p4.Eta() );
120 inclResp->Fill(ptgen, etarec,
response, recWgt);
121 if (irho) rhobinsResp.at(*irho)->Fill(ptgen, etarec,
response, recWgt);
122 if (imu ) mubinsResp.at(*imu )->Fill(ptgen, etarec,
response, recWgt);
130 inclResp->SetDirectory(fOut);
131 inclResp->SetTitle(
"Response (inclusive)");
132 inclResp->Write(
"Response");
134 for (
int irho = 0; irho < nRhoBins; ++irho) {
140 TDirectory * d_rho = fOut->mkdir(Form(
"rhobin%d", irho+1), title);
142 rhobinsResp.at(irho)->SetDirectory(fOut);
143 rhobinsResp.at(irho)->SetTitle(
"Response (" + title +
")");
144 rhobinsResp.at(irho)->Write(
"Response");
146 title = Form(
"%d < #mu < %d", irho*10, (irho+1)*10);
148 TDirectory * d_mu = fOut->mkdir(Form(
"mubin%d", irho+1), title);
150 mubinsResp.at(irho)->SetDirectory(fOut);
151 mubinsResp.at(irho)->SetTitle(
"Response (" + title +
")");
152 mubinsResp.at(irho)->Write(
"Response");
155 metainfo.Set<
bool>(
"git",
"complete",
true);
157 cout << __func__ <<
' ' << slice <<
" end" << endl;
◆ getJMEtable()
void DAS::JetEnergy::getJMEtable |
( |
const fs::path |
input, |
|
|
const fs::path |
output, |
|
|
const pt::ptree & |
config, |
|
|
const int |
steering |
|
) |
| |
Obtain corrections and uncertainties in ROOT format for plotting and comparison.
- Todo:
- Compound corrections?
- Todo:
- Use
THn
to get one more dimension
- Todo:
- check binning
- Todo:
- check binning
- Todo:
- check binning
- Parameters
-
input | JetMET JES tables in JSON format |
output | name of output root file |
config | config file from `DTOptions` |
steering | steering parameters from `DTOptions` |
39 cout << __func__ <<
" start" << endl;
44 auto R =
config.get<
int>(
"flags.R");
48 cout << R <<
' ' <<
year <<
' ' <<
algo << endl;
51 auto level =
config.get<
string>(
"corrections.JME.level");
52 string suffix =
"AK" + to_string(R) +
"PF" +
algo;
54 auto cset = correction::CorrectionSet::from_file(
input.string());
55 auto sf = GetCorrection<correction::Correction::Ref>(
59 cout <<
"Declaring histogram" << endl;
61 map<string, correction::Variable::Type> named_inputs;
62 if (level ==
"L1FastJet" || level ==
"PtResolution") {
63 h = make_unique<TH3D>(level.c_str(), level.c_str(),
68 named_inputs[
"JetA"] = M_PI * pow(R*0.1, 2);
69 named_inputs[
"JetPhi"] = 0;
71 else if (level ==
"L2Relative" || level ==
"L2L3Residual")
72 h = make_unique<TH2D>(level.c_str(), level.c_str(),
75 else if (level ==
"PtResolution") {
100 h = make_unique<TH3D>(level.c_str(), level.c_str(),
106 else if (level ==
"ScaleFactor") {
133 h = make_unique<TH1D>(level.c_str(), level.c_str(),
136 h = make_unique<TH2D>(level.c_str(), level.c_str(),
139 named_inputs[
"systematic"] =
"nom";
142 h = make_unique<TH2D>(level.c_str(), level.c_str(),
146 cout <<
"dim = " << h->GetDimension() << endl;
147 for (
int etabin = 1; etabin <= h->GetNbinsX(); ++etabin)
148 for (
int ptbin = 1; ptbin <= h->GetNbinsY(); ++ ptbin)
149 for (
int rhobin = 1; rhobin <= h->GetNbinsZ(); ++rhobin) {
154 named_inputs[
"JetEta"] = h->GetXaxis()->GetBinCenter(etabin);
155 named_inputs[
"JetPt" ] = h->GetYaxis()->GetBinCenter( ptbin);
156 named_inputs[
"Rho" ] = h->GetZaxis()->GetBinCenter(rhobin);
159 vector<correction::Variable::Type>
inputs;
160 for (
const correction::Variable&
input:
sf->inputs())
165 double value =
sf->evaluate(
inputs);
167 cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
169 if (level ==
"L1FastJet" || level ==
"PtResolution")
170 h->SetBinContent(etabin, ptbin, rhobin, value);
171 else if (level ==
"ScaleFactor")
172 h->SetBinContent(etabin, value);
174 h->SetBinContent(etabin, ptbin, value);
178 cout << __func__ <<
" end" << endl;
◆ GetLogBinning()
std::vector<float> DAS::JetEnergy::GetLogBinning |
( |
float |
minpt, |
|
|
float |
maxpt, |
|
|
int |
nbins |
|
) |
| |
|
inline |
Make equidistant binning in log pt.
- Todo:
- this accumulates imprecision --> improve
80 std::vector<float> edges;
85 std::cout << edges.size() << std::endl;
◆ GetMergedBinning()
vector<double> DAS::JetEnergy::GetMergedBinning |
( |
TH1 * |
h, |
|
|
const int |
threshold = 30 |
|
) |
| |
Provide a new binning ensuring a sufficient number of entries in each bin to assume Gaussian uncertainties.
For a gaussian distribution, as it is assumed for the number of events in each bin, the error is given as sqrt(n) where n the number of events in the bin. Assuming both the bin content and error have been normalised, by taking the division of the (BinContent/BinError)^2 ~ n, one can estimate the number of events in that bin.
There is no special treatment for negative weights, as the samples used to develop the code have no such weight. A warning is however displayed in case negative bins are found.
51 const int N = h->GetNbinsX();
52 vector<double> edges(1,0);
55 for (
int i = 1; i <= N; ++i) {
56 auto content = h->GetBinContent(i);
57 if (content == 0)
continue;
60 << h->GetName() <<
def <<
'\n';
61 auto error = h->GetBinError(i);
62 sum += pow(content/error,2);
65 edges.push_back(h->GetBinLowEdge(i+1));
67 if (edges.back() != h->GetBinLowEdge(N+1))
68 edges.push_back(h->GetBinLowEdge(N+1));
◆ GetNRhoBins()
int DAS::JetEnergy::GetNRhoBins |
( |
int |
year | ) |
|
|
inline |
◆ GetR()
Determine the R values for which JetMET provides corrections.
230 auto year = metainfo.
Get<
int>(
"flags",
"year"),
231 R = metainfo.
Get<
int>(
"flags",
"R");
233 std::string warning =
"Not a standard jet size.\n";
235 if (R != 4 && R != 8)
237 return R < 6 ? 4 : 8;
240 if (R != 5 && R != 7)
242 return R < 6 ? 5 : 7;
◆ GetScaleVariations()
Get scale variations
- Note
- one may want to get the Rochester variations as well for Z+jet balance
-
this function may be redundant with
DAS::Unfolding::GetVariations()
33 vector<DAS::Uncertainties::Variation> variations;
34 variations.emplace_back(
"nominal",
"nominal");
35 if (!applySyst)
return variations;
39 BOOST_THROW_EXCEPTION( invalid_argument(groupContents->GetName()) );
42 for (
const TObject * obj: *groupContents) {
43 const auto os =
dynamic_cast<const TObjString*
>(obj);
44 if (!os) BOOST_THROW_EXCEPTION( invalid_argument(obj->GetName()) );
45 TString s = os->GetString();
46 s.ReplaceAll(
SysUp,
"");
◆ GetShortCampaign()
std::string DAS::JetEnergy::GetShortCampaign |
( |
const std::string & |
campaign | ) |
|
Extracts for isntance Summer19UL18
from Summer19UL18_RunA
111 regex r(
"^(Summer|Fall|Autumn|Winter|Spring)[0-9]{2}[A-Za-z0-9]*");
112 smatch campaign_short;
113 if (!regex_search(
campaign, campaign_short, r))
114 BOOST_THROW_EXCEPTION( invalid_argument(
"The campaign could not be identified"));
115 return campaign_short.str();
◆ loopDirsFromFitResponse()
void DAS::JetEnergy::loopDirsFromFitResponse |
( |
TDirectory * |
dIn, |
|
|
TDirectory * |
dOut, |
|
|
const int |
steering |
|
) |
| |
Copy directory structure. Given an input directory dIn
recursively loop over all directory hierarchy/contents and reproduce it in the directory dOut
of the output file. For each TH1
instance found containing the resolution, perform a NSC fit by calling FitResolution()
function.
198 for (
const auto&& obj: *(dIn->GetListOfKeys())) {
199 auto const key =
dynamic_cast<TKey*
>(obj);
200 if (
key->IsFolder() ) {
201 auto ddIn =
dynamic_cast<TDirectory*
>(
key->ReadObj() );
202 if ( TString( ddIn->GetName() ).Contains(
"ptbin") )
continue;
203 TDirectory * ddOut = dOut->mkdir(ddIn->GetName(), ddIn->GetTitle());
205 cout <<
bold << ddIn->GetPath() <<
def << endl;
209 else if ( TString(
key->ReadObj()->GetName() ).Contains(
"sigmaCore") ) {
210 auto res = unique_ptr<TH1>(
dynamic_cast<TH1*
>(
key->ReadObj()) );
211 res->SetDirectory(dOut);
215 if ( TString(res->GetName()).Contains(
"Fail") )
continue;
220 else cout <<
orange <<
"Ignoring " <<
key->ReadObj()->GetName() <<
" of `"
221 <<
key->ReadObj()->ClassName() <<
"` type" <<
def << endl;
223 cout <<
def << flush;
◆ loopDirsFromGetResponse()
void DAS::JetEnergy::loopDirsFromGetResponse |
( |
TDirectory * |
dIn, |
|
|
TDirectory * |
dOut, |
|
|
const int |
steering, |
|
|
const DT::Slice |
slice |
|
) |
| |
Copy directory structure. Given an input file dIn recursively loop over all directory hierarchy/contents and reproduce it in output file dOut. For each TH3 instance found containing the response perform a DCB fit by calling FitResponseByBin() function.
469 for (
const auto&& obj: *(dIn->GetListOfKeys())) {
470 auto const key =
dynamic_cast<TKey*
>(obj);
471 if (
key->IsFolder() ) {
472 auto ddIn =
dynamic_cast<TDirectory*
>(
key->ReadObj() );
473 if (ddIn ==
nullptr)
continue;
474 TDirectory * ddOut = dOut->mkdir(ddIn->GetName(), ddIn->GetTitle());
476 cout <<
bold << ddIn->GetPath() <<
def << endl;
480 else if ( TString(
key->ReadObj()->ClassName() ).Contains(
"TH3") ) {
482 auto hist3D = unique_ptr<TH3>(
dynamic_cast<TH3*
>(
key->ReadObj()) );
483 hist3D->SetDirectory(0);
484 bool isAbsEta = hist3D->GetYaxis()->GetBinLowEdge(1) >= 0;
485 TDirectory * ddOut = dOut->mkdir(
key->GetName(),
key->ReadObj()->GetTitle());
486 for (
int etabin = 1; etabin <= hist3D->GetNbinsY(); ++etabin) {
489 if (slice.second != ieta % slice.first)
continue;
491 float lowedge = hist3D->GetYaxis()->GetBinLowEdge(etabin),
492 upedge = hist3D->GetYaxis()->GetBinUpEdge (etabin);
493 TString title = isAbsEta ? Form(
"%.3f < |#eta| < %.3f", lowedge, upedge)
494 : Form(
"%.3f < #eta < %.3f" , lowedge, upedge);
495 hist3D->GetYaxis()->SetRange(etabin, etabin);
496 auto hist2D = unique_ptr<TH2>(
dynamic_cast<TH2*
>(hist3D->Project3D(
"ZX")) );
497 hist2D->SetTitle( TString( hist2D->GetTitle() ).ReplaceAll(
") ",
", " + title +
")") );
499 TDirectory * dddOut = ddOut->mkdir(Form(
"etabin%d", etabin), title);
500 cout <<
bold << dddOut->GetPath() <<
def << endl;
505 else cout <<
orange <<
"Ignoring " <<
key->ReadObj()->GetName() <<
" of `"
506 <<
key->ReadObj()->ClassName() <<
"` type" <<
def << endl;
◆ passDijetSelection()
bool DAS::JetEnergy::passDijetSelection |
( |
const auto & |
dijet, |
|
|
const auto & |
jets, |
|
|
const Variation & |
v, |
|
|
const double |
alpha = 0.3 |
|
) |
| |
Apply dijet topology selection independently from the gen or rec level. The balance is part of the standard jet calibration procedure at CMS.
- Parameters
-
dijet | dijet system |
jets | full jet vector |
v | systematic variation |
alpha | reject significant extra jets |
53 if (!dijet)
return false;
54 if (
jets.size() < 2)
return false;
57 if (dijet.DeltaPhi() < 2.7 )
return false;
60 if (std::abs(dijet.first ->p4.Eta()) >= 1.3 &&
61 std::abs(dijet.second->p4.Eta()) >= 1.3)
return false;
64 if (
jets.size() > 2) {
65 const auto pt2 =
jets.at(2).CorrPt(v);
66 if (pt2 > alpha * dijet.HT())
return false;
◆ rhoBins()
std::vector<TString> DAS::JetEnergy::rhoBins |
( |
int |
year | ) |
|
|
inline |
43 [](
double v) { return Form(
"%.2f", v); });
◆ round_to()
double DAS::JetEnergy::round_to |
( |
double |
value, |
|
|
double |
precision = 1.0 |
|
) |
| |
|
inline |
Round number value to precision. Recover "lost" precision.
102 return round(value / precision) * precision;
◆ ScanCorrections()
std::string DAS::JetEnergy::ScanCorrections |
( |
const auto & |
corrections, |
|
|
const std::string & |
key = "" |
|
) |
| |
Prints the available corrections to the standard output. If a key is given as 2nd argument, it is highlighted in the list.
- Parameters
-
corrections | corrections from `correctionlib` |
key | key to highlight in the list |
125 s <<
"Available corrections:";
126 for (
const auto& correction: corrections) {
127 bool found = correction.first ==
key;
129 s <<
' ' << correction.first;
◆ abseta_edges
const std::vector<double> abseta_edges = { 0.000, 0.261, 0.522, 0.783, 1.044, 1.305, 1.566, 1.740, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 2.964, 3.139, 3.489, 3.839, 5.191 } |
|
inlinestatic |
JERC binning (taken from JERCProtoLab repo, /macros/common_info/common_binning.hpp)
◆ absetaBins
const std::vector<TString> absetaBins = MakeTitle(abseta_edges, "|#eta|", false, true, [](double v) { return Form("%.3f", v);}) |
|
inlinestatic |
◆ deps
const auto deps = numeric_limits<double>::epsilon() |
|
static |
◆ eta_unc_edges
const std::vector<double> eta_unc_edges = {-5.4, -5.0, -4.4, -4.0, -3.5, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.2 , 2.4 , 2.6 , 2.8 , 3.0, 3.5 , 4.0 , 4.4 , 5.0 , 5.4} |
|
inlinestatic |
binning for JES uncertainties (taken from JES files)
◆ JES_variations
const std::vector<std::string> JES_variations {"AbsoluteStat", "AbsoluteScale", "AbsoluteMPFBias", "Fragmentation", "SinglePionECAL", "SinglePionHCAL", "FlavorQCD", "RelativeJEREC1", "RelativeJEREC2", "RelativeJERHF", "RelativePtBB", "RelativePtEC1", "RelativePtEC2", "RelativePtHF", "RelativeBal", "RelativeSample", "RelativeFSR", "RelativeStatFSR", "RelativeStatEC", "RelativeStatHF", "PileUpDataMC", "PileUpPtRef", "PileUpPtBB", "PileUpPtEC1", "PileUpPtEC2", "PileUpPtHF"} |
|
inlinestatic |
◆ nAbsEtaBins
const int nAbsEtaBins = abseta_edges.size()-1 |
|
inlinestatic |
◆ nEtaUncBins
const int nEtaUncBins = eta_unc_edges.size()-1 |
|
inlinestatic |
◆ nPtJERCbins
const int nPtJERCbins = pt_JERC_edges.size()-1 |
|
inlinestatic |
◆ pt_JERC_edges
const std::vector<double> pt_JERC_edges = {15, 17, 20, 23, 27, 30, 35, 40, 45, 57, 72, 90, 120, 150, 200, 300, 400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500} |
|
inlinestatic |
◆ rho_edges
const std::map<int,std::vector<double> > rho_edges |
|
inlinestatic |
Initial value:= {
{ 2016, {0, 6.69, 12.39, 18.09, 23.79, 29.49, 35.19, 40.9, 999} },
{ 2017, {0, 7.47, 13.49, 19.52, 25.54, 31.57, 37.59, 999} },
{ 2018, {0, 7.35, 13.26, 19.17, 25.08, 30.99, 36.9, 999} },
{ 2024, {0.0, 8.23, 14.49, 20.75, 27.01, 33.27, 39.53, 45.79, 52.05, 999} },
}
- Todo:
- check exact list of variations with JME group
name
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:48
cerr
Definition: Ntupliser_cfg.py:101
bins
Definition: metPhiCorrectionExample.py:89
static const std::map< int, std::vector< double > > rho_edges
Definition: common.h:31
void loopDirsFromGetResponse(TDirectory *dIn, TDirectory *dOut, const int steering, const DT::Slice slice)
Definition: fitJetResponse.cc:464
pt
Definition: jmarExample.py:19
static bool verbose
Definition: Step.h:40
static const char * def
Definition: Step.h:36
float Evaluate(const auto &correction, std::ostream &cout, const RecJet &recJet, const std::optional< float > &rho={}, const std::optional< GenJet > &genJet={}, const std::optional< RecEvent > &recEvt={}, const std::optional< std::string > &systematic={}, const std::optional< float > &JER={}, const std::optional< float > &JERSF={})
Wrapper to evaluate scale factors with DAS objects from correctionlib with exceptions.
Definition: common.h:172
static const char *const ScaleVar
Name of jet rec scales in variations.
Definition: Jet.h:38
static const std::vector< double > eta_unc_edges
binning for JES uncertainties (taken from JES files)
Definition: common.h:48
vector< double > GetMergedBinning(TH1 *h, const int threshold=30)
Definition: fitJetResponse.cc:49
string algo
Definition: jercExample.py:60
Di< GenJet, GenJet > GenDijet
Definition: Di.h:77
const std::string SysUp
Suffix used for "up" uncertainties. Follows the Combine convention.
Definition: Format.h:8
int year
Definition: Ntupliser_cfg.py:63
cset
Definition: jercExample.py:90
std::vector< TString > MakeTitle(const std::vector< double > &edges, const char *v, bool lowEdge, bool highEdge, std::function< const char *(double)> format)
Make a vector of properly formated titles from bin edges.
Definition: binnings.h:13
source
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:39
static const std::vector< double > pt_JERC_edges
Definition: common.h:22
std::string ScanCorrections(const auto &corrections, const std::string &key="")
Definition: common.h:121
p
Definition: Ntupliser_cfg.py:174
f
Definition: Ntupliser_cfg.py:309
TH1 * CumulativeResponse(TH1 *res1D)
Definition: fitJetResponse.cc:306
string key
Definition: jercExample.py:109
void FitResponseByBin(TDirectory *dir, unique_ptr< TH2 > res2D, const int steering)
Definition: fitJetResponse.cc:328
std::string GetShortCampaign(const std::string &campaign)
Extracts for isntance Summer19UL18 from Summer19UL18_RunA
Definition: common.h:108
unique_ptr< TH3 > makeRespHist(TString name, TString study)
Create TH3 histograms for the Jet response and HLT response.
Definition: MakeResponseHistos.h:20
bool pt_sort(const PhysicsObject &j1, const PhysicsObject &j2)
Definition: toolbox.h:18
static const double minpt
Definition: getMNobservables.cc:39
static const std::vector< std::string > JES_variations
Definition: common.h:51
genjets
Definition: Ntupliser_cfg.py:325
string jets
Definition: Ntupliser_cfg.py:41
DAS::Weight weight
Definition: classes.h:11
config
Definition: Ntupliser_cfg.py:317
int GetNRhoBins(int year)
Definition: common.h:39
static const double maxpt
Definition: getMNobservables.cc:39
vector< DAS::Uncertainties::Variation > GetScaleVariations(const DT::MetaInfo &metainfo, bool applySyst)
Definition: Balance.h:30
static const char * orange
Definition: colours.h:6
input
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:35
recjets
Definition: Ntupliser_cfg.py:326
static const int nEtaUncBins
Definition: common.h:49
static const std::vector< double > abseta_edges
JERC binning (taken from JERCProtoLab repo, /macros/common_info/common_binning.hpp)
Definition: common.h:27
static const char * green
Definition: Step.h:33
int GetR(const Darwin::Tools::UserInfo &metainfo)
Determine the R values for which JetMET provides corrections.
Definition: common.h:228
bool passDijetSelection(const auto &dijet, const auto &jets, const Variation &v, const double alpha=0.3)
Definition: getDijetBalance.cc:48
const std::string SysDown
Suffix used for "down" uncertainties. Follows the Combine convention.
Definition: Format.h:10
dictionary isMC
Definition: Ntupliser_cfg.py:61
Di< RecJet, RecJet > RecDijet
Definition: Di.h:80
const float threshold
Definition: sigmoid.h:14
sf
Definition: jercExample.py:112
def inputs
Definition: jercExample.py:118
Generic exception for problematic event (during event loop).
Definition: exceptions.h:63
void loopDirsFromFitResponse(TDirectory *dIn, TDirectory *dOut, const int steering)
Definition: fitJetResolution.cc:193
static const int nPtJERCbins
Definition: common.h:23
campaign
Definition: generate_html.py:86
void FitResolution(TDirectory *dir, unique_ptr< TH1 > h, const int steering)
Fit the resolution obtained from fitJetResponse in a given eta/rho bin.
Definition: fitJetResolution.cc:167
std::string GetAlgo(const Darwin::Tools::UserInfo &metainfo)
Determine from metainfo if CHS or PUPPI has been used to reconstruct the jets.
Definition: common.h:92
Definition: Variation.h:22
rho
Definition: Ntupliser_cfg.py:320
static const int nAbsEtaBins
Definition: common.h:28
static const char * bold
Definition: Step.h:35