DAS  3.0
Das Analysis System
DAS::JetEnergy Namespace Reference

Classes

struct  AbstractFit
 
class  Balance
 
struct  ResolutionFit
 
struct  ResponseFit
 

Functions

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::VariationGetScaleVariations (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)
 

Variables

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()
 

Function Documentation

◆ 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
inputsinput ROOT files (n-tuples)
outputname of output ROOT files (n-tuples)
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
sliceslices for running
51  {1,0}
52  )
53 {
54  cout << __func__ << ' ' << slice << " start" << endl;
55 
56  DT::Flow flow(steering, inputs);
57  auto tIn = flow.GetInputTree(slice);
58  auto [fOut, tOut] = flow.GetOutput(output);
59 
60  DT::MetaInfo metainfo(tOut);
61  metainfo.Check(config);
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;
68 
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");
74 
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"); // correction campaign
81  metainfo.Set<string>("corrections", "JER", "campaign", campaign);
82  const auto& core_width = config.get<float>("corrections.JER.core_width");
83 
84  metainfo.Set<float>("corrections", "JER", "core_width", core_width);
85  string algo = GetAlgo(metainfo);
86 
88  ControlPlots raw("raw");
89  vector<ControlPlots> calib { ControlPlots("nominal") };
90 
91  const bool applySyst = steering & DT::syst;
92  vector<string> vars {"nom"}; // keyword for correctionlib
94  if (applySyst) {
95  metainfo.Set<string>("variations", RecJet::ScaleVar, "JER" + SysDown);
96  metainfo.Set<string>("variations", RecJet::ScaleVar, "JER" + SysUp);
97 
98  calib.push_back(ControlPlots("JER" + SysDown));
99  calib.push_back(ControlPlots("JER" + SysUp));
100 
101  vars.push_back("up");
102  vars.push_back("down");
103  }
104 
105  // preparing correctionlib object
106  auto cset = correction::CorrectionSet::from_file(tables.string());
107  string suffix = "AK" + to_string(R) + "PF" + algo;
108 
109  auto SF = GetCorrection<correction::Correction::Ref>(
110  *cset, campaign, "MC", "ScaleFactor", suffix),
111  JER = GetCorrection<correction::Correction::Ref>(
112  *cset, campaign, "MC", "PtResolution", suffix);
113 
114  const auto& smear_tables = config.get<fs::path>("corrections.JER.smear");
115  metainfo.Set<fs::path>("corrections", "JER", "smear", smear_tables);
116 
117  // preparing the JER smearing from JME
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";
123  if (JME_smearing) {
124  auto cset_smear = correction::CorrectionSet::from_file(smear_tables.string());
125  smear = cset_smear->at("JERSmear");
126  }
127 
128  // a few lambda functions
129  auto getSFs = [year,&SF,&vars](const RecJet& recJet, ostream& cout) {
130  vector<float> sfs;
131  sfs.reserve(vars.size());
132  for (const string& var: vars)
133  sfs.push_back(Evaluate(SF, cout, recJet, {}, {}, {}, var));
134  return sfs;
135  };
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;
143 
144  if (applySyst) {
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;
149  }
150  };
151 
152  for (DT::Looper looper(tIn); looper(); ++looper) {
153  [[ maybe_unused ]]
154  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
155  cout << "===============" << endl;
156 
157  raw(*genJets, genEvt->weights.front());
158  raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
159 
160  cout << "Running matching" << endl;
161  JMEmatching matching(*recJets, *genJets);
162  auto fake_its = matching.fake_its; // we extract it because we will modify it
163 
164  cout << "GenJets:";
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();
173  cout << "\nMisses:";
174  for (const auto& genJet: matching.miss_its)
175  cout << '\t' << genJet->p4.Pt();
176  cout << "\nFakes:";
177  for (const auto& recJet: matching.fake_its)
178  cout << '\t' << recJet->CorrPt();
179  cout << endl;
180 
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;
184 
185  float genpt = gen_it->p4.Pt(),
186  recpt = rec_it->CorrPt();
187 
188  float response = recpt / genpt,
189  resolution = Evaluate(JER, cout, *rec_it, pu->rho);;
190 
191  if (abs(response - 1) > core_width * resolution) {
192  cout << "match in the tail" << endl;
193  fake_its.push_back(rec_it); // left for later, with stochastic smearing
194  continue;
195  }
196  cout << "match in the core" << endl;
197 
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;
203  }
204 
205  applySFs(*rec_it, sfs, cout);
206  }
207 
208  cout << "Fakes after identification of matches in tails: ";
209  for (const auto& recJet: fake_its)
210  cout << '\t' << recJet->CorrPt();
211  cout << endl;
212 
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;
217 
218  float resolution = Evaluate(JER, cout, *fake_it, pu->rho);;
219  vector<float> sfs = getSFs(*fake_it, cout);
220 
221  for (float& sf: sfs) {
222  if (JME_smearing) // recJet, rho, genJet, recEvt, syst, JER, JERSF
223  sf = Evaluate(smear, cout, *fake_it, pu->rho, {}, *recEvt, {}, resolution, sf);
224  else {
225  normal_distribution<float> d{0,resolution};
226  static mt19937 m_random_generator(metainfo.Seed<20394242>(slice));
227  float response = d(m_random_generator);
228  sf = max(0.f, 1+response*sqrt(max(0.f, sf*sf-1)));
229  cout << "smearing "
230  << "RecPt=" << fake_it->CorrPt() << ' '
231  << "JER=" << resolution << ' '
232  << "response = " << response << ' '
233  << "-> " << sf << endl;
234  }
235  }
236 
237  applySFs(*fake_it, sfs, cout);
238  }
239 
240  cout << "Sorting by descending pt" << endl;
241  sort(recJets->begin(), recJets->end(), pt_sort);
242 
243  cout << "RecJets after smearing:";
244  for (const auto& recJet: *recJets)
245  cout << '\t' << recJet.CorrPt();
246  cout << endl;
247 
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);
251  }
252 
253  cout << "Filling" << endl;
254  if (steering & DT::fill) tOut->Fill();
255  }
256 
257  metainfo.Set<bool>("git", "complete", true);
258  raw.Write(fOut);
259  for (auto& c: calib)
260  c.Write(fOut);
261 
262  cout << __func__ << ' ' << slice << " end" << endl;
263 }

