DAS  3.0
Das Analysis System
Classes | Functions | Variables
DAS::PUstaub Namespace Reference

Classes

struct  Plots
 

Functions

unordered_map< unsigned long long, float > GetMaxGenPts (const fs::path &input)
 
void applyPUcleaning (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
unordered_map< unsigned long long, float > GetMaxGenHts (const fs::path &input)
 
void applyPUcleaningHT (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 

Variables

static const std::vector< double > pthat_edges = {0, 1e-6, 15, 30, 50, 80, 120, 170, 300, 470, 600, 800, 1000, 1400, 1800, 2400, 3200, 6500}
 
static const int nPtHatBins = pthat_edges.size()-1
 

Function Documentation

◆ applyPUcleaning()

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

Remove PU staub (i.e. high-weight events due to PU sampling)

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
72  {1,0}
73  )
74 {
75  cout << __func__ << ' ' << slice << " start" << endl;
76 
77  DT::Flow flow(steering, inputs);
78  auto tIn = flow.GetInputTree(slice);
79  auto [fOut, tOut] = flow.GetOutput(output);
80 
81  DT::MetaInfo metainfo(tOut);
82  metainfo.Check(config);
83  auto isMC = metainfo.Get<bool>("flags", "isMC");
84  if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
85  make_unique<TFile>(inputs.front().c_str() )) );
86  auto R = metainfo.Get<int>("flags", "R");
87  const auto maxDR = R/10.;
88 
89  auto fMBmaxgenpt = config.get<fs::path>("corrections.PUcleaning.MBmaxgenpt");
90  const auto& maxgenpts = GetMaxGenPts(fMBmaxgenpt);
91  metainfo.Set<fs::path>("corrections", "PUcleaning", "MBmaxgenpt", fMBmaxgenpt);
92 
93  auto artifacts_minpt = config.get<float>("corrections.PUcleaning.artifacts_minpt");
94  metainfo.Set<float>("corrections", "PUcleaning", "artifacts_minpt", artifacts_minpt);
95 
96  // branches
97  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
98  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
99  auto pileup = flow.GetBranchReadWrite<PileUp>("pileup");
100  auto genJets = flow.GetBranchReadOnly <vector<GenJet>>("genJets");
101  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
102 
103  auto pt_logw = make_unique<TH2F>("pt_logw", ";Jet p_{T}^{rec} [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20),
104  mjj_logw = make_unique<TH2F>("mjj_logw", ";m_{jj}^{rec} [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20);
105 
106  ControlPlots::isMC = isMC;
107  ControlPlots raw("raw"),
108  genveto("genveto"),
109  corrected("corrected");
110 
111  Plots PUbefore("PUbefore"),
112  PUafter("PUafter"),
113  PUgenveto("PUgenveto");
114 
115  for (DT::Looper looper(tIn); looper(); ++looper) {
116  [[ maybe_unused ]]
117  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
118 
119  raw(*genJets, genEvt->weights.front());
120  raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
121  PUbefore(*genJets, *recJets, *genEvt, *recEvt, *pileup);
122 
123  // 1- removal based on the hard scale
124  auto pthats_it = max_element(pileup->pthats.begin(), pileup->pthats.end());
125  const float maxpthatPU = pthats_it != pileup->pthats.end() ? *pthats_it : 0.;
126  bool bad_rec = (maxpthatPU > genEvt->hard_scale);
127 
128  // 2- removal based on the leading gen pt in the pileup
129  if (!bad_rec)
130  for (size_t i = 0; i < pileup->pthats.size(); ++i) {
131  if (pileup->pthats.at(i) > 0) continue;
132  unsigned long long evtID = pileup->MBevents.at(i);
133  auto it = maxgenpts.find(evtID);
134  if (it == maxgenpts.end()) continue;
135  bad_rec |= it->second > genJets->front().p4.Pt();
136  cout << recEvt->evtNo << ' ' << evtID << ' ' << it->second << endl;
137  if (bad_rec) break;
138  }
139 
140  // 3- removal of reconstruction artefacts happening at excessively high pt
141  if (recJets->size() > 0 && artifacts_minpt > 0) {
142  const auto& recjet = recJets->begin();
143  bool matched = false;
144  for (const auto& genjet: *genJets)
145  matched |= DeltaR(genjet.p4, recjet->p4) < maxDR;
146  float maxgenpt = genJets->size() > 0 ? genJets->front().p4.Pt()
147  : genEvt->hard_scale;
148  if (recjet->p4.Pt() > max(2*maxgenpt, artifacts_minpt))
149  recJets->erase(recjet);
150  }
151 
152  // clean up collections
153  if (bad_rec) {
154  recJets->clear();
155  recEvt->weights *= 0;
156  genveto(*genJets, genEvt->weights.front());
157  genveto(*recJets, genEvt->weights.front() * recEvt->weights.front());
158  PUgenveto(*genJets, *recJets, *genEvt, *recEvt, *pileup);
159  }
160  pileup->pthats.clear();
161  pileup->MBevents.clear();
162 
163  corrected(*genJets, genEvt->weights.front());
164  corrected(*recJets, genEvt->weights.front() * recEvt->weights.front());
165  PUafter(*genJets, *recJets, *genEvt, *recEvt, *pileup);
166 
167  if (steering & DT::fill) tOut->Fill();
168  }
169 
170  metainfo.Set<bool>("git", "complete", true);
171  fOut->cd();
172  raw.Write(fOut);
173  genveto.Write(fOut);
174  corrected.Write(fOut);
175  PUbefore.Write(fOut);
176  PUgenveto.Write(fOut);
177  PUafter.Write(fOut);
178 
179  cout << __func__ << ' ' << slice << " stop" << endl;
180 }

◆ applyPUcleaningHT()

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

Remove PU staub (i.e. high-weight events due to PU sampling)

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
72  {1,0}
73  )
74 {
75  cout << __func__ << ' ' << slice << " start" << endl;
76 
77  DT::Flow flow(steering, inputs);
78  auto tIn = flow.GetInputTree(slice);
79  auto [fOut, tOut] = flow.GetOutput(output);
80 
81  DT::MetaInfo metainfo(tOut);
82  metainfo.Check(config);
83  auto isMC = metainfo.Get<bool>("flags", "isMC");
84  if (!isMC) BOOST_THROW_EXCEPTION( DE::BadInput("Only MC may be used as input.",
85  make_unique<TFile>(inputs.front().c_str() )) );
86  auto R = metainfo.Get<int>("flags", "R");
87  const auto maxDR = R/10.;
88 
89  auto fMBmaxgenht = config.get<fs::path>("corrections.PUcleaning.MBmaxgenht");
90  const auto& maxgenhts = GetMaxGenHts(fMBmaxgenht);
91  metainfo.Set<fs::path>("corrections", "PUcleaning", "MBmaxgenht", fMBmaxgenht);
92 
93  auto artifacts_minpt = config.get<float>("corrections.PUcleaning.artifacts_minpt");
94  metainfo.Set<float>("corrections", "PUcleaning", "artifacts_minpt", artifacts_minpt);
95 
96  // branches
97  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
98  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
99  auto pileup = flow.GetBranchReadWrite<PileUp>("pileup");
100  auto genJets = flow.GetBranchReadOnly <vector<GenJet>>("genJets");
101  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
102 
103  auto pt_logw = make_unique<TH2F>("pt_logw", ";Jet p_{T}^{rec} [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20),
104  mjj_logw = make_unique<TH2F>("mjj_logw", ";m_{jj}^{rec} [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20);
105 
106  ControlPlots::isMC = isMC;
107  ControlPlots raw("raw"),
108  genveto("genveto"),
109  corrected("corrected");
110 
111  Plots PUbefore("PUbefore"),
112  PUafter("PUafter"),
113  PUgenveto("PUgenveto");
114 
115  for (DT::Looper looper(tIn); looper(); ++looper) {
116  [[ maybe_unused ]]
117  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
118 
119  raw(*genJets, genEvt->weights.front());
120  raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
121  PUbefore(*genJets, *recJets, *genEvt, *recEvt, *pileup);
122 
123  bool bad_rec = false;
124  float ht = 0;
125  for (const auto& jet : *genJets) {
126  ht += jet.p4.Pt();
127  }
128 
129  // 2- removal based on the leading gen pt in the pileup
130  for (size_t i = 0; i < pileup->pthats.size(); ++i) {
131  unsigned long long evtID = pileup->MBevents.at(i);
132  auto it = maxgenhts.find(evtID);
133  if (it == maxgenhts.end()) continue;
134  bad_rec |= it->second > ht;
135  cout << recEvt->evtNo << ' ' << evtID << ' ' << it->second << endl;
136  if (bad_rec) break;
137  }
138 
139  // 3- removal of reconstruction artefacts happening at excessively high pt
140  if (recJets->size() > 0 && artifacts_minpt > 0) {
141  const auto& recjet = recJets->begin();
142  bool matched = false;
143  for (const auto& genjet: *genJets)
144  matched |= DeltaR(genjet.p4, recjet->p4) < maxDR;
145  float maxgenpt = genJets->size() > 0 ? genJets->front().p4.Pt()
146  : genEvt->hard_scale;
147  if (recjet->p4.Pt() > max(2*maxgenpt, artifacts_minpt))
148  genEvt->weights.front() = 0;
149  }
150 
151  // clean up collections
152  if (bad_rec) {
153  recJets->clear();
154  recEvt->weights *= 0;
155  genveto(*genJets, genEvt->weights.front());
156  genveto(*recJets, genEvt->weights.front() * recEvt->weights.front());
157  PUgenveto(*genJets, *recJets, *genEvt, *recEvt, *pileup);
158  }
159  pileup->pthats.clear();
160  pileup->MBevents.clear();
161 
162  corrected(*genJets, genEvt->weights.front());
163  corrected(*recJets, genEvt->weights.front() * recEvt->weights.front());
164  PUafter(*genJets, *recJets, *genEvt, *recEvt, *pileup);
165 
166  if (steering & DT::fill) tOut->Fill();
167  }
168 
169  metainfo.Set<bool>("git", "complete", true);
170  fOut->cd();
171  raw.Write(fOut);
172  genveto.Write(fOut);
173  corrected.Write(fOut);
174  PUbefore.Write(fOut);
175  PUgenveto.Write(fOut);
176  PUafter.Write(fOut);
177 
178  cout << __func__ << ' ' << slice << " stop" << endl;
179 }

◆ GetMaxGenHts()

unordered_map<unsigned long long, float> DAS::PUstaub::GetMaxGenHts ( const fs::path &  input)

Build a map from 2-column file with max gen pt and event IDs from MB sample.

41 {
42  unordered_map<unsigned long long, float> maxgenhts;
43  if (input == "/dev/null")
44  return maxgenhts;
45 
46  if (!fs::exists(input))
47  BOOST_THROW_EXCEPTION(fs::filesystem_error("Missing 2-column file with max gen ht.",
48  input, make_error_code(errc::no_such_file_or_directory)));
49 
50  // use property tree to read the file
51  pt::ptree prop_tree;
52  read_info(input.c_str(), prop_tree);
53 
54  // but convert to unordered map for event loop
55  // --> better optimised
56  for (const auto& evtID_ht: prop_tree) {
57  auto evtID = stoull(evtID_ht.first);
58  auto ht = evtID_ht.second.get_value<float>();
59  maxgenhts[evtID] = ht;
60  }
61 
62  return maxgenhts;
63 }

◆ GetMaxGenPts()

unordered_map<unsigned long long, float> DAS::PUstaub::GetMaxGenPts ( const fs::path &  input)

Build a map from 2-column file with max gen pt and event IDs from MB sample.

41 {
42  unordered_map<unsigned long long, float> maxgenpts;
43  if (input == "/dev/null")
44  return maxgenpts;
45 
46  if (!fs::exists(input))
47  BOOST_THROW_EXCEPTION(fs::filesystem_error("Missing 2-column file with max gen pt.",
48  input, make_error_code(errc::no_such_file_or_directory)));
49 
50  // use property tree to read the file
51  pt::ptree prop_tree;
52  read_info(input.c_str(), prop_tree);
53 
54  // but convert to unordered map for event loop
55  // --> better optimised
56  for (const auto& evtID_pt: prop_tree) {
57  auto evtID = stoull(evtID_pt.first);
58  auto pt = evtID_pt.second.get_value<float>();
59  maxgenpts[evtID] = pt;
60  }
61 
62  return maxgenpts;
63 }

Variable Documentation

◆ nPtHatBins

const int nPtHatBins = pthat_edges.size()-1
static

◆ pthat_edges

const std::vector<double> pthat_edges = {0, 1e-6, 15, 30, 50, 80, 120, 170, 300, 470, 600, 800, 1000, 1400, 1800, 2400, 3200, 6500}
static
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
DAS::PhysicsObject::p4
FourVector p4
raw four-momentum directly after reconstruction
Definition: PhysicsObject.h:50
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:22
DAS::pt_edges
static const std::vector< double > pt_edges
Definition: binnings.h:33
DAS::nMjjBins
static const int nMjjBins
Definition: binnings.h:40
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
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:297
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
DAS::Mjj_edges
static const std::vector< double > Mjj_edges
Definition: binnings.h:34
DAS::PUstaub::GetMaxGenPts
unordered_map< unsigned long long, float > GetMaxGenPts(const fs::path &input)
Build a map from 2-column file with max gen pt and event IDs from MB sample.
Definition: applyPUcleaning.cc:40
DAS::nPtBins
static const int nPtBins
Definition: binnings.h:39
pileup
DAS::PileUp pileup
Definition: classes.h:27
Ntupliser_cfg.isMC
string isMC
Definition: Ntupliser_cfg.py:59
genjet
DAS::GenJet genjet
Definition: classes.h:14
DAS::PUstaub::GetMaxGenHts
unordered_map< unsigned long long, float > GetMaxGenHts(const fs::path &input)
Build a map from 2-column file with max gen pt and event IDs from MB sample.
Definition: applyPUcleaningHT.cc:40
DAS::PileUp::pthats
std::vector< float > pthats
all hard scales found in PU
Definition: Event.h:106
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
DAS::PileUp::MBevents
std::vector< unsigned long long > MBevents
event IDs in MB sample
Definition: Event.h:107
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83