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

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

◆ 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
jmarExample.pt
pt
Definition: jmarExample.py:19
Step::verbose
static bool verbose
Definition: Step.h:40
DAS::PhysicsObject::p4
FourVector p4
raw four-momentum directly after reconstruction
Definition: PhysicsObject.h:47
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:33
Darwin::Tools::GetChain
std::unique_ptr< TChain > GetChain(std::vector< std::filesystem::path > inputs, const char *name="events")
Load chain from a list of files.
Definition: FileUtils.cc:67
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:65
Ntupliser_cfg.genjets
genjets
Definition: Ntupliser_cfg.py:271
Darwin::Tools::fill
@ fill
activate -f to fill the tree
Definition: Options.h:28
recjet
DAS::RecJet recjet
Definition: classes.h:15
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:263
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
Ntupliser_cfg.recjets
recjets
Definition: Ntupliser_cfg.py:272
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
DT_GetOutput
#define DT_GetOutput(output)
Definition: FileUtils.h:96
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:68