|
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) |
|
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) |
|
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::map< int, int > | nRhoBins |
|
static const std::map< int, std::vector< TString > > | rhoBins |
|
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(
inputs, 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 JMEmatching<vector<RecJet>, vector<GenJet>>::maxDR = (R/10.)/2;
68 auto genEvt = flow.GetBranchReadOnly<GenEvent>(
"genEvent");
69 auto recEvt = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
70 auto pu = flow.GetBranchReadOnly<PileUp >(
"pileup" );
71 auto genJets = flow.GetBranchReadOnly <vector<GenJet>>(
"genJets");
72 auto recJets = flow.GetBranchReadWrite<vector<RecJet>>(
"recJets");
74 const auto& tables =
config.get<fs::path>(
"corrections.JER.tables");
75 metainfo.Set<fs::path>(
"corrections",
"JER",
"tables", tables);
76 if (!fs::exists(tables))
77 BOOST_THROW_EXCEPTION( fs::filesystem_error(
"File does not exists.",
78 tables, make_error_code(errc::no_such_file_or_directory)) );
79 const auto campaign =
config.get<
string>(
"corrections.JER.campaign");
80 metainfo.Set<
string>(
"corrections",
"JER",
"campaign",
campaign);
81 const auto& core_width =
config.get<
float>(
"corrections.JER.core_width");
83 metainfo.Set<
float>(
"corrections",
"JER",
"core_width", core_width);
87 ControlPlots raw(
"raw");
88 vector<ControlPlots> calib { ControlPlots(
"nominal") };
90 const bool applySyst = steering &
DT::syst;
91 vector<string> vars {
"nom"};
94 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
"JER" +
SysDown);
95 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
"JER" +
SysUp);
97 calib.push_back(ControlPlots(
"JER" +
SysDown));
98 calib.push_back(ControlPlots(
"JER" +
SysUp));
100 vars.push_back(
"up");
101 vars.push_back(
"down");
105 auto cset = correction::CorrectionSet::from_file(tables.string());
106 string suffix =
"AK" + to_string(R) +
"PF" +
algo;
108 auto SF = GetCorrection<correction::Correction::Ref>(
110 JER = GetCorrection<correction::Correction::Ref>(
114 auto getSFs = [&SF,&vars](
const RecJet&
recjet) {
116 for (
const string var: vars)
117 sfs.push_back(SF->evaluate({recjet.p4.Eta(), var}));
120 auto applySFs = [applySyst](RecJet&
recjet,
const vector<float>& sfs) {
123 scale *= sfs.front();
130 auto getJER = [pu,&JER](
const RecJet&
recjet) {
134 for (
DT::Looper looper(tIn); looper(); ++looper) {
137 cout <<
"===============" << endl;
139 raw(*genJets, genEvt->weights.front());
140 raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
142 cout <<
"Running matching" << endl;
143 JMEmatching matching(*recJets, *genJets);
144 auto fake_its = matching.fake_its;
146 cout <<
"Applying scaling smearing for jets matched in the core of the response" << endl;
147 for (
const auto& [rec_it, gen_it]: matching.match_its) {
148 cout <<
"---------------" << endl;
150 float genpt = gen_it->p4.Pt(),
151 recpt = rec_it->CorrPt();
154 resolution = getJER(*rec_it);
156 cout <<
"|- genpt = " << genpt <<
'\n'
157 <<
"|- recpt = " << recpt <<
'\n'
158 <<
"|- response = " <<
response <<
'\n'
159 <<
"|- resolution = " << resolution << endl;
161 vector<float> sfs = getSFs(*rec_it);
162 if (abs(
response - 1) > core_width * resolution) {
163 cout <<
"|- match in the tail" << endl;
164 fake_its.push_back(rec_it);
167 cout <<
"|- match in the core" << endl;
169 for (
float&
sf: sfs) {
171 sf = max(0.
f, 1+(
sf-1)*(recpt-genpt)/recpt);
172 cout <<
" -> " <<
sf << endl;
175 cout <<
"|- applying SFs" << endl;
176 applySFs(*rec_it, sfs);
179 cout <<
"Applying stochastic smearing for all other jets "
180 "(either unmatched, or matched in the tails of the response)" << endl;
181 for (
auto& fake_it: fake_its) {
183 float resolution = getJER(*fake_it);
184 normal_distribution<float> d{0,resolution};
185 static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
186 float response = d(m_random_generator);
188 cout <<
"|- recpt = " << fake_it->CorrPt() <<
'\n'
189 <<
"|- response = " <<
response <<
'\n'
190 <<
"|- resolution = " << resolution << endl;
192 vector<float> sfs = getSFs(*fake_it);
193 for (
float&
sf: sfs) {
196 cout <<
" -> " <<
sf << endl;
199 cout <<
"|- applying SFs" << endl;
201 applySFs(*fake_it, sfs);
204 cout <<
"Sorting by descending pt" << endl;
205 sort(recJets->begin(), recJets->end(),
pt_sort);
207 for (
size_t i = 0; i < calib.size(); ++i) {
208 calib.at(i)(*genJets, genEvt->weights.front());
209 calib.at(i)(*recJets, genEvt->weights.front() * recEvt->weights.front(), i);
212 cout <<
"Filling" << endl;
213 if (steering &
DT::fill) tOut->Fill();
216 metainfo.Set<
bool>(
"git",
"complete",
true);
221 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)
- Todo:
- CorrPt?
- 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(
inputs, slice);
48 auto [fOut, tOut] = flow.GetOutput(output);
52 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
53 auto R =
GetR(metainfo);
55 auto genEvt =
isMC ? flow.GetBranchReadOnly<GenEvent>(
"genEvent") :
nullptr;
56 auto recEvt = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
57 auto pu = flow.GetBranchReadOnly<PileUp >(
"pileup" );
58 auto recJets = flow.GetBranchReadWrite<vector<RecJet>>(
"recJets");
60 const auto tables =
config.get<fs::path>(
"corrections.JES.tables");
61 metainfo.Set<fs::path>(
"corrections",
"JES",
"tables", tables);
62 if (!fs::exists(tables))
63 BOOST_THROW_EXCEPTION( fs::filesystem_error(
"File does not exists.",
64 tables, make_error_code(errc::no_such_file_or_directory)) );
66 metainfo.Set<
string>(
"corrections",
"JES",
"campaign",
campaign);
71 ControlPlots raw(
"raw");
72 vector<ControlPlots> calib { ControlPlots(
"nominal") };
74 const bool applySyst = steering &
DT::syst;
78 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
source +
SysDown);
79 metainfo.Set<
string>(
"variations", RecJet::ScaleVar,
source +
SysUp);
86 auto cset = correction::CorrectionSet::from_file(tables.string());
88 string suffix =
"AK" + to_string(R) +
"PF" +
algo;
89 auto nomSF = GetCorrection<correction::CompoundCorrection::Ref>(
102 vector<correction::Correction::Ref> varSFs;
108 varSFs.push_back(GetCorrection<correction::Correction::Ref>(
109 *
cset, short_campaign,
"MC",
source, suffix));
111 for (
DT::Looper looper(tIn); looper(); ++looper) {
115 double weight = recEvt->weights.front();
121 for (
auto& recJet: *recJets) {
124 if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
128 auto corr = nomSF->evaluate({recJet.area, recJet.p4.Eta(), recJet.p4.Pt(), pu->rho});
129 cout << recJet.p4.Pt() <<
' ' << corr <<
'\n';
130 recJet.scales.front() = corr;
132 if (!applySyst)
continue;
136 for (
const correction::Correction::Ref& varSF: varSFs) {
137 auto var = varSF->evaluate({recJet.p4.Eta(), recJet.p4.Pt()});
138 recJet.scales.push_back(corr*(1-var));
139 recJet.scales.push_back(corr*(1+var));
145 for (
size_t i = 0; i < calib.size(); ++i)
146 calib.at(i)(*recJets,
weight, i);
148 sort(recJets->begin(), recJets->end(),
pt_sort);
149 if (steering &
DT::fill) tOut->Fill();
152 metainfo.Set<
bool>(
"git",
"complete",
true);
157 cout << __func__ <<
' ' << slice <<
" end" << endl;
◆ assertValidBinning()
void DAS::JetEnergy::assertValidBinning |
( |
const std::vector< double > & |
v | ) |
|
|
inline |
64 for (
size_t i = 1; i<v.size(); ++i) {
65 if (v[i] > v[i-1])
continue;
66 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);
◆ 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.
92 if (metainfo.
Find(
"flags",
"labels",
"CHS" ))
return "chs";
93 if (metainfo.
Find(
"flags",
"labels",
"PUPPI"))
return "Puppi";
95 std::cerr <<
orange <<
"Couldn't identify CHS or PUPPI. Running default "
96 "CHS for Run 2 (PUPPI foir Run 3)\n" <<
def;
98 auto year = metainfo.
Get<
int>(
"flags",
"year");
99 if (
year > 2019)
return "Puppi";
100 if (
year > 2014)
return "chs";
◆ getBinning()
std::vector<double> DAS::JetEnergy::getBinning |
( |
int |
nBins, |
|
|
float |
first, |
|
|
float |
last |
|
) |
| |
|
inline |
55 std::vector<double>
bins(nBins+1);
56 for(
int i = 0; i <= nBins; ++i)
57 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 |
142 namespace al = boost::algorithm;
143 for (
const auto& correction: corrections) {
144 if (!al::starts_with(correction.first,
campaign) ||
145 !al::contains (correction.first, type) ||
146 !al::contains (correction.first, level) ||
147 !al::ends_with (correction.first, suffix))
150 return correction.second;
153 BOOST_THROW_EXCEPTION(
155 "No `"s +
campaign +
"*"s + type +
"*" + level +
"*"s + suffix
156 +
"` 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 |
81 cout << __func__ <<
' ' << slice <<
" start" << endl;
84 auto tIn = flow.GetInputTree(
inputs, slice);
85 auto [fOut, tOut] = flow.GetOutput(output);
89 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
90 const auto alpha =
config.get<
float>(
"skim.dijet.alpha");
91 metainfo.Set<
float>(
"skim",
"dijet",
"alpha", alpha);
93 auto gEv =
isMC ? flow.GetBranchReadOnly<GenEvent>(
"genEvent"):
nullptr;
94 auto rEv = flow.GetBranchReadOnly<RecEvent>(
"recEvent");
96 auto genJets =
isMC ? flow.GetBranchReadOnly<vector<GenJet>>(
"genJets") :
nullptr;
97 auto recJets = flow.GetBranchReadOnly<vector<RecJet>>(
"recJets");
99 auto genDijet =
isMC ? flow.GetBranchReadOnly<
GenDijet>(
"genDijet") :
nullptr;
100 auto recDijet = flow.GetBranchReadOnly<
RecDijet>(
"recDijet");
106 vector<Balance> recBalances;
107 TDirectory * rec = fOut->mkdir(
"rec");
112 vector<Balance> genBalances;
113 TDirectory * gen =
nullptr;
115 gen = fOut->mkdir(
"gen");
119 TRandom3 r3(metainfo.Seed<39856>(slice));
120 for (
DT::Looper looper(tIn); looper(); ++looper) {
125 auto r = r3.Uniform();
126 int itag = (r < 0.5),
129 for (Balance& bal: recBalances) {
130 auto gw =
isMC ? gEv->Weight(bal.v) : Weight{1.,0},
131 rw = rEv->Weight(bal.v);
134 bal((itag ? *recDijet->first : *recDijet->second),
135 (iprobe ? *recDijet->first : *recDijet->second), gw*rw);
140 for (Balance& bal: genBalances) {
141 auto gw = gEv->Weight(bal.v);
143 bal((itag ? *genDijet->first : *genDijet->second),
144 (iprobe ? *genDijet->first : *genDijet->second), gw);
149 for (
auto&
p: recBalances)
p.Write();
150 for (
auto&
p: genBalances)
p.Write();
152 metainfo.Set<
bool>(
"git",
"complete",
true);
154 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(
inputs, 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");
82 for (
int rhobin = 1; rhobin <=
nRhoBins.at(
year); ++rhobin) {
83 rhobinsResp.push_back(
makeRespHist( Form(
"rhobin%d", rhobin) ,
"JER") );
84 mubinsResp.push_back(
makeRespHist( Form(
"mubin%d", rhobin) ,
"JER") );
88 for (
DT::Looper looper(tIn); looper(); ++looper) {
92 auto recWgt = rEv->weights.front();
93 optional<size_t> irho, imu;
95 static const size_t nrho =
static_cast<size_t>(
nRhoBins.at(
year));
96 for (irho = 0; pileUp->rho >
rho_edges.at(
year).at(*irho+1) && *irho < nrho; ++*irho);
100 static const size_t imuMax = mubinsResp.size()-1;
101 imu =
static_cast<size_t>(pileUp->GetTrPU())/10;
102 *imu = min(*imu,imuMax);
108 matching.match_its.resize(min(matching.match_its.size(),3lu));
110 for (
auto& [rec_it,gen_it]: matching.match_its) {
111 auto ptrec = rec_it->CorrPt();
112 auto ptgen = gen_it->p4.Pt();
116 auto etarec = abs( rec_it->p4.Eta() );
118 inclResp->Fill(ptgen, etarec,
response, recWgt);
119 if (irho) rhobinsResp.at(*irho)->Fill(ptgen, etarec,
response, recWgt);
120 if (imu ) mubinsResp.at(*imu )->Fill(ptgen, etarec,
response, recWgt);
128 inclResp->SetDirectory(fOut);
129 inclResp->SetTitle(
"Response (inclusive)");
130 inclResp->Write(
"Response");
138 TDirectory * d_rho = fOut->mkdir(Form(
"rhobin%d", irho+1), title);
140 rhobinsResp.at(irho)->SetDirectory(fOut);
141 rhobinsResp.at(irho)->SetTitle(
"Response (" + title +
")");
142 rhobinsResp.at(irho)->Write(
"Response");
144 title = Form(
"%d < #mu < %d", irho*10, (irho+1)*10);
146 TDirectory * d_mu = fOut->mkdir(Form(
"mubin%d", irho+1), title);
148 mubinsResp.at(irho)->SetDirectory(fOut);
149 mubinsResp.at(irho)->SetTitle(
"Response (" + title +
")");
150 mubinsResp.at(irho)->Write(
"Response");
153 metainfo.Set<
bool>(
"git",
"complete",
true);
155 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:
- identify algo
- 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` |
40 cout << __func__ <<
" start" << endl;
44 auto R =
config.get<
int>(
"flags.R");
47 cout << R <<
' ' <<
year <<
' ' <<
algo << endl;
49 const auto tag =
config.get<
string>(
"corrections.JES.tag"),
50 level =
config.get<
string>(
"corrections.JES.level");
51 string key = tag +
"_" + level +
"_AK" + to_string(R) +
"PF" +
algo;
54 auto cset = correction::CorrectionSet::from_file(
input.string());
55 correction::Correction::Ref
sf;
59 catch (std::out_of_range& e) {
60 BOOST_THROW_EXCEPTION( std::invalid_argument(
"Nothing corresponding to that key") );
63 cout <<
"Declaring histogram" << endl;
65 map<string, correction::Variable::Type> named_inputs;
66 if (level ==
"L1FastJet" || level ==
"PtResolution") {
67 h = make_unique<TH3D>(level.c_str(), level.c_str(),
72 named_inputs[
"JetA"] = M_PI * pow(R*0.1, 2);
74 else if (level ==
"L2Relative" || level ==
"L2L3Residual")
75 h = make_unique<TH2D>(level.c_str(), level.c_str(),
78 else if (level ==
"ScaleFactor") {
79 h = make_unique<TH1D>(level.c_str(), level.c_str(),
81 named_inputs[
"systematic"] =
"nom";
84 h = make_unique<TH2D>(level.c_str(), level.c_str(),
88 cout <<
"dim = " << h->GetDimension() << endl;
89 for (
int etabin = 1; etabin <= h->GetNbinsX(); ++etabin)
90 for (
int ptbin = 1; ptbin <= h->GetNbinsY(); ++ ptbin)
91 for (
int rhobin = 1; rhobin <= h->GetNbinsZ(); ++rhobin) {
93 named_inputs[
"JetEta"] = h->GetXaxis()->GetBinCenter(etabin);
94 named_inputs[
"JetPt" ] = h->GetYaxis()->GetBinCenter( ptbin);
95 named_inputs[
"Rho" ] = h->GetZaxis()->GetBinCenter(rhobin);
98 vector<correction::Variable::Type>
inputs;
99 for (
const correction::Variable&
input:
sf->inputs())
103 double value =
sf->evaluate(
inputs);
105 cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
107 if (level ==
"L1FastJet")
108 h->SetBinContent(etabin, ptbin, rhobin, value);
109 else if (level ==
"ScaleFactor")
110 h->SetBinContent(etabin, value);
112 h->SetBinContent(etabin, ptbin, value);
116 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
78 std::vector<float> edges;
83 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));
◆ GetR()
Determine the R values for which JetMET provides corrections.
166 auto year = metainfo.
Get<
int>(
"flags",
"year"),
167 R = metainfo.
Get<
int>(
"flags",
"R");
169 std::string warning =
"Not a standard jet size.\n";
171 if (R != 4 && R != 8)
173 return R < 6 ? 4 : 8;
176 if (R != 5 && R != 7)
178 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
109 regex r(
"^(Summer|Fall|Autumn|Winter|Spring)[0-9]{2}[A-Za-z0-9]*");
110 smatch campaign_short;
111 if (!regex_search(
campaign, campaign_short, r))
112 BOOST_THROW_EXCEPTION( invalid_argument(
"The campaign could not be identified"));
113 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 |
52 if (!dijet)
return false;
53 if (
jets.size() < 2)
return false;
56 if (dijet.DeltaPhi() < 2.7 )
return false;
59 if (std::abs(dijet.first ->p4.Eta()) >= 3.0 ||
60 std::abs(dijet.second->p4.Eta()) >= 3.0)
return false;
63 if (
jets.size() > 2) {
64 const auto pt2 =
jets.at(2).CorrPt(v);
65 if (pt2 > alpha * dijet.HT())
return false;
◆ 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 |
123 s <<
"Available corrections:";
124 for (
const auto& correction: corrections) {
125 bool found = correction.first ==
key;
127 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 |
◆ nRhoBins
const std::map<int,int> nRhoBins |
|
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} }
}
◆ rhoBins
const std::map<int,std::vector<TString> > rhoBins |
|
inlinestatic |
Initial value:= {
{ 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);}) }
}
- 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:93
bins
Definition: metPhiCorrectionExample.py:89
static const std::vector< double > y_edges
Definition: binnings.h:36
static const std::map< int, std::vector< double > > rho_edges
Definition: common.h:26
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
FourVector p4
raw four-momentum directly after reconstruction
Definition: PhysicsObject.h:50
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:46
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:17
std::string ScanCorrections(const auto &corrections, const std::string &key="")
Definition: common.h:119
p
Definition: Ntupliser_cfg.py:362
static const std::vector< double > pt_edges
Definition: binnings.h:33
f
Definition: Ntupliser_cfg.py:256
TH1 * CumulativeResponse(TH1 *res1D)
Definition: fitJetResponse.cc:306
std::vector< float > scales
energy scale corrections and variations
Definition: PhysicsObject.h:51
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:106
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:42
static const std::vector< std::string > JES_variations
Definition: common.h:49
genjets
Definition: Ntupliser_cfg.py:272
string jets
Definition: Ntupliser_cfg.py:41
DAS::RecJet recjet
Definition: classes.h:15
DAS::Weight weight
Definition: classes.h:11
config
Definition: Ntupliser_cfg.py:264
static const double maxpt
Definition: getMNobservables.cc:42
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:273
static const int nEtaUncBins
Definition: common.h:47
static const std::vector< double > abseta_edges
JERC binning (taken from JERCProtoLab repo, /macros/common_info/common_binning.hpp)
Definition: common.h:22
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:164
float CorrPt(size_t i=0) const
corrected transverse momentum
Definition: PhysicsObject.h:55
bool passDijetSelection(const auto &dijet, const auto &jets, const Variation &v, const double alpha=0.3)
Definition: getDijetBalance.cc:47
const std::string SysDown
Suffix used for "down" uncertainties. Follows the Combine convention.
Definition: Format.h:10
string isMC
Definition: Ntupliser_cfg.py:59
Di< RecJet, RecJet > RecDijet
Definition: Di.h:80
static const std::map< int, int > nRhoBins
Definition: common.h:32
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:18
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:90
Definition: Variation.h:22
static const int nAbsEtaBins
Definition: common.h:23
static const char * bold
Definition: Step.h:35