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

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

◆ 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  unique_ptr<TChain> tIn = DT::GetChain(inputs);
55  unique_ptr<TFile> fOut(DT_GetOutput(output));
56  auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
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  RecEvent * recEvt = nullptr;
72  GenEvent * genEvt = nullptr;
73  PileUp * pileup = nullptr;
74  tIn->SetBranchAddress("recEvent", &recEvt);
75  if (isMC)
76  tIn->SetBranchAddress("genEvent", &genEvt);
77  tIn->SetBranchAddress("pileup", &pileup);
78 
79  vector<RecJet> * recJets = nullptr;
80  tIn->SetBranchAddress("recJets", &recJets);
81 
82  TList * ev_wgts = metainfo.List("variations", RecEvent::WeightVar);
83  const size_t NEvWgts = max(1,ev_wgts->GetSize()),
84  NVar = isMC ? NEvWgts : NEvWgts + 2; // the variations obtained with PileUp::GetTrPU()
85  cout << NEvWgts << ' ' << NVar << endl;
86 
87  cout << "Declaring histograms" << endl;
88 
89  vector<PUprofile::Plots> global;
90  TH2D recptLogWgt("recptLogWgt", ";p_{T}^{rec};log(w)", nPtBins, pt_edges.data(), 400, -30, 20);
91  map<int, vector<PUprofile::Plots>> diffs;
93  const vector<string> vars {"nominal", "PU" + SysUp, "PU" + SysDown};
94  for (string var: vars) {
95  auto name = var + "_global";
96  global.push_back(PUprofile::Plots(name));
97  for (auto& turnon: turnons) {
98  name = var + Form("_HLT_PFJet%d", atoi(turnon.first.c_str()));
99  int pt = turnon.second.get_value<int>();
100  diffs[pt].push_back(PUprofile::Plots(name));
101  }
102  name = var + "_ZB";
103  diffs[0].push_back(PUprofile::Plots(name));
104  }
105 
106  PileUp::r3.SetSeed(metainfo.Seed<78923845>(slice));
107 
108  for (DT::Looper looper(tIn, slice); looper(); ++looper) {
109  [[ maybe_unused ]]
110  static auto& cout = (steering & DT::verbose) == DT::verbose ? ::cout : DT::dev_null;
111 
112  auto weight = recEvt->weights.front();
113  if (isMC) weight *= genEvt->weights.front();
114 
115  // find leading pt in tracker acceptance
116  // NOTE: this is biased to high-pt analyses
117  auto leadingInTk = recJets->begin();
118  while (leadingInTk != recJets->end()
119  && leadingInTk->AbsRap() >= 3.0)
120  ++leadingInTk;
121 
122  auto rit = prev(diffs.rend()); // zero-bias trigger by default
123  if (leadingInTk != recJets->end()) { // single jet triggers
124  auto leading_pt = leadingInTk->CorrPt();
125 
126  // find the corresponding trigger
127  rit = diffs.rbegin();
128  while (rit != diffs.rend()) {
129  if (rit->first < leading_pt) break;
130  ++rit;
131  }
132  if (rit == diffs.rend()) { // in case no single-jet trigger has been found, then it should return to ZB case
133  cerr << red << "No single-jet trigger has been found with `leading_pt = " << leading_pt
134  << "` in event " << recEvt->evtNo << ", therefore will consider ZeroBias trigger.\n" << normal;
135  rit = prev(diffs.rend());
136  }
137  }
138 
139  // fill the plots
140  global.at(0)(*pileup, weight, '0', isMC);
141  global.at(1)(*pileup, weight, '+', isMC);
142  global.at(2)(*pileup, weight, '-', isMC);
143  auto& diff = rit->second;
144  diff .at(0)(*pileup, weight, '0', isMC);
145  diff .at(1)(*pileup, weight, '+', isMC);
146  diff .at(2)(*pileup, weight, '-', isMC);
147 
148 
149  if (recJets->size() > 0)
150  recptLogWgt.Fill(recJets->front().CorrPt(), log(weight));
151  }
152 
153  cout << "Writing" << endl;
154  for (size_t i = 0; i < 3; ++i) {
155  fOut->cd();
156  auto var = vars.at(i);
157  auto d = fOut->mkdir(var.c_str());
158  global.at(i).Write(d);
159  for (auto& diff: diffs) {
160  TString name(diff.second.at(i).name);
161  name.ReplaceAll(var + "_", "");
162  auto dd = d->mkdir(name);
163  diff.second.at(i).Write(dd);
164  }
165  }
166 
167  fOut->cd();
168  recptLogWgt.Write("recptLogWgt");
169  tOut->Write();
170 
171  cout << __func__ << ' ' << slice << " stop" << endl;
172 }
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
jmarExample.pt
pt
Definition: jmarExample.py:19
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:30
Step::verbose
static bool verbose
Definition: Step.h:40
DAS::PhysicsObject::p4
FourVector p4
raw four-momentum directly after reconstruction
Definition: PhysicsObject.h:47
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: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::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:65
Darwin::Tools::fill
@ fill
activate -f to fill the tree
Definition: Options.h:28
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:263
normal
static const char * normal
Definition: colours.h:8
DAS::nPtBins
static const int nPtBins
Definition: binnings.h:39
pileup
DAS::PileUp pileup
Definition: classes.h:27
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
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