DAS  3.0
Das Analysis System
Namespaces | Classes | Typedefs | Enumerations | Functions | Variables
DAS Namespace Reference

Namespaces

 BTagging
 
 Chebyshev
 
 DoubleCrystalBall
 
 JetEnergy
 
 Jets
 
 JetVeto
 
 MissingET
 
 MN
 
 Muon
 
 Muons
 
 Normalisation
 
 Photon
 
 Prefiring
 
 PUprofile
 
 PUstaub
 
 TPS
 
 Uncertainties
 
 Unfolding
 

Classes

class  AbstractEvent
 
struct  AbstractPhysicsObject
 
struct  ControlPlots
 
class  CorrectionlibSFApplier
 
struct  Di
 
class  GenericHistSFApplier
 
struct  GenEvent
 
struct  GenJet
 
struct  GenMuon
 
struct  GenPhoton
 
struct  Greta
 
struct  GretaPolynomial
 
class  GretaPolynomialReader
 
struct  Helper
 
class  JMEmatching
 
struct  MET
 
struct  Parameters
 
class  PhysicsObject
 
struct  PileUp
 
struct  PrimaryVertex
 
struct  RecEvent
 
struct  RecJet
 
struct  RecMuon
 
struct  RecPhoton
 
class  Teddy
 
struct  Thunberg
 
struct  Trigger
 
struct  TriggerLumi
 
struct  Weight
 

Typedefs

using GenDijet = Di< GenJet, GenJet >
 
using GenDimuon = Di< GenMuon, GenMuon >
 
using GenZJet = Di< GenDimuon, GenJet >
 
using GenMmg = Di< GenDimuon, GenPhoton >
 
using RecDijet = Di< RecJet, RecJet >
 
using RecDimuon = Di< RecMuon, RecMuon >
 
using RecZJet = Di< RecDimuon, RecJet >
 
using RecMmg = Di< RecDimuon, RecPhoton >
 
typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< float > > FourVector
 
typedef std::vector< Weight > Weights
 

Enumerations

enum  BinningOptions { default_def = 0b00, logarithmic = 0b01, step = 0b10 }
 

Functions

