DAS  3.0
Das Analysis System
Classes | Functions
DAS::PUprofile Namespace Reference

Classes

struct  Correction
 
struct  Plots
 

Functions

void applyBinnedPUprofCorrection (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void applyPUprofCorrection (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void getPUprofile (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 

Function Documentation

◆ applyBinnedPUprofCorrection()

void DAS::PUprofile::applyBinnedPUprofCorrection ( const vector< fs::path >  inputs,
const fs::path  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Apply the PU profile reweighting to the simulation.

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
44  {1,0}
45  )
46 {
47  cout << __func__ << ' ' << slice << " start" << endl;
48 
49  DT::Flow flow(steering);
50  auto tIn = flow.GetInputTree(inputs, slice);
51  auto [fOut, tOut] = flow.GetOutput(output);
52 
53  DT::MetaInfo metainfo(tOut);
54  metainfo.Check(config);
55  auto isMC = metainfo.Get<bool>("flags", "isMC");
56  if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
57  make_unique<TFile>(inputs.front().c_str() )) );
58 
59  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
60  auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
61  auto pu = flow.GetBranchReadOnly<PileUp >("pileup" );
62  auto genJets = flow.GetBranchReadWrite<vector<GenJet>>("genJets");
63  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
64 
65  ControlPlots::isMC = isMC;
66  ControlPlots raw("raw");
67  bool applySyst = steering & DT::syst;
68  vector<ControlPlots> calib { ControlPlots("nominal") };
69  if (applySyst) {
70  metainfo.Set<string>("variations", RecEvent::WeightVar, "PU" + SysUp);
71  metainfo.Set<string>("variations", RecEvent::WeightVar, "PU" + SysDown);
72 
73  calib.push_back(ControlPlots("PU" + SysUp));
74  calib.push_back(ControlPlots("PU" + SysDown));
75  }
76 
77  auto turnon_file = config.get<fs::path>("corrections.normalisation.turnons");
78  auto DataProf = config.get<fs::path>("corrections.PUprofile.Data");
79  auto MCprof = config.get<fs::path>("corrections.PUprofile.MC");
80  for (auto& file: {turnon_file, DataProf, MCprof})
81  if (!fs::exists(file))
82  BOOST_THROW_EXCEPTION( fs::filesystem_error("Input file could not be found",
83  file, make_error_code(errc::no_such_file_or_directory)));
84 
85  pt::ptree turnons;
86  pt::read_info(turnon_file.c_str(), turnons);
87 
88  auto fMCprof = TFile::Open(MCprof.c_str(), "READ"),
89  fDataProf = TFile::Open(DataProf.c_str(), "READ");
90  map<float, PUprofile::Correction> corrections;
91  auto maxWeight = config.get<float>("corrections.PUprofile.maxWeight");
92  for (auto& turnon: turnons) {
93  auto threshold = atoi(turnon.first.c_str());
94  auto trigger = [threshold](const char * variation, bool global) {
95  return global ? Form("%s/intPU" , variation)
96  : Form("%s/HLT_PFJet%d/intPU", variation, threshold);
97  };
98  auto pt = turnon.second.get_value<float>();
99  PUprofile::Correction correction(fMCprof, fDataProf, trigger, maxWeight);
100  corrections.insert( {pt, correction} );
101  }
102  const bool hasZB = fDataProf->GetDirectory("nominal/ZB") && dynamic_cast<TH1*>(fDataProf->Get("nominal/ZB/intPU"))->Integral() > 0;
103  assert(!hasZB);
104  if (hasZB) {
105  auto trigger = [](const char * variation, bool global) {
106  return global ? Form("%s/intPU" , variation)
107  : Form("%s/ZB/intPU", variation);
108  };
109  PUprofile::Correction correction(fMCprof, fDataProf, trigger, maxWeight);
110  corrections.insert( {0, correction} );
111  }
112  fDataProf->Close();
113  fMCprof->Close();
114 
115  metainfo.Set<fs::path>("corrections", "normalisation", "turnons", turnon_file);
116  metainfo.Set<fs::path>("corrections", "PUprofile", "Data", DataProf);
117  metainfo.Set<fs::path>("corrections", "PUprofile", "MC", MCprof);
118  metainfo.Set<float>("corrections", "PUprofile", "maxWeight", maxWeight);
119 
120  PileUp::r3.SetSeed(metainfo.Seed<756347956>(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  int intpu = pu->intpu;
126  auto weight_raw = recEvt->weights.front(); // get current nominal weight before correction
127 
128  // 1) find leading jet pt in tracker acceptance (if it exists)
129  auto pt0 = corrections.begin()->first; // if ZB, then by default ZB, otherwise it will be trig40
130  auto InTk = [](const auto& jet) { return jet.AbsRap() < 3.0; };
131  auto leadRecJetInTk = find_if(recJets->begin(), recJets->end(), InTk);
132  if (leadRecJetInTk != recJets->end())
133  pt0 = leadRecJetInTk->CorrPt();
134  else if (!hasZB) {
135  auto leadGenJetInTk = find_if(genJets->begin(), genJets->end(), InTk);
136  if (leadGenJetInTk != genJets->end())
137  pt0 = leadGenJetInTk->p4.Pt();
138  }
139 
140  auto it = find_if(corrections.rbegin(), prev(corrections.rend()),
141  [&pt0](const auto& th_corr) { return pt0 > th_corr.first; });
142  const auto& correction = it->second;
143  cout << pt0 << ' ' << it->first << endl;
144 
145  // loop over existing weights and correct them with the nominal value
146  auto nominal_correction = correction(intpu, '0');
147  recEvt->weights *= nominal_correction;
148 
149  auto weight_nominal = recEvt->weights.front();
150 
151  // filling at gen level
152  raw (*genJets, genEvt->weights.front());
153  calib.front()(*genJets, genEvt->weights.front());
154 
155  // filling at rec level
156  raw (*recJets, genEvt->weights.front() * weight_raw );
157  calib.front()(*recJets, genEvt->weights.front() * weight_nominal);
158 
159  // systematics
160  if (applySyst) {
161  auto weight_upper = weight_raw*correction(intpu, '+'),
162  weight_lower = weight_raw*correction(intpu, '-');
163  recEvt->weights.push_back(Weight{weight_upper});
164  recEvt->weights.push_back(Weight{weight_lower});
165  calib.at(1)(*genJets, genEvt->weights.front());
166  calib.at(2)(*genJets, genEvt->weights.front());
167  calib.at(1)(*recJets, genEvt->weights.front() * weight_upper);
168  calib.at(2)(*recJets, genEvt->weights.front() * weight_lower);
169  }
170 
171  // filling tree
172  if (steering & DT::fill) tOut->Fill();
173  }
174 
175  auto d = fOut->mkdir("corrections");
176  for (auto& correction: corrections) {
177  d->cd();
178  int threshold = correction.first;
179  auto it = turnons.begin();
180  while (it != turnons.end()) {
181  if (it->second.get_value<int>() == threshold) break;
182  ++it;
183  }
184  //auto it = find_if(turnons.begin(), turnons.end(), [threshold](auto& turnon)
185  // { return threshold == turnon.second.get_value<int>(); });
186  if (it != turnons.end()) {
187  auto dd = d->mkdir(Form("HLT_PFJet%d", atoi(it->first.c_str())));
188  bool reset = slice.second > 0;
189  correction.second.Write(dd, reset);
190  }
191  else { // assuming ZB
192  auto dd = d->mkdir("ZB");
193  bool reset = slice.second > 0;
194  correction.second.Write(dd, reset);
195  }
196  }
197 
198  raw.Write(fOut);
199  for (auto& c: calib)
200  c.Write(fOut);
201 
202  metainfo.Set<bool>("git", "complete", true);
203 
204  cout << __func__ << ' ' << slice << " end" << endl;
205 }

◆ applyPUprofCorrection()

void DAS::PUprofile::applyPUprofCorrection ( const vector< fs::path >  inputs,
const fs::path  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Apply the PU profile reweighting to the simulation.

Todo:
Use raw pt??
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);
47  auto tIn = flow.GetInputTree(inputs, 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  if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
54  make_unique<TFile>(inputs.front().c_str() )) );
55 
56  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
57  auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
58  auto pu = flow.GetBranchReadOnly<PileUp >("pileup" );
59  auto genJets = flow.GetBranchReadWrite<vector<GenJet>>("genJets");
60  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
61 
62  ControlPlots::isMC = isMC;
63  ControlPlots raw("raw");
64  vector<ControlPlots> calib { ControlPlots("nominal") };
65  if (steering & DT::syst) {
66  metainfo.Set<string>("variations", RecEvent::WeightVar, "PU" + SysUp);
67  metainfo.Set<string>("variations", RecEvent::WeightVar, "PU" + SysDown);
68 
69  calib.push_back(ControlPlots("PU" + SysUp));
70  calib.push_back(ControlPlots("PU" + SysDown));
71  }
72 
73  TH2 * recptLogWgt = new TH2D("recptLogWgt", "recptLogWgt;Jet p_{T}^{rec} (GeV);log(w)", nPtBins, pt_edges.data(), 400, -30, 20);
74  TH2 * recptPUwoWeights = new TH2D("recptPUwoWeights", "recptPUwoWeights;Jet p_{T}^{rec} (GeV);intPU", nPtBins, pt_edges.data(), 121, -0.5, 120.5);
75  TH2 * recptPUwWeights = new TH2D("recptPUwWeights", "recptPUwWeights;Jet p_{T}^{rec} (GeV);intPU", nPtBins, pt_edges.data(), 121, -0.5, 120.5);
76 
77  auto DataProf = config.get<fs::path>("corrections.PUprofile.Data");
78  auto MCprof = config.get<fs::path>("corrections.PUprofile.MC");
79  for (auto& file: {DataProf, MCprof})
80  if (!fs::exists(file))
81  BOOST_THROW_EXCEPTION( fs::filesystem_error("Input file could not be found",
82  file, make_error_code(errc::no_such_file_or_directory)));
83  auto fMCprof = TFile::Open(MCprof.c_str(), "READ"),
84  fDataProf = TFile::Open(DataProf.c_str(), "READ");
85  auto global = [](const char * variation, bool) { return Form("%s/intPU", variation); };
86  auto maxWeight = config.get<float>("corrections.PUprofile.maxWeight");
87  PUprofile::Correction correction(fMCprof, fDataProf, global, maxWeight);
88  fDataProf->Close();
89  fMCprof->Close();
90 
91  metainfo.Set<fs::path>("corrections", "PUprofile", "Data", DataProf);
92  metainfo.Set<fs::path>("corrections", "PUprofile", "MC", MCprof);
93  metainfo.Set<float>("corrections", "PUprofile", "maxWeight", maxWeight);
94 
95  PileUp::r3.SetSeed(metainfo.Seed<756347956>(slice));
96  for (DT::Looper looper(tIn); looper(); ++looper) {
97  [[ maybe_unused ]]
98  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
99 
100  int intpu = pu->intpu;
101  auto weight_raw = recEvt->weights.front(); // get current nominal weight before correction
102 
103  // loop over existing weights and correct them with the nominal value
104  auto nominal_correction = correction(intpu, '0');
105  recEvt->weights *= nominal_correction;
106 
107  for (auto& recjet: *recJets) {
108  recptPUwoWeights->Fill(recjet.p4.Pt(), intpu);
109  recptPUwWeights->Fill(recjet.p4.Pt(), intpu, nominal_correction);
110  }
111 
112  for (const auto& recJet: *recJets)
113  recptLogWgt->Fill(recJet.CorrPt(), log(nominal_correction));
114 
115  auto weight_nominal = recEvt->weights.front();
116 
117  // filling at gen level
118  raw (*genJets, genEvt->weights.front());
119  calib.front()(*genJets, genEvt->weights.front());
120 
121  // filling at rec level
122  raw (*recJets, genEvt->weights.front() * weight_raw );
123  calib.front()(*recJets, genEvt->weights.front() * weight_nominal);
124 
125  // systematics
126  if (steering & DT::syst) {
127  auto weight_upper = weight_raw*correction(intpu, '+'),
128  weight_lower = weight_raw*correction(intpu, '-');
129  recEvt->weights.push_back(Weight{weight_upper});
130  recEvt->weights.push_back(Weight{weight_lower});
131  calib.at(1)(*genJets, genEvt->weights.front());
132  calib.at(2)(*genJets, genEvt->weights.front());
133  calib.at(1)(*recJets, genEvt->weights.front() * weight_upper);
134  calib.at(2)(*recJets, genEvt->weights.front() * weight_lower);
135  }
136 
137  // filling tree
138  if (steering & DT::fill) tOut->Fill();
139  }
140 
141  TDirectory * dir = fOut->mkdir("corrections");
142  bool reset = slice.second > 0;
143  correction.Write(dir, reset);
144 
145  raw.Write(fOut);
146  for (auto& c: calib)
147  c.Write(fOut);
148 
149  fOut->cd();
150  recptPUwoWeights->Write("recptPUwoWeights");
151  recptPUwWeights->Write("recptPUwWeights");
152  recptLogWgt->Write("recptLogWgt");
153 
154  metainfo.Set<bool>("git", "complete", true);
155 
156  cout << __func__ << ' ' << slice << " end" << endl;
157 }