◆ 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
inputsinput ROOT files (n-tuples)
outputname of output root file containing the histograms
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
sliceslices for running
41  {1,0}
42  )
43 {
44  cout << __func__ << ' ' << slice << " start" << endl;
45 
46  DT::Flow flow(steering, inputs);
47  auto tIn = flow.GetInputTree(slice);
48  auto [fOut, tOut] = flow.GetOutput(output);
49 
50  DT::MetaInfo metainfo(tOut);
51  metainfo.Check(config);
52  auto isMC = metainfo.Get<bool>("flags", "isMC");
53  auto R = GetR(metainfo);
54  auto year = metainfo.Get<int>("flags", "year");
55 
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");
60 
61  const auto tables = config.get<fs::path>("corrections.JES.tables"); // path to directory containing corrections in JSON format
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)) );
66  auto campaign = config.get<string>("corrections.JES.campaign"); // calibration campaign
67  metainfo.Set<string>("corrections", "JES", "campaign", campaign); // e.g. `Summer19UL18_RunA` (everything before "MC" or "DATA")
68  string algo = GetAlgo(metainfo);
69 
72  ControlPlots raw("raw");
73  vector<ControlPlots> calib { ControlPlots("nominal") };
74 
75  const bool applySyst = steering & DT::syst;
76 
77  if (applySyst)
78  for (string source: JES_variations) {
79  metainfo.Set<string>("variations", RecJet::ScaleVar, source + SysDown);
80  metainfo.Set<string>("variations", RecJet::ScaleVar, source + SysUp);
81 
82  calib.push_back(ControlPlots(source + SysDown));
83  calib.push_back(ControlPlots(source + SysUp));
84  }
85 
86  // prepare correctionlib
87  auto cset = correction::CorrectionSet::from_file(tables.string());
88 
89  string suffix = "AK" + to_string(R) + "PF" + algo;
90  auto nomSF = GetCorrection<correction::CompoundCorrection::Ref>(
91  cset->compound(), campaign, isMC ? "MC" : "DATA", "L1L2L3Res", suffix);
92  cout << ScanCorrections(cset->compound(), nomSF->name()) << endl;
93  // Example:
94  // - Summer19UL18_RunA_V6_DATA_L1L2L3Res_AK4PFchs
95  // - Summer19UL18_RunB_V6_DATA_L1L2L3Res_AK4PFchs
96  // - Summer19UL18_RunC_V6_DATA_L1L2L3Res_AK4PFchs
97  // - Summer19UL18_RunD_V6_DATA_L1L2L3Res_AK4PFchs
98  // - Summer19UL18_V5_MC_L1L2L3Res_AK4PFchs
99 
100  if (steering & DT::verbose)
101  cout << ScanCorrections(*cset) << endl;
102 
103  vector<correction::Correction::Ref> varSFs;
104  varSFs.reserve(JES_variations.size());
105 
106  string short_campaign = GetShortCampaign(campaign);
107  if (applySyst)
108  for (string source: JES_variations)
109  varSFs.push_back(GetCorrection<correction::Correction::Ref>(
110  *cset, short_campaign, "MC", source, suffix));
111 
112  for (DT::Looper looper(tIn); looper(); ++looper) {
113  [[ maybe_unused ]]
114  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
115 
116  double weight = recEvt->weights.front();
117  if (isMC) weight *= genEvt->weights.front();
118 
119  // CP before JES
120  raw(*recJets, weight);
121 
122  for (auto& recJet: *recJets) {
123 
124  // sanity check
125  if (recJet.scales.size() != 1 || recJet.scales.front() != 1)
126  BOOST_THROW_EXCEPTION( DE::AnomalousEvent("Unexpected JES corrections", tIn) );
127 
128  // nominal value
129  float corr = Evaluate(nomSF, cout, recJet, pu->rho, {}, *recEvt);
130  recJet.scales.front() = corr;
131 
132  if (!applySyst) continue;
133 
134  // load uncertainties in JECs
135  recJet.scales.reserve(JES_variations.size()*2+1);
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));
140  }
141  }
142 
143  // CP after corrections
144  for (size_t i = 0; i < calib.size(); ++i)
145  calib.at(i)(*recJets, weight, i);
146 
147  sort(recJets->begin(), recJets->end(), pt_sort); // sort again the jets by pt
148  if (steering & DT::fill) tOut->Fill();
149  }
150 
151  metainfo.Set<bool>("git", "complete", true);
152  raw.Write(fOut);
153  for (auto& c: calib)
154  c.Write(fOut);
155 
156  cout << __func__ << ' ' << slice << " end" << endl;
157 }

◆ assertValidBinning()

void DAS::JetEnergy::assertValidBinning ( const std::vector< double > &  v)
inline
64 {
65  assert(v.size() > 1);
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';
69  }
70 }

◆ 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.

307 {
308  double cont = 0;
309  const int NptBins = res1D->GetNbinsX();
310  for (int n = 1; n <= NptBins; ++n) {
311  cont+=res1D->GetBinContent(n);
312  res1D->SetBinContent(n, cont);
313  }
314  return res1D;
315 }

◆ 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
Astd::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
correctionsimple or compound correction
coutoutput stream, possible redirected
recJetfor all corrections
rhoe.g. for smearing
genJetfor stoch smearing
recEvtfor the run or event number
systematice.g. for JER
JERe.g. for stoch smearing
JERSFe.g. for stoch smearing
175  {},
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 = {}
181  )
182 try {
183  using namespace std;
184 
185  cout << correction->name();
186  vector<correction::Variable::Type> inputs;
187  inputs.reserve(10);
188 
189  auto push_back = [&cout,&inputs](const auto& v) {
190  inputs.push_back(v);
191  cout << '=' << v;
192  };
193 
194  for (const correction::Variable& input: correction->inputs()) {
195  const string& n = input.name();
196  cout << ' ' << n;
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) );
214  }
215  auto corr = correction->evaluate(inputs);
216  cout << " -> " << corr << endl;
217  return corr;
218 }
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);
224 }

