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)

Todo:
Hard-coded cut currently tuned for Pythia UL18
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);
78  auto tIn = flow.GetInputTree(inputs, 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  // branches
94  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
95  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
96  auto pileup = flow.GetBranchReadWrite<PileUp>("pileup");
97  auto genJets = flow.GetBranchReadOnly <vector<GenJet>>("genJets");
98  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
99 
100  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),
101  mjj_logw = make_unique<TH2F>("mjj_logw", ";m_{jj}^{rec} [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20);
102 
103  ControlPlots::isMC = isMC;
104  ControlPlots raw("raw"),
105  genveto("genveto"),
106  corrected("corrected");
107 
108  Plots PUbefore("PUbefore"),
109  PUafter("PUafter"),
110  PUgenveto("PUgenveto");
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  raw(*genJets, genEvt->weights.front());
117  raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
118  PUbefore(*genJets, *recJets, *genEvt, *recEvt, *pileup);
119 
120  // 1- removal based on the hard scale
121  auto pthats_it = max_element(pileup->pthats.begin(), pileup->pthats.end());
122  const float maxpthatPU = pthats_it != pileup->pthats.end() ? *pthats_it : 0.;
123  bool bad_rec = (maxpthatPU > genEvt->hard_scale);
124 
125  // 2- removal based on the leading gen pt in the pileup
126  if (!bad_rec)
127  for (size_t i = 0; i < pileup->pthats.size(); ++i) {
128  if (pileup->pthats.at(i) > 0) continue;
129  unsigned long long evtID = pileup->MBevents.at(i);
130  auto it = maxgenpts.find(evtID);
131  if (it == maxgenpts.end()) continue;
132  bad_rec |= it->second > genJets->front().p4.Pt();
133  cout << recEvt->evtNo << ' ' << evtID << ' ' << it->second << endl;
134  if (bad_rec) break;
135  }
136 
137  // 3- removal of reconstruction artefacts happening at excessively high pt
138  if (recJets->size() > 0) {
139  const auto& recjet = recJets->begin();
140  bool matched = false;
141  for (const auto& genjet: *genJets)
142  matched |= DeltaR(genjet.p4, recjet->p4) < maxDR;
143  float maxgenpt = genJets->size() > 0 ? genJets->front().p4.Pt()
144  : genEvt->hard_scale;
145  static const float m = 156;
146  if (recjet->p4.Pt() > max(2*maxgenpt, m))
147  recJets->erase(recjet);
148  }
149 
150  // clean up collections
151  if (bad_rec) {
152  recJets->clear();
153  recEvt->weights *= 0;
154  genveto(*genJets, genEvt->weights.front());
155  genveto(*recJets, genEvt->weights.front() * recEvt->weights.front());
156  PUgenveto(*genJets, *recJets, *genEvt, *recEvt, *pileup);
157  }
158  pileup->pthats.clear();
159  pileup->MBevents.clear();
160 
161  corrected(*genJets, genEvt->weights.front());
162  corrected(*recJets, genEvt->weights.front() * recEvt->weights.front());
163  PUafter(*genJets, *recJets, *genEvt, *recEvt, *pileup);
164 
165  if (steering & DT::fill) tOut->Fill();
166  }
167 
168  metainfo.Set<bool>("git", "complete", true);
169  fOut->cd();
170  raw.Write(fOut);
171  genveto.Write(fOut);
172  corrected.Write(fOut);
173  PUbefore.Write(fOut);
174  PUgenveto.Write(fOut);
175  PUafter.Write(fOut);
176 
177  cout << __func__ << ' ' << slice << " stop" << endl;
178 }

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

Todo:
Hard-coded cut currently tuned for Pythia UL18
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);
78  auto tIn = flow.GetInputTree(inputs, 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  // branches
94  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
95  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
96  auto pileup = flow.GetBranchReadWrite<PileUp>("pileup");
97  auto genJets = flow.GetBranchReadOnly <vector<GenJet>>("genJets");
98  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
99 
100  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),
101  mjj_logw = make_unique<TH2F>("mjj_logw", ";m_{jj}^{rec} [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20);
102 
103  ControlPlots::isMC = isMC;
104  ControlPlots raw("raw"),
105  genveto("genveto"),
106  corrected("corrected");
107 
108  Plots PUbefore("PUbefore"),
109  PUafter("PUafter"),
110  PUgenveto("PUgenveto");
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  raw(*genJets, genEvt->weights.front());
117  raw(*recJets, genEvt->weights.front() * recEvt->weights.front());
118  PUbefore(*genJets, *recJets, *genEvt, *recEvt, *pileup);
119 
120  bool bad_rec = false;
121  float ht = 0;
122  for (const auto& jet : *genJets) {
123  ht += jet.p4.Pt();
124  }
125 
126  // 2- removal based on the leading gen pt in the pileup
127  for (size_t i = 0; i < pileup->pthats.size(); ++i) {
128  unsigned long long evtID = pileup->MBevents.at(i);
129  auto it = maxgenhts.find(evtID);
130  if (it == maxgenhts.end()) continue;
131  bad_rec |= it->second > ht;
132  cout << recEvt->evtNo << ' ' << evtID << ' ' << it->second << endl;
133  if (bad_rec) break;
134  }
135 
136  // 3- removal of reconstruction artefacts happening at excessively high pt
137  if (recJets->size() > 0) {
138  const auto& recjet = recJets->begin();
139  bool matched = false;
140  for (const auto& genjet: *genJets)
141  matched |= DeltaR(genjet.p4, recjet->p4) < maxDR;
142  float maxgenpt = genJets->size() > 0 ? genJets->front().p4.Pt()
143  : genEvt->hard_scale;
144  static const float m = 156;
145  if (recjet->p4.Pt() > max(2*maxgenpt, m))
146  genEvt->weights.front() = 0;
147  }
148 
149  // clean up collections
150  if (bad_rec) {
151  recJets->clear();
152  recEvt->weights *= 0;
153  genveto(*genJets, genEvt->weights.front());
154  genveto(*recJets, genEvt->weights.front() * recEvt->weights.front());
155  PUgenveto(*genJets, *recJets, *genEvt, *recEvt, *pileup);
156  }
157  pileup->pthats.clear();
158  pileup->MBevents.clear();
159 
160  corrected(*genJets, genEvt->weights.front());
161  corrected(*recJets, genEvt->weights.front() * recEvt->weights.front());
162  PUafter(*genJets, *recJets, *genEvt, *recEvt, *pileup);
163 
164  if (steering & DT::fill) tOut->Fill();
165  }
166 
167  metainfo.Set<bool>("git", "complete", true);
168  fOut->cd();
169  raw.Write(fOut);
170  genveto.Write(fOut);
171  corrected.Write(fOut);
172  PUbefore.Write(fOut);
173  PUgenveto.Write(fOut);
174  PUafter.Write(fOut);
175 
176  cout << __func__ << ' ' << slice << " stop" << endl;
177 }

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