◆ getPUprofile()

void DAS::PUprofile::getPUprofile ( const vector< fs::path > &  inputs,
const fs::path &  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Get the PU profile from the n-tuple (data and MC)

Todo:
throw if PU correction already applied
Todo:
also plot without weight
Todo:
Use forward triggers?
Parameters
inputsinput ROOT files (n-tuple)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::Options`
steeringparameters obtained from explicit options
slicenumber and index of slice
49  {1,0}
50  )
51 {
52  cout << __func__ << ' ' << slice << " start" << endl;
53 
54  DT::Flow flow(steering);
55  auto tIn = flow.GetInputTree(inputs, slice);
56  auto [fOut, tOut] = flow.GetOutput(output);
57 
58  DT::MetaInfo metainfo(tOut);
59  metainfo.Check(config);
60  auto isMC = metainfo.Get<bool>("flags", "isMC");
62  //if (!isMC) assert(metainfo.HasCorrection("PUlatest"));
63 
64  auto turnon_file = config.get<fs::path>("corrections.normalisation.turnons");
65  if (!fs::exists(turnon_file))
66  BOOST_THROW_EXCEPTION( fs::filesystem_error("Can't find turnons",
67  turnon_file, make_error_code(errc::no_such_file_or_directory)) );
68  pt::ptree turnons;
69  pt::read_info(turnon_file.c_str(), turnons);
70 
71  auto genEvt = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
72  auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
73  auto pu = flow.GetBranchReadOnly<PileUp >("pileup" );
74  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
75 
76  TList * ev_wgts = metainfo.List("variations", RecEvent::WeightVar);
77  const size_t NEvWgts = max(1,ev_wgts->GetSize()),
78  NVar = isMC ? NEvWgts : NEvWgts + 2; // the variations obtained with PileUp::GetTrPU()
79  cout << NEvWgts << ' ' << NVar << endl;
80 
81  cout << "Declaring histograms" << endl;
82 
83  vector<PUprofile::Plots> global;
84  TH2D recptLogWgt("recptLogWgt", ";p_{T}^{rec};log(w)", nPtBins, pt_edges.data(), 400, -30, 20);
85  map<int, vector<PUprofile::Plots>> diffs;
87  const vector<string> vars {"nominal", "PU" + SysUp, "PU" + SysDown};
88  for (string var: vars) {
89  auto name = var + "_global";
90  global.push_back(PUprofile::Plots(name));
91  for (auto& turnon: turnons) {
92  name = var + Form("_HLT_PFJet%d", atoi(turnon.first.c_str()));
93  int pt = turnon.second.get_value<int>();
94  diffs[pt].push_back(PUprofile::Plots(name));
95  }
96  name = var + "_ZB";
97  diffs[0].push_back(PUprofile::Plots(name));
98  }
99 
100  PileUp::r3.SetSeed(metainfo.Seed<78923845>(slice));
101 
102  for (DT::Looper looper(tIn); looper(); ++looper) {
103  [[ maybe_unused ]]
104  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
105 
106  auto weight = recEvt->weights.front();
107  if (isMC) weight *= genEvt->weights.front();
108 
109  // find leading pt in tracker acceptance
110  // NOTE: this is biased to high-pt analyses
111  auto leadingInTk = recJets->begin();
112  while (leadingInTk != recJets->end()
113  && leadingInTk->AbsRap() >= 3.0)
114  ++leadingInTk;
115 
116  auto rit = prev(diffs.rend()); // zero-bias trigger by default
117  if (leadingInTk != recJets->end()) { // single jet triggers
118  auto leading_pt = leadingInTk->CorrPt();
119 
120  // find the corresponding trigger
121  rit = diffs.rbegin();
122  while (rit != diffs.rend()) {
123  if (rit->first < leading_pt) break;
124  ++rit;
125  }
126  if (rit == diffs.rend()) { // in case no single-jet trigger has been found, then it should return to ZB case
127  cerr << red << "No single-jet trigger has been found with `leading_pt = " << leading_pt
128  << "` in event " << recEvt->evtNo << ", therefore will consider ZeroBias trigger.\n" << normal;
129  rit = prev(diffs.rend());
130  }
131  }
132 
133  // fill the plots
134  global.at(0)(*pu, weight, '0', isMC);
135  global.at(1)(*pu, weight, '+', isMC);
136  global.at(2)(*pu, weight, '-', isMC);
137  auto& diff = rit->second;
138  diff .at(0)(*pu, weight, '0', isMC);
139  diff .at(1)(*pu, weight, '+', isMC);
140  diff .at(2)(*pu, weight, '-', isMC);
141 
142 
143  if (recJets->size() > 0)
144  recptLogWgt.Fill(recJets->front().CorrPt(), log(weight));
145  }
146 
147  cout << "Writing" << endl;
148  for (size_t i = 0; i < 3; ++i) {
149  fOut->cd();
150  auto var = vars.at(i);
151  auto d = fOut->mkdir(var.c_str());
152  global.at(i).Write(d);
153  for (auto& diff: diffs) {
154  TString name(diff.second.at(i).name);
155  name.ReplaceAll(var + "_", "");
156  auto dd = d->mkdir(name);
157  diff.second.at(i).Write(dd);
158  }
159  }
160 
161  fOut->cd();
162  recptLogWgt.Write("recptLogWgt");
163 
164  cout << __func__ << ' ' << slice << " stop" << endl;
165 }
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:93
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:69
Step::verbose
static bool verbose
Definition: Step.h:40
DAS::PhysicsObject::p4
FourVector p4
raw four-momentum directly after reconstruction
Definition: PhysicsObject.h:50
DAS::SysUp
const std::string SysUp
Suffix used for "up" uncertainties. Follows the Combine convention.
Definition: Format.h:8
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:32
DAS::pt_edges
static const std::vector< double > pt_edges
Definition: binnings.h:33
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:29
DAS::Normalisation::reset
void reset(Weights &wgts, const double v=0)
Definition: applyDataNormalisation.h:28
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
recjet
DAS::RecJet recjet
Definition: classes.h:15
weight
DAS::Weight weight
Definition: classes.h:11
Step::red
static const char * red
Definition: Step.h:34
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:264
normal
static const char * normal
Definition: colours.h:8
DAS::nPtBins
static const int nPtBins
Definition: binnings.h:39
DAS::SysDown
const std::string SysDown
Suffix used for "down" uncertainties. Follows the Combine convention.
Definition: Format.h:10
Ntupliser_cfg.isMC
string isMC
Definition: Ntupliser_cfg.py:59
DAS::Normalisation::threshold
const float threshold
Definition: sigmoid.h:14
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83