◆ 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).

79 {
80  const int m = DLog->FindBin(mu);
81 
82  // find maximum on the LHS
83  int imax = m;
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);
87 
88  // find minimum on the RHS
89  int imin = m;
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);
93 
94  return {x1,x2};
95 }

◆ 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
inputinput ROOT file (histograms)
outputname of output ROOT (histograms)
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
234 {
235  cout << __func__ << " start" << endl;
236 
237  auto fIn = make_unique<TFile>(input .c_str(), "READ"),
238  fOut = make_unique<TFile>(output.c_str(), "RECREATE"); // TODO
239  JetEnergy::loopDirsFromFitResponse(fIn.get(), fOut.get(), steering);
240 
241  cout << __func__ << " stop" << endl;
242 }

◆ 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
inputinput ROOT file (histograms)
outputname of output ROOT (histograms)
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
slicenumber and index of slice
517  {1,0}
518  )
519 {
520  cout << __func__ << ' ' << slice << " start" << endl;
521 
522  auto fIn = make_unique<TFile>(input .c_str(), "READ"),
523  fOut = make_unique<TFile>(output.c_str(), "RECREATE");
524  JetEnergy::loopDirsFromGetResponse(fIn.get(), fOut.get(), steering, slice);
525 
526  cout << __func__ << ' ' << slice << " stop" << endl;
527 }

◆ 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

170 {
171  dir->cd();
172  TString stdPrint = h->GetTitle();
173  stdPrint.ReplaceAll(TString( h->GetName() ), "");
174  cout << h->GetName() << '\t' << stdPrint << endl;
175  auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
176 
177  ResolutionFit res(h, cout);
178 
179  res.NSC();
180  res.NSCd();
181 
185 }

◆ 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:

  1. Fit the Gaussian core to get a better estimation of the mean and sigma.
  2. Extend the fit right range, transposing to a Single-Sided Crystal Ball fit, using the mean and sigma better define the fit parameters.
  3. 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.

331 {
332  dir->cd();
333 
334  auto h_ptBins = unique_ptr<TH1>( dynamic_cast<TH1*>(res2D->ProjectionX("h_ptBins", 1, -1)) ); // Only for pt binnning extraction
335  h_ptBins->GetXaxis()->SetTitle("p_{T}^{gen}");
336 
337  // Declaring histograms (control plots)
338  map<TString, unique_ptr<TH1>> hists;
339  for (TString name: {"meanCore", "sigmaCore", // Obtained from DCB fit
340  "kLtail", "kRtail", // DCB tails transitions
341  "nLtail", "nRtail", // DCB tails slope
342  "Lcore", "Rcore", // estimates from derivative
343  "chi2ndf", "Median"}) {
344  // for successful fits
345  hists.insert( {name, unique_ptr<TH1>(dynamic_cast<TH1*>( h_ptBins->Clone(name) )) } );
346  hists[name]->Reset();
347  hists[name]->SetTitle(name);
348  hists[name]->SetDirectory(nullptr);
349 
350  // for failed fits
351  name += "Fail";
352  hists.insert( {name, unique_ptr<TH1>(dynamic_cast<TH1*>( h_ptBins->Clone(name) )) } );
353  hists[name]->Reset();
354  hists[name]->SetTitle(name);
355  hists[name]->SetDirectory(nullptr);
356  }
357 
358  int Npt = res2D->GetXaxis()->GetNbins();
359  for (int ptbin = 1; ptbin <= Npt; ++ptbin) { // Repeat fit procedure for all response distributions
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;
370  static auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
371 
372  auto original = unique_ptr<TH1>( dynamic_cast<TH1*>(res2D->ProjectionY("original", ptbin, ptbin)) ); // Response distribution
373  if (original->GetEntries() < 100) {
374  cout << orange << "Skipped due to low number of events!" << def << endl;
375  continue;
376  }
377 
378  vector<double> edges = GetMergedBinning(original.get());
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));
382  //h->Rebin(2);
383 
384  TDirectory * d_pt = dir->mkdir(Form("ptbin%d", ptbin), title);
385  d_pt->cd();
386  h->SetTitle(title);
387  h->SetDirectory(d_pt);
388  float normalisation = abs( h->Integral() );
389  h->Scale(1.f/normalisation, "width");
390  //cout << "normalisation = " << h->Integral() << endl;
391  h->Write();
392 
393  // Calculating median reponse and Interquartile range (IQR) error
394  // While the mean has standard error estimates, the median uses the
395  // interquartile range to estimate the robustness of the median against outliers or skewed data.
396  // https://en.wikipedia.org/wiki/Interquartile_range
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());
403 
404  TString hlt_title = Form("%.0f < p_{T}^{%s} < %.0f", lowEdge, hlt_level, upEdge);
405 
406  // Calculating Cumulative response histogram and storing it to directory
407  auto hc_tmp = dynamic_cast<TH1*>( h->Clone("hc_tmp") );
408  auto hc = CumulativeResponse(std::move(hc_tmp)); // Cumulative response
409  hc->SetNameTitle("h_cumulative", hlt_title);
410  hc->SetDirectory(d_pt);
411  hc->Write();
412 
413  // the complex fit is split over several steps using a ResponseFit object
414  ResponseFit resp(h, cout);
415 
416  hists.at("Lcore")->SetBinContent(ptbin, resp.p[ResponseFit::KL]);
417  hists.at("Rcore")->SetBinContent(ptbin, resp.p[ResponseFit::KR]);
418 
419  resp.gausCore();
420  resp.SSCB();
421  resp.DSCB();
422 
423  TString suffix = resp.good() ? "" : "Fail";
424 
425  // Save fit parameters into control plots
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]);
430 
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]);
439 
440  hists.at("chi2ndf" + suffix)->SetBinContent(ptbin, resp.chi2ndf .value_or(0));
441  hists.at("chi2ndf" + suffix)->SetBinError (ptbin, resp.chi2ndfErr.value_or(0));
442 
443  hists.at("Median" + suffix)->SetBinContent(ptbin, median);
444  hists.at("Median" + suffix)->SetBinError(ptbin, IQRError);
445  }
446 
447  dir->cd();
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);
454  h.second->Write();
455  }
456 }