void mergeNtuples (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
std::vector< TString > MakeTitle (const std::vector< double > &edges, const char *v, bool lowEdge, bool highEdge, std::function< const char *(double)> format)
 
Darwin::Tools::Options Options (const char *, int=Darwin::Tools::none)
 
std::unique_ptr< TH1 > getHistSafe (std::unique_ptr< TFile > &f, const std::string &name)
 
double identity (double x)
 
template<const size_t d>
TF1 * GetSmoothFit (TH1 *h, int ifm, int ifM, int nSigmaStop=0, const char *options="Q0SNRE", bool autoStop=false)
 
template<const size_t d>
TF1 * GetSmoothFit (TH1 *h, double m, double M, int nSigmaStop=0, const char *options="Q0SNRE", bool autoStop=false)
 
bool pt_sort (const PhysicsObject &j1, const PhysicsObject &j2)
 
template<typename T = double>
std::vector< T > MakeBinning (std::pair< T, T > range, double param, int opt=BinningOptions::default_def)
 
template<typename TTreePtr >
bool branchExists (const TTreePtr &tree, TString brName)
 
template<typename T >
std::vector< T * > GetObjects (TDirectory *dir)
 
TDirectory * GetDirectory (TDirectory *dir, const std::vector< const char * > &names)
 
void EventLoop (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void SimpleExec (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering)
 
void getDPSobservable (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
vector< double > GetTriggerTurnons (const fs::path &f)
 
void getMETfraction (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void applyZeroBiasNormalisation (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
std::map< int, TriggerLumi > GetLumiFromFiles (const std::filesystem::path &lumi_file, const std::filesystem::path &turnon_file)
 
map< int, TriggerLumi > GetLumiFromFiles (const fs::path &lumi_file, const fs::path &turnon_file)
 
bool operator== (const PhysicsObject &l, const PhysicsObject &r)
 
bool operator< (const PhysicsObject &l, const PhysicsObject &r)
 
bool operator> (const PhysicsObject &l, const PhysicsObject &r)
 
bool operator== (const Weight &w, const int v)
 
bool operator== (const Weight &w, const float v)
 
bool operator== (const Weight &w, const double v)
 
bool operator== (const Weight &l, const Weight &r)
 
double operator* (const Weight &w, const int v)
 
double operator* (const int v, const Weight &w)
 
double operator* (const Weight &w, const float v)
 
double operator* (const float v, const Weight &w)
 
double operator* (const Weight &w, const double v)
 
double operator* (const double v, const Weight &w)
 
Weight & operator*= (Weight &w, const int v)
 
Weight & operator/= (Weight &w, const int v)
 
Weight & operator*= (Weight &w, const float v)
 
Weight & operator/= (Weight &w, const float v)
 
Weight & operator*= (Weight &w, const double v)
 
Weight & operator/= (Weight &w, const double v)
 
double operator* (const Weight &w1, const Weight &w2)
 
Weights & operator*= (Weights &wgts, const int v)
 
Weights & operator/= (Weights &wgts, const int v)
 
Weights & operator*= (Weights &wgts, const float v)
 
Weights & operator/= (Weights &wgts, const float v)
 
Weights & operator*= (Weights &wgts, const double v)
 
Weights & operator/= (Weights &wgts, const double v)
 
void getHighScalePUeventIDs (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void getHighScalePUeventIDsHT (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
pair< int, int > GetBinRange (TH1 *h)
 
void ScaleWidth (TH1 *h, TH2 *cov)
 
void getSmoothFits (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
Eigen::VectorXd identity (const Eigen::VectorXd &x)
 

Variables

static const std::vector< double > pthat_edges {30,34,39,44,50,56,63,71,80,89,98,108,120,131,143,156,170,196,226,260,300,336,375,420,470,500,531,564,600,645,693,744,800,846,894,946,1000,1088,1183,1287,1400,1491,1587,1690,1800,1934,2078,2233,2400,2579,2771,2978,3200,3578,4000,4472,5000}
 
static const std::vector< double > pt_edges {9, 12, 15,18,21,24,28,32,37,43,49,56,64,74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832,4037}
 
static const std::vector< double > Mjj_edges {160, 200, 249, 306, 372, 449, 539, 641, 756, 887, 1029, 1187, 1361, 1556, 1769, 2008, 2273, 2572, 2915, 3306, 3754, 4244, 4805, 5374, 6094, 6908, 7861, 8929, 10050 }
 
static const std::vector< double > pt_av_edges {147, 175, 207, 243, 284, 329, 380, 437, 499, 569, 646, 732, 827, 931, 1046, 1171, 1307, 1458, 1621, 1806, 2003, 2217, 2453, 2702}
 
static const std::vector< double > y_edges {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}
 
static const std::vector< double > n_edges {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}
 
static const int nYbins = y_edges.size()-1
 
static const int nPtBins = pt_edges.size()-1
 
static const int nMjjBins = Mjj_edges.size()-1
 
static const int nPtDijetBins = pt_av_edges.size()-1
 
static const int maxMult = n_edges.size()-1
 
static const std::vector< TString > yBins = MakeTitle(y_edges, "|y|", false, true, [](double v) { return Form("%.1f", v);} )
 
static const double ptmin = pt_edges.front()
 
static const double ptmax = pt_edges.back()
 
static const double ymin = y_edges.front()
 
static const double ymax = y_edges.back()
 
static const int nPUbins = 100
 
static const std::vector< double > binsY = {0., 0.25, .5, 0.75, 1., 1.25, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 4.25, 4.5, 4.75, 5., 5.25, 5.5, 5.75, 6., 6.25, 6.5, 6.75, 7., 7.25, 7.5, 7.75, 8., 8.25, 8.5, 8.75, 9., 9.25, 9.5}
 
static const std::vector< double > binsPhi = {0.,.1*M_PI, .2*M_PI, .3*M_PI, .4*M_PI, .5*M_PI, .6*M_PI, .7*M_PI, .8*M_PI, .9*M_PI, M_PI}
 
const std::string SysUp = "Up"
 
const std::string SysDown = "Down"
 

Typedef Documentation

◆ FourVector

typedef ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< float > > FourVector

◆ GenDijet

using GenDijet = Di<GenJet,GenJet>

◆ GenDimuon

using GenDimuon = Di<GenMuon,GenMuon>

◆ GenMmg

using GenMmg = Di<GenDimuon,GenPhoton>

◆ GenZJet

using GenZJet = Di<GenDimuon,GenJet>

◆ RecDijet

using RecDijet = Di<RecJet,RecJet>

◆ RecDimuon

using RecDimuon = Di<RecMuon,RecMuon>

◆ RecMmg

using RecMmg = Di<RecDimuon,RecPhoton>

◆ RecZJet

using RecZJet = Di<RecDimuon,RecJet>

◆ Weights

typedef std::vector<Weight> Weights

Enumeration Type Documentation

◆ BinningOptions

enum BinningOptions
Enumerator
default_def 

linear with param as number of bins

logarithmic 

log with param as number of bins

step 

linear with param as step for binning

23  {
24  default_def = 0b00,
25  logarithmic = 0b01,
26  step = 0b10
27 };

Function Documentation

◆ applyZeroBiasNormalisation()

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

Normalise with effective luminosity or prescales. The weight $\frac{ \text{Prescale} * \mathcal{L}^{-1} }{ \epsilon_\text{trigger} } $ is applied by event. Apply phaseSelection and remove events already covered by jet triggers.

To combine ZeroBias and JetHT together, we take the first trigger's TurnOn point(lowest pt threshold in $DARWIN_TABLES/triggers/lumi) as maxpt.

Parameters
inputsinput ROOT files (n-tuples)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
45  {1,0}
46  )
47 {
48  cout << __func__ << ' ' << slice << " start" << endl;
49 
50  DT::Flow flow(steering, inputs);
51  auto tIn = flow.GetInputTree(slice);
52  auto [fOut, tOut] = flow.GetOutput(output);
53 
54  DT::MetaInfo metainfo(tOut);
55  metainfo.Check(config);
56  auto isMC = metainfo.Get<bool>("flags", "isMC");
57  if (isMC)
58  BOOST_THROW_EXCEPTION( DE::BadInput("Only real data may be used as input.", metainfo) );
59 
60  auto total_lumi = config.get<float>("corrections.normalisation.total_lumi");
61  auto maxpt = config.get<float>("corrections.normalisation.maxpt");
62  auto method = config.get<string>("corrections.normalisation.method" );
63 
64  auto inv_total_lumi = 1. / total_lumi;
65 
66  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets", DT::facultative);
67  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
68  auto trigger = flow.GetBranchReadOnly<Trigger>("zbTrigger");
69 
70  for (DT::Looper looper(tIn); looper(); ++looper) {
71  [[ maybe_unused]]
72  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
73 
74  // we find the leading jet in tracker acceptance
75  auto leadingInTk = DAS::Normalisation::phaseSel(*recJets);
76 
77  // removed events alreadu covered by jet triggers
78  if (leadingInTk != recJets->end() && leadingInTk->CorrPt() >= maxpt) continue;
79 
80  auto preHLT = trigger->PreHLT.front();
81  auto preL1min = trigger->PreL1min.front();
82  auto preL1max = trigger->PreL1max.front();
83  if (preL1min != preL1max)
84  BOOST_THROW_EXCEPTION( invalid_argument("Different preL1min and preL1max value.") );
85 
86  auto prescale = preHLT * preL1min;
87 
88  // sanity checks:
89  // we assume that the same prescale is indeed constant for a given LS
90  {
91  static vector<map<pair<int,int>,int>> prescales(1); // this could be written in a simpler way...
92  pair<int,int> runlum = {recEvt->runNo, recEvt->lumi};
93  if (prescales.front().count(runlum))
94  assert(prescales.front().at(runlum) == prescale);
95  else
96  prescales.front()[runlum] = prescale;
97  }
98 
99  if (method == "prescales")
100  recEvt->weights *= prescale * inv_total_lumi * 1;
101  else if (method == "lumi")
102  recEvt->weights *= inv_total_lumi;
103  else {
104  BOOST_THROW_EXCEPTION( invalid_argument("method: " + method + " is not defined.") );
105  }
106 
107  if (steering & DT::fill) tOut->Fill();
108  }
109 
110  metainfo.Set<bool>("git", "complete", true);
111 
112  cout << __func__ << ' ' << slice << " stop" << endl;
113 }

◆ branchExists()

bool DAS::branchExists ( const TTreePtr &  tree,
TString  brName 
)

Check if branch exists

Loop over branches of tree and check if brName exists

75 {
76  auto brList = tree->GetListOfBranches();
77  bool brFound = false;
78  for (auto it = brList->begin(); it != brList->end(); ++it) {
79  TString name = (*it)->GetName();
80  if (name == brName) brFound = true;
81  }
82  return brFound;
83 }

◆ EventLoop()

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

Template for function (TODO)

Parameters
inputsinput ROOT files (n-tuples)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
39  {1,0}
40  )
41 {
42  cout << __func__ << ' ' << slice << " start" << endl;
43 
44  DT::Flow flow(steering, inputs);
45  auto tIn = flow.GetInputTree(slice);
46  auto [fOut, tOut] = flow.GetOutput(output);
47 
48  DT::MetaInfo metainfo(tOut);
49  metainfo.Check(config);
50  auto isMC = metainfo.Get<bool>("flags", "isMC");
51  // TODO: fetch info from MetaInfo and/or fill it with new information
52 
53  auto genJets = isMC ? flow.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
54  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets", DT::facultative);
55  // TODO: load any branch you may need
56 
57  for (DT::Looper looper(tIn); looper(); ++looper) {
58  [[ maybe_unused]]
59  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
60 
61  // TODO: modify the branches
62 
63  if (steering & DT::fill) tOut->Fill();
64  }
65 
66  metainfo.Set<bool>("git", "complete", true);
67 
68  cout << __func__ << ' ' << slice << " stop" << endl;
69 }

◆ GetBinRange()

pair<int,int> DAS::GetBinRange ( TH1 *  h)

Find range for non-empty bins in an input histogram.

33 {
34  int N = h->GetNbinsX();
35  int im, iM;
36  for (im = 1; h->GetBinContent(im) <= 0 && im <= N; ++im);
37  for (iM = im; h->GetBinContent(iM) > 0 && iM <= N; ++iM);
38  --iM;
39  return {im,iM};
40 }

◆ GetDirectory()

TDirectory* DAS::GetDirectory ( TDirectory *  dir,
const std::vector< const char * > &  names 
)
inline

Get ((...)sub)subdirectory of dir.

103 {
104  for (const char * name: names) {
105  TDirectory * subdir = dir->GetDirectory(name);
106  if (!subdir)
107  BOOST_THROW_EXCEPTION( std::invalid_argument(
108  Form("`%s` could not be found in `%s`", name, dir->GetName())) );
109  dir = subdir;
110  }
111  return dir;
112 }

◆ getDPSobservable()

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

Calculates and histograms DPS observables from mixed events.

Parameters
inputsinput ROOT file(s) or directory
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
40  {1,0}
41  )
42 {
43  cout << __func__ << " " << slice << " start" << endl;
44 
45  DT::Flow flow(steering, inputs);
46  auto tIn = flow.GetInputTree("events");
47 
48  auto [fOut, tOut] = flow.GetOutput(output);
49  DT::MetaInfo metainfo(tOut);
50  metainfo.Check(config);
51  auto isMC = metainfo.Get<bool>("flags", "isMC");
52 
53  // --- Declare histograms for classic DPS observables (inspired by PhysRevD.93.054019) ---
54  TH1F* h_delta_phi = new TH1F("h_delta_phi", "Azimuthal Angle between Dimuon and Leading Jet;#Delta#phi(Z, j_{1});Events", 10, 0, M_PI);
55  TH1F* h_delta_pt_rel = new TH1F("h_delta_pt_rel", "Relative pT Balance;#Delta p_{T}^{rel}(j_{1}, j_{2});Events", 10, 0, 1);
56  TH1F* h_delta_S = new TH1F("h_delta_S", "Azimuthal Separation between Z and Dijet;#Delta S (Z, dijet);Events", 10, 0, M_PI);
57 
58  // --- Declare histograms for high-purity MPI studies (inspired by PhysRevLett.132.041901) ---
59  // For "Average Leading Jet pT vs. Z pT" plot
60  TProfile* h_avg_jet_pt_vs_z_pt = new TProfile("h_avg_jet_pt_vs_z_pt", "Average Leading Jet pT vs. Z pT;<p_{T,Z}> [GeV];<p_{T,jet1}> [GeV]", 50, 0, 100);
61 
62  // For "Inclusive Jet Spectrum" plot
63  TH1F* h_jet_pt_spectrum = new TH1F("h_jet_pt_spectrum",
64  "Jet pT Spectrum (for cumulative studies);p_{T,jet} [GeV];Jets", 200, 0, 1000);
65 
66  // --- Branches to read from the n-tuple ---
67  auto recEvent = flow.GetBranchReadOnly<RecEvent>("recEvent");
68  auto genEvent = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
69  auto recMuons = flow.GetBranchReadOnly<vector<RecMuon>>("recMuons");
70  auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
71 
72  // --- Event Loop ---
73  for (DT::Looper looper(tIn); looper(); ++looper) {
74 
75  // Basic object selection
76  if (recMuons->size() < 2) {
77  continue;
78  }
79 
80  // Reconstruct the Z boson from the two leading muons
81  auto dimuon = make_di(recMuons->at(0), recMuons->at(1));
82 
83  // --- Fill plots for classic DPS observables ---
84  if (recJets->size() >= 2) {
85  const auto& jet1 = recJets->at(0);
86  const auto& jet2 = recJets->at(1);
87  auto dijet = make_di(jet1, jet2);
88 
89  double w = recEvent->Weight() * dimuon.Weight() * dijet.Weight();
90  if (isMC) w *= genEvent->Weight();
91 
92  // 1. Delta Phi (dimuon, leading jet)
93  double delta_phi = abs(ROOT::Math::VectorUtil::DeltaPhi(dimuon.CorrP4(), jet1.CorrP4()));
94  h_delta_phi->Fill(delta_phi, w);
95 
96  // 2. Delta pT rel (jet1, jet2)
97  double delta_pt_rel = dijet.CorrP4().Pt() / (jet1.CorrPt() + jet2.CorrPt());
98  h_delta_pt_rel->Fill(delta_pt_rel, w);
99 
100  // 3. Delta S (Z, dijet)
101  if (dijet.CorrPt() > 0) {
102  double delta_S = acos((dimuon.CorrP4().Vect().Dot(dijet.CorrP4().Vect())) / (dimuon.CorrPt() * dijet.CorrPt()));
103  h_delta_S->Fill(delta_S, w);
104  }
105  }
106 
107  // --- Selection cuts for high-purity MPI studies (PhysRevLett.132.041901) ---
108  bool passes_common_cuts = (dimuon.CorrP4().M() > 66.0 && dimuon.CorrP4().M() < 116.0 &&
109  recMuons->at(0).CorrPt() > 27.0 && recMuons->at(1).CorrPt() > 27.0 &&
110  recMuons->at(0).AbsRap() < 2.4 &&
111  recMuons->at(1).AbsRap() < 2.4);
112 
113  if (!passes_common_cuts) {
114  continue;
115  }
116 
117  double w_common = recEvent->Weight() * dimuon.Weight();
118  if (isMC) w_common *= genEvent->Weight();
119 
120  // --- Fill "Average Leading Jet pT vs. Z pT" plot ---
121  if (!recJets->empty()) {
122  const auto& leadingJet = recJets->at(0);
123  h_avg_jet_pt_vs_z_pt->Fill(dimuon.CorrPt(), leadingJet.CorrPt(), w_common * leadingJet.Weight());
124  }
125 
126  // --- Fill temporary histogram for "Inclusive Jet Spectrum" ---
127  if (dimuon.CorrPt() < 2.0) {
128  for (const auto& jet : *recJets | views::filter([](const RecJet& x) { return x.AbsRap() < 2.0; })) {
129  h_jet_pt_spectrum->Fill(jet.CorrPt(), w_common * jet.Weight());
130  }
131  }
132  }
133 
134  // --- Save Histograms ---
135  fOut->cd();
136  h_delta_phi->Write();
137  h_delta_pt_rel->Write();
138  h_delta_S->Write();
139  h_avg_jet_pt_vs_z_pt->Write();
140  h_jet_pt_spectrum->Write();
141 
142  cout << __func__ << " " << slice << " stop" << endl;
143 }

◆ getHighScalePUeventIDs()

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

Make a 2-column list of MB event IDs with the max genpt (if any).

Unlike most commands, this command should be run on the output of the n-tupliser directly, and only makes sense (even though it would technically work) for MB samples without pileup.

Todo:
include test in CI
Parameters
inputsinput ROOT files (n-tuples)
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, inputs);
55  auto tIn = flow.GetInputTree(slice, "ntupliser/events");
56 
57  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
58  auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
59  auto genJets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
60 
61  const auto minpt = config.get<float>("corrections.PUstaub.minpt");
62 
63  unordered_map<unsigned long long, float> maxgenpts;
64 
65  for (DT::Looper looper(tIn); looper(); ++looper) {
66  [[ maybe_unused]]
67  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
68 
69  if (genJets->size() == 0) continue;
70  float pthat = genEvt->hard_scale,
71  maxgenpt = genJets->front().p4.Pt();
72 
73  if (pthat > 0 ? (maxgenpt > 3*pthat) // refinement to remove crazy parton shower
74  : (maxgenpt > minpt)) // refinment to remove diffractive events
75  maxgenpts[recEvt->evtNo] = maxgenpt;
76  }
77 
78  pt::ptree prop_tree;
79  for (const auto& evtID_pt: maxgenpts) {
80  auto evtID = Form("%llu", evtID_pt.first);
81  auto pt = evtID_pt.second;
82  prop_tree.put<float>(evtID, pt);
83  }
84  write_info(output.c_str(), prop_tree);
85 
86  cout << __func__ << ' ' << slice << " stop" << endl;
87 }

◆ getHighScalePUeventIDsHT()

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

Make a 2-column list of MB event IDs with the max genht (if any).

Unlike most commands, this command should be run on the output of the n-tupliser directly, and only makes sense (even though it would technically work) for MB samples without pileup.

Todo:
include test in CI
Parameters
inputsinput ROOT files (n-tuples)
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, inputs);
55  auto tIn = flow.GetInputTree(slice, "ntupliser/events");
56 
57  auto genEvt = flow.GetBranchReadOnly<GenEvent>("genEvent");
58  auto recEvt = flow.GetBranchReadOnly<RecEvent>("recEvent");
59  auto genJets = flow.GetBranchReadOnly<vector<GenJet>>("genJets");
60 
61  const auto minht = config.get<float>("corrections.PUstaub.minht");
62 
63  unordered_map<unsigned long long, float> hts;
64 
65  for (DT::Looper looper(tIn); looper(); ++looper) {
66  [[ maybe_unused]]
67  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
68 
69  if (genJets->size() == 0) continue;
70 
71  float ht = 0;
72  for (const auto& jet : *genJets) {
73  ht += jet.p4.Pt();
74  }
75 
76  if (ht > minht)
77  hts[recEvt->evtNo] = ht;
78 
79  }
80 
81  pt::ptree prop_tree;
82  for (const auto& evtID_ht: hts) {
83  auto evtID = Form("%llu", evtID_ht.first);
84  auto ht = evtID_ht.second;
85  prop_tree.put<float>(evtID, ht);
86  }
87  write_info(output.c_str(), prop_tree);
88 
89  cout << __func__ << ' ' << slice << " stop" << endl;
90 }

◆ getHistSafe()

std::unique_ptr<TH1> DAS::getHistSafe ( std::unique_ptr< TFile > &  f,
const std::string &  name 
)
inline

Gets an histogram from a TFile.

Exceptions
Darwin::Exceptions::BadInputif the histogram doesn't exist
Todo:
Merge with Flow
26 {
28  auto h = std::unique_ptr<TH1>(f->Get<TH1>(name.c_str()));
29  if (!h) {
30  std::cerr << name << std::endl;
31  BOOST_THROW_EXCEPTION( Darwin::Exceptions::BadInput("Can't find required histogram", f) );
32  }
33  h->SetDirectory(nullptr);
34  return h;
35 }

◆ GetLumiFromFiles() [1/2]

map<int, TriggerLumi> DAS::GetLumiFromFiles ( const fs::path &  lumi_file,
const fs::path &  turnon_file 
)

GetLumiFromFiles returns a map containing key value pairs corresponding to the trigger threshold at HLT (e.g. 40, 60, etc.) and the object containing the turn-on point and its effective luminosity.

37 {
38  map<int, TriggerLumi> trigger_luminosities;
39 
40  ifstream infile; // will be used twice consecutively
41 
42  auto checkFile = [&infile]() {
43  if (!infile.fail()) return;
44  cerr << "Failure at opening the file\n";
45  exit(EXIT_FAILURE);
46  };
47 
48  // 1) get effective luminosity
49  assert(fs::exists(lumi_file));
50  infile.open(lumi_file);
51  cout << "Opening lumi file: " << lumi_file << endl;
52  checkFile();
53  while (infile.good()) {
54  int trigger;
55  float lumi;
56  infile >> trigger >> lumi;
57  if (infile.eof()) break;
58  TriggerLumi trigger_lumi(trigger);
59  trigger_lumi.weight = 1./lumi;
60  cout << trigger << ' ' << lumi << ' ' << trigger_lumi.weight << '\n';
61  trigger_luminosities.insert({trigger, trigger_lumi});
62  }
63  infile.close();
64 
65  // 2) get turn-on points
66  cout << "Opening turn-on file: " << turnon_file << endl;
67  assert(fs::exists(turnon_file));
68  infile.open(turnon_file);
69  checkFile();
70  while (infile.good()) {
71  int trigger;
72  double turnon;
73  infile >> trigger >> turnon;
74  if (infile.eof()) break;
75  cout << trigger << ' ' << turnon << '\n';
76  TriggerLumi& trigger_lumi = trigger_luminosities.at(trigger);
77  trigger_lumi.turnon = ceil(turnon); //goes to int
78  }
79  infile.close();
80 
81  // then remove pairs with dummy entry
82  for (auto trigger_lumi = trigger_luminosities.begin();
83  trigger_lumi != trigger_luminosities.end(); /* nothing */) {
84 
85  // remove or increment
86  if ( trigger_lumi->second.turnon == DUMMY
87  || trigger_lumi->second.weight == DUMMY)
88  trigger_lumi = trigger_luminosities.erase(trigger_lumi);
89  else ++trigger_lumi;
90 
91  }
92 
93  // display
94  for (auto trigger_lumi: trigger_luminosities)
95  cout << trigger_lumi.first << '\t'
96  << trigger_lumi.second.turnon << '\t'
97  << trigger_lumi.second.weight << '\n';
98 
99  cout << flush;
100 
101  return trigger_luminosities;
102 }

◆ GetLumiFromFiles() [2/2]

std::map<int, TriggerLumi> DAS::GetLumiFromFiles ( const std::filesystem::path &  lumi_file,
const std::filesystem::path &  turnon_file 
)

path to text file with turn-on points

Parameters
lumi_filepath to text file with effective luminosities

◆ getMETfraction()

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

Get MET fraction for each jet trigger separately.

Parameters
inputsinput ROOT file (n-tuple)
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::Options`
steeringparameters obtained from explicit options
slicenumber and index of slice
64  {1,0}
65  )
66 {
67  cout << __func__ << ' ' << slice << " start" << endl;
68 
69  DT::Flow flow(steering, inputs);
70  auto tIn = flow.GetInputTree(slice);
71  auto tOut = flow.GetOutputTree(output);
72 
73  DT::MetaInfo metainfo(tOut);
74  auto isMC = metainfo.Get<bool>("flags", "isMC");
75 
76  auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
77  auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
78  auto met = flow.GetBranchReadOnly<MET>("met");
79  auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
80 
81  vector<double> MET_edges;
82  for (float edge = 0; edge <= 1; edge += 0.01)
83  MET_edges.push_back(edge);
84  int nMETbins = MET_edges.size()-1;
85 
86  auto turnon_file = config.get<fs::path>("corrections.normalisation.turnons");
87  vector<double> trig_edges = GetTriggerTurnons(turnon_file);
88  int nTrigBins = trig_edges.size() - 1;
89 
90  auto METfraction = new TH2F("METfraction", ";MET fraction;p^{leading}_{T} (GeV)",
91  nMETbins, MET_edges.data(),
92  nTrigBins, trig_edges.data());
93 
94  for (DT::Looper looper(tIn); looper(); ++looper) {
95  [[ maybe_unused ]]
96  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
97 
98  if (recJets->empty()) continue;
99  static auto rapAcc = y_edges.back();
100  if (recJets->front().AbsRap() >= rapAcc) continue;
101  if (std::abs(met->SumEt) < feps) continue;
102 
103  auto pt = recJets->front().CorrPt();
104  auto frac = met->Et/met->SumEt;
105  double w = (isMC ? gEv->weights.front().v : 1) * rEv->weights.front().v;
106 
107  cout << setw(15) << pt << setw(15) << frac << setw(15) << w << endl;
108 
109  // we do not use the jet weight,
110  // because the pt is only determined
111  // to find the right trigger window
112  // -> only the event weight matters
113  METfraction->Fill(frac, pt, w);
114  }
115 
116  METfraction->Write();
117 
118  metainfo.Set<bool>("git", "complete", true);
119 
120  cout << __func__ << ' ' << slice << " end" << endl;
121 }

◆ GetObjects()

std::vector<T*> DAS::GetObjects ( TDirectory *  dir)

Find all objects of type T directly in a TDirectory (i.e. not recursive)

88 {
89  if (!dir)
90  BOOST_THROW_EXCEPTION( std::runtime_error("Trying to access inexistent `TDirectory`") );
91  std::vector<T*> objects;
92  for (const auto&& element: *(dir->GetListOfKeys())) {
93  auto key = dynamic_cast<TKey*>(element);
94  if (auto obj = dynamic_cast<T*>(key->ReadObj()); obj)
95  objects.push_back(obj);
96  }
97  return objects;
98 }

◆ GetSmoothFit() [1/2]

TF1* DAS::GetSmoothFit ( TH1 *  h,
double  m,
double  M,
int  nSigmaStop = 0,
const char *  options = "Q0SNRE",
bool  autoStop = false 
)
Parameters
hhistogram to fit (NB: assumed to be normalised to bin width)
mmin
Mmax
nSigmaStopstop automatically before maxdegree if good chi2/ndf +/- nSigma ~ 1
optionsfit options
375 {
376  assert(m < M);
377 
378  int ifm = h->FindBin(m+feps),
379  ifM = h->FindBin(M-feps);
380  cout << "Axis interval (" << m << ',' << M << ") to bin interval (" << ifm << ',' << ifM << ')' << endl;
381 
382  assert(ifm < ifM);
383 
384  return GetSmoothFit<d>(h, ifm, ifM, nSigmaStop, options, autoStop);
385 }

◆ GetSmoothFit() [2/2]

TF1* DAS::GetSmoothFit ( TH1 *  h,
int  ifm,
int  ifM,
int  nSigmaStop = 0,
const char *  options = "Q0SNRE",
bool  autoStop = false 
)
Parameters
hhistogram to fit (NB: assumed to be normalised to bin width)
ifmmin index
ifMmax index
nSigmaStopstop automatically before maxdegree if good chi2/ndf +/- nSigma ~ 1
optionsfit options
273 {
274  assert(ifm > 0);
275  assert(ifM <= h->GetNbinsX());
276  assert(h->GetDimension() == 1);
277 
278  double contentfm = h->GetBinContent(ifm),
279  contentfM = h->GetBinContent(ifM);
280 
281  double m = h->GetBinLowEdge(ifm), M = h->GetBinLowEdge(ifM+1);
282 
283  cout << "From bin " << ifm << " (m = " << m << ", content = " << contentfm << ")"
284  << " to bin " << ifM << " (M = " << M << ", content = " << contentfM << ")" << endl;
285 
286  if (abs(contentfm) <= deps || abs(contentfM) <= deps) {
287  cerr << "Fit cannot start with empty bins\n";
288  //h->Print("all");
289  return nullptr;
290  }
291 
292  double errorfm = h->GetBinError(ifm),
293  errorfM = h->GetBinError(ifM);
294 
295  if (abs(errorfm) <= deps || abs(errorfM) <= deps) {
296  cerr << "There must be uncertainties for the fit\n";
297  //h->Print("all");
298  return nullptr;
299  }
300 
301  double fm = log(contentfm),
302  fM = log(contentfM);
303 
304  double x2n = dmax, error = dmax;
305 
306  static int iname = 0; ++iname;
307  //Greta greta(maxdegree, m, M);
308  Thunberg<d, log, exp> greta(m, M);
309  const size_t maxdegree = greta.maxdegree;
310  TF1 * f = new TF1(Form("greta%d", iname), greta, m, M, greta.npars);
311 
312  // setting initial parameters
313  f->SetParameter(0,(fM+fm)/2);
314  f->SetParameter(1,(fM-fm)/2);
315  for (size_t degree = 2; degree <= maxdegree; ++degree)
316  f->FixParameter(degree,0);
317 
318  auto dump = [f](int degree) {
319  cout << "parameters:\t";
320  for (int i = 0; i <= degree; ++i) cout << f->GetParameter(i) << ' ';
321  cout << endl;
322  };
323 
324  //running
325  for (size_t degree = 2; degree <= maxdegree; ++degree) {
326 
327  dump(degree-1);
328  f->ReleaseParameter(degree);
329 
330  TFitResultPtr r = h->Fit(f, options);
331 
332  double chi2 = r->Chi2();
333  if (!isnormal(chi2)) {
334  h->Print("all");
335  cerr << "Fit aborted (chi2 = " << chi2 << ")\n";
336  continue;
337  }
338  double ndf = r->Ndf();
339 
340  int status = r->Status();
341  if (status != 0) {
342  h->Print("all");
343  cerr << "Fit aborted (status = " << status << ")\n";
344  continue;
345  }
346 
347  auto lastx2n = x2n;
348  x2n = chi2/ndf;
349  error = sqrt(2./ndf);
350  cout << "\e[1m" << degree << ' ' << status << '\t' << chi2 << ' ' << ndf << '\t' << x2n << "\u00B1" << error << "\e[0m" << endl;
351 
352  if (ndf == 2) break;
353  if (autoStop && x2n > lastx2n) {
354  cerr << "Chi2/ndf has increased -- stepping back and stopping!\n";
355  f->FixParameter(degree,0);
356  --degree;
357  h->Fit(f, options);
358  break;
359  }
360  if (abs(x2n-1) <= nSigmaStop*error) break;
361  }
362  dump(maxdegree);
363 
364  return f;
365 }

◆ getSmoothFits()

void DAS::getSmoothFits ( const fs::path &  input,
const fs::path &  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Smooth relative uncertainties on the output of getUnfPhysDist.

Parameters
inputinput ROOT files (histograms)
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
66  {1,0}
67  )
68 {
69  cout << __func__ << " start" << endl;
70 
71  auto fIn = make_unique<TFile>(input.c_str(), "READ");
72 
73  auto fOut = DT::GetOutputFile(output);
74 
75  const TString nDist = config.get<string>("uncertainties.smoothness.distribution"),
76  nCov = config.get<string>("uncertainties.smoothness.covariance");
77  cout << nDist << ' ' << nCov << endl;
78 
79  Step::verbose = steering & DT::verbose;
80  for (const auto& var: GetObjects<TDirectory>(fIn.get())) {
81  // parallelisation (option `-j`)
82  {
83  static int ivar = 0; // used for parallelisation
84  ++ivar;
85  if (slice.second != ivar % slice.first) continue;
86  }
87 
88  TString varname = var->GetName();
89  if (varname.Contains("fake")) continue; // TODO: calculate cov matrix
90  if (varname.Contains("miss")) continue; // TODO: calculate cov matrix
91 
92  // looping over all available observables
93  for (const auto& obs: GetObjects<TDirectory>(var->GetDirectory(nDist))) {
94  TString obsname = obs->GetName();
95  cout << varname << ' ' << obsname << endl;
96 
97  [[ maybe_unused ]] // trick to only activate `cout` in the loop if option `-v` has been given
98  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
99 
100  var->cd();
101  TDirectory * dInDist = GetDirectory(var, {nDist, obsname}),
102  * dInCov = GetDirectory(var, {nCov, obsname, obsname}),
103  * dOutDist = fOut->mkdir(var->GetName(), "", true)
104  ->mkdir(nDist, "", true)
105  ->mkdir(obsname);
106  for (const auto& h: GetObjects<TH1>(dInDist)) {
107  if (h->Integral() <= 0) {
108  delete h;
109  continue;
110  }
111 
112  auto [im, iM] = GetBinRange(h);
113  if (iM - im <= 2) {
114  delete h;
115  continue;
116  }
117 
118  TString n = h->GetName(),
119  covName = n + "__" + n;
120 
121  auto cov = dInCov->Get<TH2>(covName);
122  if (!cov) {
123  h->Print("all");
124  BOOST_THROW_EXCEPTION( DE::BadInput("Covariance matrix "s + covName
125  + " could not be found for "s + varname, *dInCov) );
126  }
127 
128  ScaleWidth(h, cov);
129 
130  using namespace Step;
131  auto f = GetSmoothFit<Basis::Chebyshev,log,exp>(h, cov, im, iM, 10,
132  EarlyStopping::Chi2ndf, 2, 0, cout);
133  // NOTE: degree 10 and 2 sigmas are inspired from personal experience with inclusive jet
134  // --> we may want to make options from the command line in the future
135  dOutDist->cd();
136  for (int i = im; i <= iM; ++i) {
137  double err2 = cov->GetBinContent(i,i);
138  h->SetBinError(i, sqrt(err2));
139  }
140  h->SetDirectory(dOutDist);
141  h->Write();
142  f->Write(Form("fit_%s", h->GetName()));
143 
144  delete h;
145  delete cov;
146  delete f;
147  }
148  }
149  }
150 
151  cout << __func__ << " stop" << endl;
152 }

◆ GetTriggerTurnons()

vector<double> DAS::GetTriggerTurnons ( const fs::path &  f)

Load the two-column file with the HLT and corresponding offline PF thresholds.

41 {
42  if (!fs::exists(f))
43  BOOST_THROW_EXCEPTION(fs::filesystem_error("Bad input", f,
44  make_error_code(errc::no_such_file_or_directory)));
45 
46  pt::ptree thresholds;
47  read_info(f.string(), thresholds);
48 
49  vector<double> edges(1, 30);
50  for (const auto& threshold: thresholds)
51  edges.push_back(threshold.second.get_value<double>());
52  edges.push_back(6500);
53 
54  return edges;
55 }

◆ identity() [1/2]

Eigen::VectorXd DAS::identity ( const Eigen::VectorXd &  x)

identity (Eigen)

Trivial function returning argument, used as default function in the constructor of Teddy

◆ identity() [2/2]

double DAS::identity ( double  x)
inline
23 { return x; };

◆ MakeBinning()

std::vector<T> DAS::MakeBinning ( std::pair< T, T >  range,
double  param,
int  opt = BinningOptions::default_def 
)

Define a regular binning.

Returns
vector of bin edges

We use Boost counting iterator to create a vector out of a sequence

Parameters
rangelowest and highest bin edges
parameither number of bins or bin width (depending on `opt`)
optby default linear binning where `params` is the number of bins
38 {
39  using namespace std;
40 
41  if (range.second <= range.first)
42  BOOST_THROW_EXCEPTION(invalid_argument("Invalid binning range"));
43  if (param <= 0)
44  BOOST_THROW_EXCEPTION(invalid_argument("Invalid binning parameter"));
45  if ((opt & BinningOptions::logarithmic) && range.first <= 0)
46  BOOST_THROW_EXCEPTION(invalid_argument("Logarithmic binning requires range.first > 0"));
47  if ((opt & BinningOptions::logarithmic) && (opt & BinningOptions::step) && param <= 1)
48  BOOST_THROW_EXCEPTION(invalid_argument("Log-step binning requires param > 1"));
49 
50  if (opt & BinningOptions::logarithmic) {
51  range.first = log10(range.first);
52  range.second = log10(range.second);
53  if (opt & BinningOptions::step) param = log10(param);
54  }
55 
56  double w = opt & BinningOptions::step ? param : (range.second-range.first)/param;
58  using boost::counting_iterator;
59  vector<T> edges(counting_iterator<T>(range.first /w),
60  counting_iterator<T>(range.second/w+1));
61  ranges::for_each(edges, [w](T& edge) { edge *= w; });
62  if (opt & BinningOptions::logarithmic)
63  ranges::for_each(edges, [](T& edge) { edge = pow(10,edge); });
64 
65  return edges;
66 }

◆ MakeTitle()

std::vector<TString> DAS::MakeTitle ( const std::vector< double > &  edges,
const char *  v,
bool  lowEdge,
bool  highEdge,
std::function< const char *(double)>  format 
)
inline

Make a vector of properly formated titles from bin edges.

18 {
19  assert(edges.size() > 1);
20  std::vector<TString> titles;
21  for (size_t i = 1; i < edges.size(); ++i) {
22  TString title;
23  if (i > 1 || lowEdge) title += Form("%s < ", format(edges.at(i-1)));
24  title += v;
25  if (i < edges.size()-1 || highEdge) title += Form(" < %s", format(edges.at(i )));
26  titles.push_back(title);
27  }
28  return titles;
29 }

◆ mergeNtuples()

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

Get fractioned n-tuples after the n-tuplisation and merge them into a single file. Complementary information about the pile-up may be included. A minimal selection is applied, e.g. on the primary vertex (PV).

Parameters
inputsinput ROOT file (n-tuple)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
45  {1,0}
46  )
47 {
48  cout << __func__ << ' ' << slice << " start" << endl;
49 
50  DT::Flow flow(steering, inputs);
51  auto tIn = flow.GetInputTree(slice, "ntupliser/events");
52  auto tOut = flow.GetOutputTree(output);
53 
54  DT::MetaInfo metainfo(tOut);
55  metainfo.Check(config);
56  bool isMC = metainfo.Get<bool>("flags", "isMC");
57  int year = metainfo.Get<int>("flags", "year");
58 
59  // event information
60  auto recEvt = flow.GetBranchReadWrite<RecEvent>("recEvent");
61  auto vtx = flow.GetBranchReadWrite<PrimaryVertex>("primaryvertex");
62 
63  // physics object information
64  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets", DT::facultative);
65  auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
66  auto recPhotons = flow.GetBranchReadWrite<vector<RecPhoton>>("recPhotons", DT::facultative);
67 
68  // jet veto business
69  const auto& jetveto_file = config.get<fs::path>("corrections.jetvetomap.filename");
70  optional<JetVeto::Conservative> jetveto;
71  if (jetveto_file != "/dev/null")
72  jetveto = make_optional<JetVeto::Conservative>(jetveto_file);
73  else
74  cout << orange << "No jet veto file was provided." << def << endl;
75  metainfo.Set<fs::path>("corrections", "jetvetomap", "filename", jetveto_file);
76 
77  // MET filter business
78  auto met = flow.GetBranchReadOnly<MET>("met");
79  fs::path p_METfilters = config.get<fs::path>("corrections.METfilters");
80  optional<MissingET::Filters> metfilters;
81  if (p_METfilters != "/dev/null") {
82  if (!metainfo.Find("corrections", "METfilters"))
83  BOOST_THROW_EXCEPTION( DE::BadInput("No MET filters in this tree", tOut) );
84  if (metainfo.Get<fs::path>("corrections", "METfilters") != p_METfilters) {
85  cerr << orange << "You are now using a different list of MET filters w.r.t. n-tuplisation.\n" << def;
86  metainfo.Set<bool>("git", "reproducible", false);
87  }
88  metfilters = make_optional<MissingET::Filters>(p_METfilters);
89  }
90  else
91  cout << orange << "No MET filters will be applied." << def << endl;
92 
93  for (DT::Looper looper(tIn); looper(); ++looper) {
94  [[ maybe_unused ]]
95  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
96 
97  const bool badPV = abs(vtx->z) > 24 /* cm */ || abs(vtx->Rho) > 2 /* cm */ || vtx->fake;
98  if (badPV) {
99  recEvt->weights *= 0;
100  if (recJets) recJets->clear();
101  if (recMuons) recMuons->clear();
102  if (recPhotons) recPhotons->clear();
103  }
104 
105  if (recJets && jetveto)
106  for (auto& recJet: *recJets)
107  (*jetveto)(recJet);
108 
109  if (metfilters) (*metfilters)(*met, *recEvt);
110 
111  if (steering & DT::fill) tOut->Fill();
112  }
113 
114  metainfo.Set<bool>("git", "complete", true);
115 
116  cout << __func__ << ' ' << slice << " end" << endl;
117 }

◆ operator*() [1/7]

double DAS::operator* ( const double  v,
const Weight &  w 
)
inline
32 { return w.v * v; }

◆ operator*() [2/7]

double DAS::operator* ( const float  v,
const Weight &  w 
)
inline
30 { return w.v * v; }

◆ operator*() [3/7]

double DAS::operator* ( const int  v,
const Weight &  w 
)
inline
28 { return w.v * v; }

◆ operator*() [4/7]

double DAS::operator* ( const Weight &  w,
const double  v 
)
inline
31 { return w.v * v; }

◆ operator*() [5/7]

double DAS::operator* ( const Weight &  w,
const float  v 
)
inline
29 { return w.v * v; }

◆ operator*() [6/7]

double DAS::operator* ( const Weight &  w,
const int  v 
)
inline
27 { return w.v * v; }

◆ operator*() [7/7]

double DAS::operator* ( const Weight &  w1,
const Weight &  w2 
)
inline
39 { return w1.v * w2.v; }

◆ operator*=() [1/6]

Weight& DAS::operator*= ( Weight &  w,
const double  v 
)
inline
37 { w.v *= v; return w; }

◆ operator*=() [2/6]

Weight& DAS::operator*= ( Weight &  w,
const float  v 
)
inline
35 { w.v *= v; return w; }

◆ operator*=() [3/6]

Weight& DAS::operator*= ( Weight &  w,
const int  v 
)
inline
33 { w.v *= v; return w; }

◆ operator*=() [4/6]

Weights& DAS::operator*= ( Weights &  wgts,
const double  v 
)
inline
47 { for (auto& w: wgts) w *= v; return wgts; }

◆ operator*=() [5/6]

Weights& DAS::operator*= ( Weights &  wgts,
const float  v 
)
inline
45 { for (auto& w: wgts) w *= v; return wgts; }

◆ operator*=() [6/6]

Weights& DAS::operator*= ( Weights &  wgts,
const int  v 
)
inline
43 { for (auto& w: wgts) w *= v; return wgts; }

◆ operator/=() [1/6]

Weight& DAS::operator/= ( Weight &  w,
const double  v 
)
inline
38 { w.v /= v; return w; }

◆ operator/=() [2/6]

Weight& DAS::operator/= ( Weight &  w,
const float  v 
)
inline
36 { w.v /= v; return w; }

◆ operator/=() [3/6]

Weight& DAS::operator/= ( Weight &  w,
const int  v 
)
inline
34 { w.v /= v; return w; }

◆ operator/=() [4/6]

Weights& DAS::operator/= ( Weights &  wgts,
const double  v 
)
inline
48 { for (auto& w: wgts) w /= v; return wgts; }

◆ operator/=() [5/6]

Weights& DAS::operator/= ( Weights &  wgts,
const float  v 
)
inline
46 { for (auto& w: wgts) w /= v; return wgts; }

◆ operator/=() [6/6]

Weights& DAS::operator/= ( Weights &  wgts,
const int  v 
)
inline
44 { for (auto& w: wgts) w /= v; return wgts; }

◆ operator<()

bool DAS::operator< ( const PhysicsObject &  l,
const PhysicsObject &  r 
)
inline
83 {
84  return l.CorrPt() < r.CorrPt();
85 }

◆ operator==() [1/5]

bool DAS::operator== ( const PhysicsObject &  l,
const PhysicsObject &  r 
)
inline
78 {
79  return l.p4 == r.p4 && l.scales == r.scales && l.weights == r.weights;
80 }

◆ operator==() [2/5]

bool DAS::operator== ( const Weight &  l,
const Weight &  r 
)
inline
26 { return l.v == r.v && l.i == r.i; }

◆ operator==() [3/5]

bool DAS::operator== ( const Weight &  w,
const double  v 
)
inline
25 { return w.v == v; }

◆ operator==() [4/5]

bool DAS::operator== ( const Weight &  w,
const float  v 
)
inline
24 { return w.v == v; }

◆ operator==() [5/5]

bool DAS::operator== ( const Weight &  w,
const int  v 
)
inline
23 { return w.v == v; }

◆ operator>()

bool DAS::operator> ( const PhysicsObject &  l,
const PhysicsObject &  r 
)
inline
88 {
89  return l.CorrPt() > r.CorrPt();
90 }

◆ Options()

DT::Options Options ( const char *  description,
int  options = Darwin::Tools::none 
)

Constructs Darwin options with the correct commit information.

15 {
16 #ifdef DAS_IN_CMSSW
17  // Getting the commit hash is not supported by the scram build system. Pass
18  // a null hash, which means using the Darwin commit hash. This is the best
19  // we can do.
20  return DT::Options(description, options, nullptr);
21 #else
22  return DT::Options(description, options, DAS_COMMIT_HASH);
23 #endif
24 }

◆ pt_sort()

bool DAS::pt_sort ( const PhysicsObject &  j1,
const PhysicsObject &  j2 
)
inline
19 {
20  return j1.CorrPt() > j2.CorrPt();
21 }

◆ ScaleWidth()

void DAS::ScaleWidth ( TH1 *  h,
TH2 *  cov 
)

Divide the histogram by bin width and adapt the covariance matrix accordingly.

46 {
47  for (int i = 1; i <= h->GetNbinsX(); ++i) {
48  double width = h->GetBinWidth(i);
49  double content = h->GetBinContent(i);
50  h->SetBinContent(i, content/width);
51  for (int j = 1; j <= h->GetNbinsX(); ++j) {
52  double width2 = h->GetBinWidth(j);
53  double content = cov->GetBinContent(i,j);
54  cov->SetBinContent(i, j, content/width/width2);
55  }
56  }
57 }

◆ SimpleExec()

void DAS::SimpleExec ( const fs::path &  input,
const fs::path &  output,
const pt::ptree &  config,
const int  steering 
)

Template for function (TODO)

Parameters
inputinput ROOT files (histograms)
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
35 {
36  cout << __func__ << " start" << endl;
37 
38  // TODO
39 
40  cout << __func__ << " stop" << endl;
41 }

Variable Documentation

◆ binsPhi

const std::vector<double> binsPhi = {0.,.1*M_PI, .2*M_PI, .3*M_PI, .4*M_PI, .5*M_PI, .6*M_PI, .7*M_PI, .8*M_PI, .9*M_PI, M_PI}
static

◆ binsY

const std::vector<double> binsY = {0., 0.25, .5, 0.75, 1., 1.25, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 4.25, 4.5, 4.75, 5., 5.25, 5.5, 5.75, 6., 6.25, 6.5, 6.75, 7., 7.25, 7.5, 7.75, 8., 8.25, 8.5, 8.75, 9., 9.25, 9.5}
static

◆ maxMult

const int maxMult = n_edges.size()-1
static

◆ Mjj_edges

const std::vector<double> Mjj_edges {160, 200, 249, 306, 372, 449, 539, 641, 756, 887, 1029, 1187, 1361, 1556, 1769, 2008, 2273, 2572, 2915, 3306, 3754, 4244, 4805, 5374, 6094, 6908, 7861, 8929, 10050 }
static

◆ n_edges

const std::vector<double> n_edges {0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}
static

◆ nMjjBins

const int nMjjBins = Mjj_edges.size()-1
static

◆ nPtBins

const int nPtBins = pt_edges.size()-1
static

◆ nPtDijetBins

const int nPtDijetBins = pt_av_edges.size()-1
static

◆ nPUbins

const int nPUbins = 100
static

◆ nYbins

const int nYbins = y_edges.size()-1
static

◆ pt_av_edges

const std::vector<double> pt_av_edges {147, 175, 207, 243, 284, 329, 380, 437, 499, 569, 646, 732, 827, 931, 1046, 1171, 1307, 1458, 1621, 1806, 2003, 2217, 2453, 2702}
static

◆ pt_edges

const std::vector<double> pt_edges {9, 12, 15,18,21,24,28,32,37,43,49,56,64,74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832,4037}
static

◆ pthat_edges

const std::vector<double> pthat_edges {30,34,39,44,50,56,63,71,80,89,98,108,120,131,143,156,170,196,226,260,300,336,375,420,470,500,531,564,600,645,693,744,800,846,894,946,1000,1088,1183,1287,1400,1491,1587,1690,1800,1934,2078,2233,2400,2579,2771,2978,3200,3578,4000,4472,5000}
static

◆ ptmax

const double ptmax = pt_edges.back()
static

◆ ptmin

const double ptmin = pt_edges.front()
static

◆ SysDown

const std::string SysDown = "Down"
inline

Suffix used for "down" uncertainties. Follows the Combine convention.

◆ SysUp

const std::string SysUp = "Up"
inline

Suffix used for "up" uncertainties. Follows the Combine convention.

◆ y_edges

const std::vector<double> y_edges {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}
static

◆ yBins

const std::vector<TString> yBins = MakeTitle(y_edges, "|y|", false, true, [](double v) { return Form("%.1f", v);} )
static

◆ ymax

const double ymax = y_edges.back()
static

◆ ymin

const double ymin = y_edges.front()
static
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
DAS::Normalisation::phaseSel
std::vector< DAS::RecJet >::iterator phaseSel(std::vector< DAS::RecJet > &recjets)
Definition: PhaseSelection.h:9
method
method
Definition: Core-gitclone-lastrun.txt:5
metPhiCorrectionExample.infile
string infile
Definition: metPhiCorrectionExample.py:25
Ntupliser_cfg.cerr
cerr
Definition: Ntupliser_cfg.py:105
DAS::y_edges
static const std::vector< double > y_edges
Definition: binnings.h:36
Darwin::Tools::fill
@ fill
activate -f to fill the tree
Definition: Options.h:32
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.options
options
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:41
jmarExample.pt
pt
Definition: jmarExample.py:19
dmax
static const double dmax
Definition: Greta.h:19
Darwin::Tools::Flow
User-friendly handling of input and output n-tuples.
Definition: Flow.h:80
Step::verbose
static bool verbose
Definition: Step.h:40
DAS::Normalisation::thresholds
vector< double > thresholds
Definition: getTriggerCurves.cc:41
Step::def
static const char * def
Definition: Step.h:36
Ntupliser_cfg.year
int year
Definition: Ntupliser_cfg.py:67
DAS::logarithmic
@ logarithmic
log with param as number of bins
Definition: toolbox.h:25
DAS::default_def
@ default_def
linear with param as number of bins
Definition: toolbox.h:24
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:22
DAS::step
@ step
linear with param as step for binning
Definition: toolbox.h:26
Ntupliser_cfg.f
f
Definition: Ntupliser_cfg.py:324
jercExample.key
string key
Definition: jercExample.py:109
DAS::ScaleWidth
void ScaleWidth(TH1 *h, TH2 *cov)
Definition: getSmoothTests.cc:45
DAS::MET::Et
float Et
transverse energy
Definition: Event.h:85
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
DAS::GetDirectory
TDirectory * GetDirectory(TDirectory *dir, const std::vector< const char * > &names)
Get ((...)sub)subdirectory of dir.
Definition: toolbox.h:102
deps
static const double deps
Definition: Greta.h:18
DAS::MN::minpt
static const double minpt
Definition: getMNobservables.cc:39
Darwin::Tools::GetOutputFile
std::shared_ptr< TFile > GetOutputFile(const std::filesystem::path &, const std::source_location=std::source_location::current())
Shortcut to create a reproducible output file (see ROOT Doxygen for details)
Definition: FileUtils.cc:147
feps
static const float feps
Definition: Greta.h:17
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:333
DAS::MN::maxpt
static const double maxpt
Definition: getMNobservables.cc:39
orange
static const char * orange
Definition: colours.h:6
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::GetBinRange
pair< int, int > GetBinRange(TH1 *h)
Find range for non-empty bins in an input histogram.
Definition: getSmoothTests.cc:32
Step
Smoothness Test with Expansion of Polynmials.
Definition: Step.h:31
DUMMY
#define DUMMY
Definition: TriggerLumi.cc:9
met
DAS::MET met
Definition: classes.h:28
Step::Chi2ndf
@ Chi2ndf
stop as soon as Chi2ndf has reached a plateau, or as soon as it is compatible with 1
Definition: Step.h:432
Ntupliser_cfg.isMC
dictionary isMC
Definition: Ntupliser_cfg.py:62
feps
static const auto feps
Definition: getMETfraction.cc:34
DAS::GetTriggerTurnons
vector< double > GetTriggerTurnons(const fs::path &f)
Load the two-column file with the HLT and corresponding offline PF thresholds.
Definition: getMETfraction.cc:40
DAS_COMMIT_HASH
#define DAS_COMMIT_HASH
Definition: version.h:1
DAS::Normalisation::threshold
const float threshold
Definition: sigmoid.h:14
jercExample.inputs
def inputs
Definition: jercExample.py:118
DAS::Options
Darwin::Tools::Options Options(const char *, int=Darwin::Tools::none)
Constructs Darwin options with the correct commit information.
Definition: DASOptions.cc:14
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
make_di
auto make_di(Obj1 &o1, Obj2 &o2)
Definition: Di.h:99
metPhiCorrectionExample.range
range
Definition: metPhiCorrectionExample.py:100
DAS::MET::SumEt
float SumEt
sum of the transverse energies of all the components (jet, leptons, etc.) present in the event
Definition: Event.h:86
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83
Darwin::Tools::facultative
@ facultative
mounting branch is facultative
Definition: Flow.h:33
Ntupliser_cfg.dump
string dump
Definition: Ntupliser_cfg.py:46