◆ GetAlgo()

std::string DAS::JetEnergy::GetAlgo ( const Darwin::Tools::UserInfo metainfo)
inline

Determine from metainfo if CHS or PUPPI has been used to reconstruct the jets.

93 {
94  if (metainfo.Find("flags", "labels", "CHS" )) return "chs";
95  if (metainfo.Find("flags", "labels", "PUPPI")) return "Puppi";
96 
97  std::cerr << orange << "Couldn't identify CHS or PUPPI. Running default "
98  "CHS for Run 2 (PUPPI foir Run 3)\n" << def;
99 
100  auto year = metainfo.Get<int>("flags", "year");
101  if (year > 2019) return "Puppi"; // run 3
102  if (year > 2014) return "chs"; // run 2
103  return ""; // run 1
104 }

◆ getBinning()

std::vector<double> DAS::JetEnergy::getBinning ( int  nBins,
float  first,
float  last 
)
inline
56 {
57  std::vector<double> bins(nBins+1);
58  for(int i = 0; i <= nBins; ++i)
59  bins[i] = first + ((last-first)/nBins) * i;
60  return bins;
61 }

◆ GetBinning()

vector<double> DAS::JetEnergy::GetBinning ( TAxis *  axis)
41 {
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);
47  return bins;
48 }

◆ 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
correctionsfrom `correctionlib`
campaigne.g. Summer19UL18_RunA
typeMC or DATA
levele.g. L1L2L3Res, TimePtEta
suffixe.g. AK4chs
143 {
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))
150  continue;
151 
152  return correction.second;
153  }
154  using namespace std;
155  BOOST_THROW_EXCEPTION(
156  invalid_argument(
157  "No `"s + campaign + "*"s + type + "*" + level + "*"s + suffix
158  + "` correction can be found in the given tables. "s
159  + ScanCorrections(corrections)
160  )
161  );
162 }

◆ 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
inputsinput ROOT files (n-tuples)
outputname of output root file containing the histograms
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
sliceslices for running
79  {1,0}
80  )
81 {
82  cout << __func__ << ' ' << slice << " start" << endl;
83 
84  DT::Flow flow(steering, inputs);
85  auto tIn = flow.GetInputTree(slice);
86  auto [fOut, tOut] = flow.GetOutput(output);
87 
88  DT::MetaInfo metainfo(tOut);
89  metainfo.Check(config);
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);
93 
94  auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent"): nullptr;
95  auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
96 
97  auto genJets = isMC ? flow.GetBranchReadOnly<vector<GenJet>>("genJets") : nullptr;
98  auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
99 
100  auto genDijet = isMC ? flow.GetBranchReadOnly<GenDijet>("genDijet") : nullptr;
101  auto recDijet = flow.GetBranchReadOnly<RecDijet>("recDijet");
102 
103  const bool applySyst = (steering & DT::syst) && !isMC;
104  const auto variations = GetScaleVariations(metainfo, applySyst);
105 
106  // detector level
107  vector<Balance> recBalances;
108  TDirectory * rec = fOut->mkdir("rec");
109  for (const Variation& v: variations)
110  recBalances.emplace_back(rec, v, pt_JERC_edges, abseta_edges);
111 
112  // particle level
113  vector<Balance> genBalances;
114  TDirectory * gen = nullptr;
115  if (isMC) {
116  gen = fOut->mkdir("gen");
117  genBalances.emplace_back(gen, variations.front(), pt_JERC_edges, abseta_edges);
118  }
119 
120  TRandom3 r3(metainfo.Seed<39856>(slice));
121  for (DT::Looper looper(tIn); looper(); ++looper) {
122  [[ maybe_unused ]]
123  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
124 
125  // we take the probe and the tag randomly
126  auto r = r3.Uniform();
127  int itag = (r < 0.5),
128  iprobe = (r >= 0.5);
129 
130  for (Balance& bal: recBalances) {
131  auto gw = isMC ? gEv->Weight(bal.v) : Weight{1.,0},
132  rw = rEv->Weight(bal.v);
133 
134  if (passDijetSelection(*recDijet, *recJets, bal.v, alpha) &&
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);
138  }
139 
140  if (!isMC) continue;
141 
142  for (Balance& bal: genBalances) {
143  auto gw = gEv->Weight(bal.v);
144  if (passDijetSelection(*genDijet, *genJets, bal.v, alpha) &&
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);
148  }
149 
150  } // end of event loop
151 
152  for (auto& p: recBalances) p.Write();
153  for (auto& p: genBalances) p.Write();
154 
155  metainfo.Set<bool>("git", "complete", true);
156 
157  cout << __func__ << ' ' << slice << " stop" << endl;
158 }

◆ 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
inputsinput ROOT files (n-tuples)
outputname of output ROOT file (histograms)
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
sliceslices for running
47  {1,0}
48  )
49 {
50  cout << __func__ << ' ' << slice << " start" << endl;
51 
52  DT::Flow flow(steering, inputs);
53  auto tIn = flow.GetInputTree(slice);
54  auto [fOut, tOut] = flow.GetOutput(output);
55 
56  DT::MetaInfo metainfo(tOut);
57  metainfo.Check(config);
58  auto isMC = metainfo.Get<bool>("flags", "isMC");
59  if (!isMC)
60  BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.", metainfo) );
61 
62  const auto R = metainfo.Get<int>("flags", "R");
63 
64  JMEmatching<>::maxDR = R/10./2.;
65  cout << "Radius for matching: " << JMEmatching<>::maxDR << endl;
66 
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;
70 
71  auto recjets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
72  auto genjets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
73 
74  unique_ptr<TH3> inclResp = makeRespHist("inclusive", "JER"); // Response inclusive in rho
75  vector<unique_ptr<TH3>> rhobinsResp, // Response in bins of rho
76  mubinsResp; // Response in bins of mu (true pileup)
77 
78  const auto year = metainfo.Get<int>("flags", "year");
79  const size_t nRhoBins = GetNRhoBins(year);
80  if (hiPU) {
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") );
86  }
87  }
88 
89  for (DT::Looper looper(tIn); looper(); ++looper) {
90  [[ maybe_unused ]]
91  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
92 
93  auto recWgt = rEv->weights.front();
94  optional<size_t> irho, imu;
95  if (hiPU) {
96  for (irho = 0; pileUp->rho > rho_edges.at(year).at(*irho+1)
97  && *irho < nRhoBins; ++*irho);
98  if (*irho >= nRhoBins)
99  BOOST_THROW_EXCEPTION(
100  DE::AnomalousEvent(Form("irho = %lu > nRhoBins = %zu",
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);
105  }
106 
107  // Matching
108  if (genjets->size() > 3) genjets->resize(3); // Keep only three leading gen jets
109  JMEmatching matching(*recjets, *genjets);
110  matching.match_its.resize(min(matching.match_its.size(),3lu));
111 
112  for (auto& [rec_it,gen_it]: matching.match_its) {
113  auto ptrec = rec_it->CorrPt();
114  auto ptgen = gen_it->p4.Pt();
115 
116  // Fill resolution from smearing
117  auto response = ptrec/ptgen;
118  auto etarec = abs( rec_it->p4.Eta() );
119 
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);
123 
125  } // End of for (pairs) loop
126  } // End of while (event) loop
127 
128  // Creating directory hierarchy and saving histograms
129  fOut->cd();
130  inclResp->SetDirectory(fOut);
131  inclResp->SetTitle("Response (inclusive)");
132  inclResp->Write("Response");
133  if (hiPU)
134  for (int irho = 0; irho < nRhoBins; ++irho) {
135 
137 
138  TString title = Form("%.2f < #rho < %.2f", rho_edges.at(year).at(irho), rho_edges.at(year).at(irho+1));
139  fOut->cd();
140  TDirectory * d_rho = fOut->mkdir(Form("rhobin%d", irho+1), title);
141  d_rho->cd();
142  rhobinsResp.at(irho)->SetDirectory(fOut);
143  rhobinsResp.at(irho)->SetTitle("Response (" + title + ")");
144  rhobinsResp.at(irho)->Write("Response");
145 
146  title = Form("%d < #mu < %d", irho*10, (irho+1)*10);
147  fOut->cd();
148  TDirectory * d_mu = fOut->mkdir(Form("mubin%d", irho+1), title);
149  d_mu->cd();
150  mubinsResp.at(irho)->SetDirectory(fOut);
151  mubinsResp.at(irho)->SetTitle("Response (" + title + ")");
152  mubinsResp.at(irho)->Write("Response");
153  }
154 
155  metainfo.Set<bool>("git", "complete", true);
156 
157  cout << __func__ << ' ' << slice << " end" << endl;
158 }

◆ 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
inputJetMET JES tables in JSON format
outputname of output root file
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
38 {
39  cout << __func__ << " start" << endl;
40 
41  auto fOut = DT::GetOutputFile(output);
42 
43  auto isMC = config.get<bool>("flags.isMC");
44  auto R = config.get<int>("flags.R");
45  auto year = config.get<int>("flags.year");
46  DT::UserInfo userinfo(config);
47  string algo = GetAlgo(userinfo);
48  cout << R << ' ' << year << ' ' << algo << endl;
49 
50  auto campaign = config.get<string>("corrections.JME.campaign"); // calibration campaign
51  auto level = config.get<string>("corrections.JME.level"); // calibration level
52  string suffix = "AK" + to_string(R) + "PF" + algo;
53 
54  auto cset = correction::CorrectionSet::from_file(input.string());
55  auto sf = GetCorrection<correction::Correction::Ref>(
56  *cset, campaign, isMC ? "MC" : "DATA", level, suffix);
58 
59  cout << "Declaring histogram" << endl;
60  unique_ptr<TH1> h;
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(),
64  nAbsEtaBins, abseta_edges .data(),
65  nPtJERCbins, pt_JERC_edges.data(),
66  rho_edges.at(year).size()-1,
67  rho_edges.at(year).data());
68  named_inputs["JetA"] = M_PI * pow(R*0.1, 2); // Jet area estimation
69  named_inputs["JetPhi"] = 0;
70  }
71  else if (level == "L2Relative" || level == "L2L3Residual")
72  h = make_unique<TH2D>(level.c_str(), level.c_str(),
73  nAbsEtaBins, abseta_edges .data(),
74  nPtJERCbins, pt_JERC_edges.data());
75  else if (level == "PtResolution") {
76  // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_PtResolution_AK4PFPuppi (v1)
77  // │ PtResolution for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
78  // │ https://gitlab.cern.ch/cms-jetmet/jerc2json
79  // │ Node counts: Binning: 39, Formula: 305
80  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
81  // │ │ JetEta (real) │
82  // │ │ pseudorapidity of the jet │
83  // │ │ Range: [-5.191, 5.191), overflow ok │
84  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
85  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
86  // │ │ JetPt (real) │
87  // │ │ pT of the jet before specific correction (for JER and uncertainties: after all corrections │
88  // │ │ applied) │
89  // │ │ Range: [-inf, inf), overflow ok │
90  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
91  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
92  // │ │ Rho (real) │
93  // │ │ energy density rho (as measure of PU) │
94  // │ │ Range: [0.0, 52.05), overflow ok │
95  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
96  // │ ╭──── â—€ output ─────╮
97  // │ │ correction (real) │
98  // │ │ No description │
99  // │ ╰───────────────────╯
100  h = make_unique<TH3D>(level.c_str(), level.c_str(),
101  nAbsEtaBins, abseta_edges .data(),
102  nPtJERCbins, pt_JERC_edges.data(),
103  rho_edges.at(year).size()-1,
104  rho_edges.at(year).data());
105  }
106  else if (level == "ScaleFactor") {
107  // ├── 📈 Summer23BPixPrompt23_RunD_JRV1_MC_ScaleFactor_AK4PFPuppi (v1)
108  // │ ScaleFactor for AK4PFPuppi jets, created from Summer23BPixPrompt23_RunD_JRV1_MC by using
109  // │ https://gitlab.cern.ch/cms-jetmet/jerc2json
110  // │ Node counts: Binning: 43, Category: 1470
111  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
112  // │ │ JetEta (real) │
113  // │ │ pseudorapidity of the jet │
114  // │ │ Range: [-5.191, 5.191), overflow ok │
115  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
116  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
117  // │ │ JetPt (real) │
118  // │ │ pT of the jet before specific correction (for JER and uncertainties: after all corrections │
119  // │ │ applied) │
120  // │ │ Range: [10.0, 5373.0), overflow ok │
121  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
122  // │ ╭────────────────────────────────────────── â–¶ input ───────────────────────────────────────────╮
123  // │ │ systematic (string) │
124  // │ │ systematics: nom, up, down │
125  // │ │ Values: down, nom, up │
126  // │ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯
127  // │ ╭──── â—€ output ─────╮
128  // │ │ correction (real) │
129  // │ │ No description │
130  // │ ╰───────────────────╯
131  //
132  if (year < 2020) // run 2
133  h = make_unique<TH1D>(level.c_str(), level.c_str(),
134  nAbsEtaBins, abseta_edges.data());
135  else // run 3
136  h = make_unique<TH2D>(level.c_str(), level.c_str(),
137  nAbsEtaBins, abseta_edges.data(),
138  nPtJERCbins, pt_JERC_edges.data());
139  named_inputs["systematic"] = "nom";
140  }
141  else // assuming JES uncertainties
142  h = make_unique<TH2D>(level.c_str(), level.c_str(),
143  nEtaUncBins, eta_unc_edges.data(),
144  nPtJERCbins, pt_JERC_edges.data());
145 
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) {
150  [[ maybe_unused ]]
151  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
152 
153  // set values by name
154  named_inputs["JetEta"] = h->GetXaxis()->GetBinCenter(etabin);
155  named_inputs["JetPt" ] = h->GetYaxis()->GetBinCenter( ptbin);
156  named_inputs["Rho" ] = h->GetZaxis()->GetBinCenter(rhobin);
157 
158  // order the input value in the necessary order
159  vector<correction::Variable::Type> inputs;
160  for (const correction::Variable& input: sf->inputs())
161  inputs.push_back(named_inputs.at(input.name()));
162 
163  // get and fill the value
164 
165  double value = sf->evaluate(inputs);
166 
167  cout << setw(5) << ptbin << setw(5) << etabin << setw(5) << rhobin << setw(15) << value << endl;
168 
169  if (level == "L1FastJet" || level == "PtResolution")
170  h->SetBinContent(etabin, ptbin, rhobin, value);
171  else if (level == "ScaleFactor")
172  h->SetBinContent(etabin, value);
173  else
174  h->SetBinContent(etabin, ptbin, value);
175  }
176  h->Write();
177 
178  cout << __func__ << " end" << endl;
179 }

◆ 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
75 {
76  assert(maxpt > minpt);
77  assert(minpt > 0);
78  assert(nbins > 1);
79 
80  std::vector<float> edges;
81 
82  float R = pow(maxpt/minpt,1./nbins);
83  for (float pt = minpt; pt <= maxpt; pt *= R) edges.push_back(pt);
85  std::cout << edges.size() << std::endl;
86 
87  return edges;
88 }

◆ 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.

50 {
51  const int N = h->GetNbinsX();
52  vector<double> edges(1,0);
53  edges.reserve(N);
54  int sum = 0;
55  for (int i = 1; i <= N; ++i) {
56  auto content = h->GetBinContent(i);
57  if (content == 0) continue;
58  if (content < 0)
59  cerr << orange << "Found negative bins in "
60  << h->GetName() << def << '\n';
61  auto error = h->GetBinError(i);
62  sum += pow(content/error,2);
63  if (sum < threshold) continue;
64  sum = 0;
65  edges.push_back(h->GetBinLowEdge(i+1));
66  }
67  if (edges.back() != h->GetBinLowEdge(N+1))
68  edges.push_back(h->GetBinLowEdge(N+1));
69 
70  return edges;
71 }

◆ GetNRhoBins()

int DAS::JetEnergy::GetNRhoBins ( int  year)
inline
39 { return rho_edges.at(year).size()-1; }

◆ GetR()

int DAS::JetEnergy::GetR ( const Darwin::Tools::UserInfo metainfo)

Determine the R values for which JetMET provides corrections.

229 {
230  auto year = metainfo.Get<int>("flags", "year"),
231  R = metainfo.Get<int>("flags", "R");
232 
233  std::string warning = "Not a standard jet size.\n";
234  if (year > 2014) { // run 2 and run 3
235  if (R != 4 && R != 8)
236  std::cerr << orange << warning << def;
237  return R < 6 ? 4 : 8;
238  }
239  else { // run 1
240  if (R != 5 && R != 7)
241  std::cerr << orange << warning << def;
242  return R < 6 ? 5 : 7;
243  }
244 }

◆ GetScaleVariations()

vector<DAS::Uncertainties::Variation> DAS::JetEnergy::GetScaleVariations ( const DT::MetaInfo metainfo,
bool  applySyst 
)

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()
32 {
33  vector<DAS::Uncertainties::Variation> variations;
34  variations.emplace_back("nominal","nominal");
35  if (!applySyst) return variations;
36 
37  const TList * groupContents = metainfo.List("variations", DAS::RecJet::ScaleVar);
38  if (!groupContents)
39  BOOST_THROW_EXCEPTION( invalid_argument(groupContents->GetName()) );
40 
41  size_t i = 0;
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, "");
47  s.ReplaceAll(SysDown, "");
48  if (find(begin(JES_variations), end(JES_variations), s) != end(JES_variations))
49  variations.emplace_back(DAS::RecJet::ScaleVar, os->GetString(), ++i);
50  }
51  return variations;
52 }

◆ GetShortCampaign()

std::string DAS::JetEnergy::GetShortCampaign ( const std::string &  campaign)

Extracts for isntance Summer19UL18 from Summer19UL18_RunA

109 {
110  using namespace std;
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();
116 }

◆ 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.

196 {
197  auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
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());
204  ddOut->cd();
205  cout << bold << ddIn->GetPath() << def << endl;
206 
207  loopDirsFromFitResponse(ddIn, ddOut, steering);
208  }
209  else if ( TString( key->ReadObj()->GetName() ).Contains("sigmaCore") ) {
210  auto res = unique_ptr<TH1>( dynamic_cast<TH1*>(key->ReadObj()) );
211  res->SetDirectory(dOut);
212  dOut->cd();
213  res->Write();
214 
215  if ( TString(res->GetName()).Contains("Fail") ) continue;
216 
217  // Fit resolution
218  FitResolution(dOut, std::move(res), steering);
219  }
220  else cout << orange << "Ignoring " << key->ReadObj()->GetName() << " of `"
221  << key->ReadObj()->ClassName() << "` type" << def << endl;
222  } // End of for (list of keys) loop
223  cout << def << flush;
224 }

◆ 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.

468 {
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());
475  ddOut->cd();
476  cout << bold << ddIn->GetPath() << def << endl;
477 
478  loopDirsFromGetResponse(ddIn, ddOut, steering, slice);
479  }
480  else if ( TString( key->ReadObj()->ClassName() ).Contains("TH3") ) {
481  // Assuming for axis, x: pt bins, y: eta bins, z: resolution bins
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) {
487  static int ieta = 0; // Used for parallelisation
488  ++ieta;
489  if (slice.second != ieta % slice.first) continue;
490 
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")) ); // Z vs X, where Z is vertical and X is horizontal axis
497  hist2D->SetTitle( TString( hist2D->GetTitle() ).ReplaceAll(") ", ", " + title + ")") );
498  ddOut->cd();
499  TDirectory * dddOut = ddOut->mkdir(Form("etabin%d", etabin), title);
500  cout << bold << dddOut->GetPath() << def << endl;
501 
502  FitResponseByBin(dddOut, std::move(hist2D), steering);
503  }
504  }
505  else cout << orange << "Ignoring " << key->ReadObj()->GetName() << " of `"
506  << key->ReadObj()->ClassName() << "` type" << def << endl;
507  } // End of for (list of keys) loop
508 }

◆ 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
dijetdijet system
jetsfull jet vector
vsystematic variation
alphareject significant extra jets
52 {
53  if (!dijet) return false;
54  if (jets.size() < 2) return false;
55 
56  // back-to-back
57  if (dijet.DeltaPhi() < 2.7 ) return false;
58 
59  // at least one jet in barrel acceptance
60  if (std::abs(dijet.first ->p4.Eta()) >= 1.3 &&
61  std::abs(dijet.second->p4.Eta()) >= 1.3) return false;
62 
63  // impact of extra jets
64  if (jets.size() > 2) {
65  const auto pt2 = jets.at(2).CorrPt(v);
66  if (pt2 > alpha * dijet.HT()) return false;
67  }
68 
69  return true;
70 }

◆ rhoBins()

std::vector<TString> DAS::JetEnergy::rhoBins ( int  year)
inline
41 {
42  return MakeTitle(rho_edges.at(year), "#rho", false, false,
43  [](double v) { return Form("%.2f", v); });
44 }

◆ round_to()

double DAS::JetEnergy::round_to ( double  value,
double  precision = 1.0 
)
inline

Round number value to precision. Recover "lost" precision.

101 {
102  return round(value / precision) * precision;
103 }

◆ 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
correctionscorrections from `correctionlib`
keykey to highlight in the list
123 {
124  std::stringstream s;
125  s << "Available corrections:";
126  for (const auto& correction: corrections) {
127  bool found = correction.first == key;
128  if (found) s << bold << green;
129  s << ' ' << correction.first;
130  if (found) s << def;
131  }
132  return s.str();
133 }

Variable Documentation

◆ 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} },
}

◆ w

const float w = 0.8
inlinestatic
Todo:
check exact list of variations with JME group
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.name
name
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:48
Ntupliser_cfg.cerr
cerr
Definition: Ntupliser_cfg.py:101
metPhiCorrectionExample.bins
bins
Definition: metPhiCorrectionExample.py:89
DAS::JetEnergy::rho_edges
static const std::map< int, std::vector< double > > rho_edges
Definition: common.h:31
DAS::JetEnergy::loopDirsFromGetResponse
void loopDirsFromGetResponse(TDirectory *dIn, TDirectory *dOut, const int steering, const DT::Slice slice)
Definition: fitJetResponse.cc:464
Darwin::Tools::fill
@ fill
activate -f to fill the tree
Definition: Options.h:27
jmarExample.pt
pt
Definition: jmarExample.py:19
Darwin::Tools::Flow
User-friendly handling of input and output n-tuples.
Definition: Flow.h:78
Step::verbose
static bool verbose
Definition: Step.h:40
Step::def
static const char * def
Definition: Step.h:36
DAS::JetEnergy::Evaluate
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
DAS::RecJet::ScaleVar
static const char *const ScaleVar
Name of jet rec scales in variations.
Definition: Jet.h:38
DAS::JetEnergy::eta_unc_edges
static const std::vector< double > eta_unc_edges
binning for JES uncertainties (taken from JES files)
Definition: common.h:48
DAS::JetEnergy::GetMergedBinning
vector< double > GetMergedBinning(TH1 *h, const int threshold=30)
Definition: fitJetResponse.cc:49
jercExample.algo
string algo
Definition: jercExample.py:60
DAS::GenDijet
Di< GenJet, GenJet > GenDijet
Definition: Di.h:77
Darwin::Tools::UserInfo::Get
T Get(TList *mother, const char *key) const
Definition: UserInfo.h:75
DAS::SysUp
const std::string SysUp
Suffix used for "up" uncertainties. Follows the Combine convention.
Definition: Format.h:8
Ntupliser_cfg.year
int year
Definition: Ntupliser_cfg.py:63
jercExample.cset
cset
Definition: jercExample.py:90
DAS::MakeTitle
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
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.source
source
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:39
DAS::JetEnergy::pt_JERC_edges
static const std::vector< double > pt_JERC_edges
Definition: common.h:22
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:22
DAS::JetEnergy::ScanCorrections
std::string ScanCorrections(const auto &corrections, const std::string &key="")
Definition: common.h:121
Ntupliser_cfg.p
p
Definition: Ntupliser_cfg.py:174
Ntupliser_cfg.f
f
Definition: Ntupliser_cfg.py:309
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:29
DAS::JetEnergy::CumulativeResponse
TH1 * CumulativeResponse(TH1 *res1D)
Definition: fitJetResponse.cc:306
jercExample.key
string key
Definition: jercExample.py:109
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
DAS::JetEnergy::FitResponseByBin
void FitResponseByBin(TDirectory *dir, unique_ptr< TH2 > res2D, const int steering)
Definition: fitJetResponse.cc:328
DAS::JetEnergy::GetShortCampaign
std::string GetShortCampaign(const std::string &campaign)
Extracts for isntance Summer19UL18 from Summer19UL18_RunA
Definition: common.h:108
DAS::Normalisation::makeRespHist
unique_ptr< TH3 > makeRespHist(TString name, TString study)
Create TH3 histograms for the Jet response and HLT response.
Definition: MakeResponseHistos.h:20
DAS::pt_sort
bool pt_sort(const PhysicsObject &j1, const PhysicsObject &j2)
Definition: toolbox.h:18
DAS::MN::minpt
static const double minpt
Definition: getMNobservables.cc:39
DAS::JetEnergy::JES_variations
static const std::vector< std::string > JES_variations
Definition: common.h:51
Ntupliser_cfg.genjets
genjets
Definition: Ntupliser_cfg.py:325
Ntupliser_cfg.jets
string jets
Definition: Ntupliser_cfg.py:41
Darwin::Tools::GetOutputFile
std::shared_ptr< TFile > GetOutputFile(const std::filesystem::path &, const std::source_location=std::source_location::current())
Shortcut to create a reproducible output file (see ROOT Doxygen for details)
Definition: FileUtils.cc:141
weight
DAS::Weight weight
Definition: classes.h:11
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:317
DAS::JetEnergy::GetNRhoBins
int GetNRhoBins(int year)
Definition: common.h:39
DAS::MN::maxpt
static const double maxpt
Definition: getMNobservables.cc:39
DAS::JetEnergy::GetScaleVariations
vector< DAS::Uncertainties::Variation > GetScaleVariations(const DT::MetaInfo &metainfo, bool applySyst)
Definition: Balance.h:30
gitlab_post_comment.response
response
Definition: gitlab_post_comment.py:23
orange
static const char * orange
Definition: colours.h:6
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.input
input
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:35
Darwin::Tools::UserInfo::Find
bool Find(TList *mother, const char *key) const
Recursive finder.
Definition: UserInfo.h:107
Ntupliser_cfg.recjets
recjets
Definition: Ntupliser_cfg.py:326
DAS::JetEnergy::nEtaUncBins
static const int nEtaUncBins
Definition: common.h:49
DAS::JetEnergy::abseta_edges
static const std::vector< double > abseta_edges
JERC binning (taken from JERCProtoLab repo, /macros/common_info/common_binning.hpp)
Definition: common.h:27
Step::green
static const char * green
Definition: Step.h:33
DAS::JetEnergy::GetR
int GetR(const Darwin::Tools::UserInfo &metainfo)
Determine the R values for which JetMET provides corrections.
Definition: common.h:228
DAS::JetEnergy::passDijetSelection
bool passDijetSelection(const auto &dijet, const auto &jets, const Variation &v, const double alpha=0.3)
Definition: getDijetBalance.cc:48
DAS::SysDown
const std::string SysDown
Suffix used for "down" uncertainties. Follows the Combine convention.
Definition: Format.h:10
Darwin::Tools::UserInfo::List
TList * List(TList *mother, const char *key) const
Accest do daughter list from a mother list.
Definition: UserInfo.cc:33
Ntupliser_cfg.isMC
dictionary isMC
Definition: Ntupliser_cfg.py:61
DAS::RecDijet
Di< RecJet, RecJet > RecDijet
Definition: Di.h:80
DAS::Normalisation::threshold
const float threshold
Definition: sigmoid.h:14
jercExample.sf
sf
Definition: jercExample.py:112
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Exceptions::AnomalousEvent
Generic exception for problematic event (during event loop).
Definition: exceptions.h:63
DAS::JetEnergy::loopDirsFromFitResponse
void loopDirsFromFitResponse(TDirectory *dIn, TDirectory *dOut, const int steering)
Definition: fitJetResolution.cc:193
DAS::JetEnergy::nPtJERCbins
static const int nPtJERCbins
Definition: common.h:23
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
generate_html.campaign
campaign
Definition: generate_html.py:86
DAS::JetEnergy::FitResolution
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
Darwin::Tools::UserInfo
Generic meta-information for n-tuple (can be used out of Darwin).
Definition: UserInfo.h:52
DAS::JetEnergy::GetAlgo
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
DAS::Uncertainties::Variation
Definition: Variation.h:22
Ntupliser_cfg.rho
rho
Definition: Ntupliser_cfg.py:320
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83
DAS::JetEnergy::nAbsEtaBins
static const int nAbsEtaBins
Definition: common.h:28
Step::bold
static const char * bold
Definition: Step.h:35