DAS  3.0
Das Analysis System
Namespaces | Classes | Functions
DAS::Unfolding Namespace Reference

Namespaces

 Checks
 
 DijetMass
 
 DijetMass3D
 
 DrellYan
 
 ExtractHistogram
 
 InclusiveJet
 
 InclusiveJet1D
 
 Migrations
 
 MNjets
 
 Rij
 
 Tikhonov
 
 ZmmY
 

Classes

struct  DistVariation
 
struct  Filler
 
class  MyTUnfoldDensity
 
struct  Observable
 
class  Transformer
 

Functions

template<typename TH1X , typename TH2X >
void CheckSymmetry (const unique_ptr< TH1X > &h, const unique_ptr< TH2X > &cov)
 
template<typename TH1X , typename TH2X >
void CheckEntries (const unique_ptr< TH1X > &h, const unique_ptr< TH2X > &cov)
 
void getToyCalculation (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void getLmatrix (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering)
 
void getUnfBinning (const fs::path &outputGen, const fs::path &outputRec, const pt::ptree &config, const int steering)
 
void getUnfHist (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void getUnfPhysDist (const fs::path &input, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void CorrectInefficiencies (unique_ptr< TH1 > &unf, const map< TString, unique_ptr< TH1 >> &miss)
 
void unfold (const fs::path &inputData, const fs::path &inputMC, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void fillCov (const DistVariation &, const std::list< int > &)
 
std::vector< DistVariation > GetVariations (Darwin::Tools::MetaInfo &, bool=false, std::ostream &=std::cout)
 
template<class T >
auto initOptionalBranch (TTreeReader &reader, const char *name)
 
std::vector< Observable * > GetObservables (boost::property_tree::ptree)
 
template<typename TH2X = TH2>
void SaveCorr (const std::unique_ptr< TH2X > &cov, TString name, TDirectory *d)
 

Function Documentation

◆ CheckEntries()

void DAS::Unfolding::CheckEntries ( const unique_ptr< TH1X > &  h,
const unique_ptr< TH2X > &  cov 
)
59 {
60  cout << __LINE__ << "\tLooking for abnormal entries" << endl;
61  for (int i = 1; i <= h->GetNbinsX(); ++i) {
62  double content = h->GetBinContent(i);
63  if (isfinite(content)) continue;
64  cerr << __LINE__ << orange << "\tWarning: bin " << i
65  << " has abnormal content: " << content << '\n' << def;
66  h->SetBinContent(i, 0);
67  h->SetBinError (i, 0);
68  for (int j = 0; j <= h->GetNbinsX(); ++j) {
69  cov->SetBinContent(i,j,0);
70  cov->SetBinContent(j,i,0);
71  }
72  }
73 }

◆ CheckSymmetry()

void DAS::Unfolding::CheckSymmetry ( const unique_ptr< TH1X > &  h,
const unique_ptr< TH2X > &  cov 
)
41 {
42  cout << __LINE__ << "\tChecking that cov matrix is symmetric" << endl;
43  for (int i = 1; i <= h->GetNbinsX(); ++i)
44  for (int j = 1; j <= h->GetNbinsX(); ++j) {
45  double normalisation = h->GetBinContent(i) * h->GetBinContent(j);
46  if (std::abs(normalisation) < deps) continue;
47  double a_ij = cov->GetBinContent(i,j) / normalisation,
48  a_ji = cov->GetBinContent(j,i) / normalisation;
49 
50  if (a_ji > 0 && std::abs(a_ij / a_ji - 1) > 1e-6)
51  cerr << __LINE__ << orange << "\tWarning: asymmetry found for i = "
52  << i << ", j = " << j << ", getting a_ij = " << a_ij
53  << " and a_ji = " << a_ji << '\n' << def;
54  }
55 }

◆ CorrectInefficiencies()

void DAS::Unfolding::CorrectInefficiencies ( unique_ptr< TH1 > &  unf,
const map< TString, unique_ptr< TH1 >> &  miss 
)

CorrectInefficiencies

Since miss entries are not accounted for in the migrations of the RM, they must be corrected after the unfolding

59 {
60  cout << __LINE__ << '\t' << __func__ << endl;
61  for (int i = 1; i <= unf->GetNbinsX(); ++i) {
62  double content = unf->GetBinContent(i),
63  factor = 1;
64  for (auto& m: miss) factor += m.second->GetBinContent(i);
65  content *= factor;
66  unf->SetBinContent(i, content);
67  // Note: we avoid explicitely the use of TH1::Multiply()
68  // in order to have full control on the propagation
69  // of the uncertainties (done later with cov matrix)
70  }
71  cout << flush;
72 }

◆ fillCov()

void fillCov ( const DistVariation &  v,
const std::list< int > &  binIDs 
)

Fill covariance matrix for a given variation (and its tmp histogram). The list allows to directly access the non-trivial elements of the histogram.

186 {
187  for (auto x: binIDs) for (auto y: binIDs) {
188  double cCov = v.cov->GetBinContent(x,y),
189  cTmp = v.tmp->GetBinContent(x)
190  * v.tmp->GetBinContent(y);
191  v.cov->SetBinContent(x,y,cCov+cTmp);
192  }
193  for (auto x: binIDs) {
194  v.tmp->SetBinContent(x, 0);
195  v.tmp->SetBinError (x, 0);
196  }
197 }

◆ getLmatrix()

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

Make the L-matrix for a given set of observables.

NOTE:

  • the x-axis must correspond to the axis of the TUnfoldBinning binning
  • the y-axis only happens to correspond to it as well in our approach, but primarily corresponds to an a-priori arbitrary number of conditions
Todo:
implement L-matrix and bias variations?
  • the weights will change slightly at gen level
  • the rec level spectrum may slightly change
Todo:
smart pointer?
Todo:
(steering & DT::verbose)
Todo:
the empty bins may change from variation to variation
Parameters
inputinput ROOT file (histograms)
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
45 {
46  cout << __func__ << " start" << endl;
47 
48  // The matrix inversion only corrects for migrations within the phase space,
49  // therefore, we take the projection of the RM on the gen-level axis.
50  // For all systematic variations that only affect the detector level, this
51  // leads the same bias. If one wants to vary the generator level, then one
52  // has to loop over variations as well, since the bias will slightly change
53  // (although this is likely to be a negligible effect).
54  //
58  auto RM = DT::GetHist<TH2>({input}, "nominal/RM");
59  auto bias = unique_ptr<TH1>(RM->ProjectionX("bias"));
60 
61  auto fOut = DT::GetOutputFile(output);
62 
63  cout << "Setting observables" << endl;
64  Observable::isMC = true;
65  auto pt_obs = config.get_child("unfolding.observables");
66  vector<Observable *> observables = GetObservables(pt_obs);
67 
68  auto genBinning = new TUnfoldBinning("gen");
69  for (Observable * obs: observables)
70  genBinning->AddBinning(obs->genBinning);
71 
73 
74  // L-matrix
75  auto L = unique_ptr<TH2>(TUnfoldBinning::CreateHistogramOfMigrations
76  (genBinning, genBinning, "L"));
77  for (Observable * obs: observables) {
78  cout << obs->genBinning->GetName() << endl;
79  obs->setLmatrix(bias, L);
80  }
81 
82  fOut->cd();
83  TDirectory * nom = fOut->mkdir("nominal");
84  nom->cd();
85  L->SetDirectory(nom);
86  L->Write();
87  bias->SetDirectory(nom);
88  bias->Write();
89 
90  cout << __func__ << " stop" << endl;
91 }

◆ GetObservables()

vector< Observable * > GetObservables ( boost::property_tree::ptree  pt)

Get the observables to unfold.

Todo:
smart pointer?
Todo:
not all observables should be fine (e.g. getUnfHist should not accept Rij::Ratios)
Todo:
not all observables should be fine (e.g. getUnfHist should not accept Rij::Ratios)
Parameters
ptexpect labels
44 {
45  vector<Observable *> observables;
46 #define OBS_TYPES (InclusiveJet::PtY)(DijetMass::MjjYmax)(DijetMass3D::MjjYbYs)(Rij::HTn)(Rij::Ratios)(DrellYan::ZPtY)(ZmmY::BF)(ZmmY::DYJetsToMuMu)(ZmmY::DYJetsToTauTau)(ZmmY::TTTo2L2Nu)(ZmmY::QCD)(ZmmY::WW)(ZmmY::WZ)(ZmmY::ZZ)(MNjets::DEtaDPhi)
47 #define IF_OBS(r, data, TYPE) \
48  if (pt.find(BOOST_PP_STRINGIZE(TYPE)) != pt.not_found()) { \
49  observables.push_back(new TYPE); \
50  for (auto it = pt.begin(); it != pt.end(); ++it) { \
51  if (it->first != BOOST_PP_STRINGIZE(TYPE)) continue; \
52  pt.erase(it); \
53  break; \
54  } \
55  }
56  BOOST_PP_SEQ_FOR_EACH(IF_OBS, _, OBS_TYPES)
57 #undef OBS_TYPES
58 #undef IF_OBS
59  if (!pt.empty()) {
60  TString unknown;
61  for (auto it = pt.begin(); it != pt.end(); ++it) {
62  unknown += it->first;
63  unknown += ' ';
64  }
65  BOOST_THROW_EXCEPTION(invalid_argument(Form("Unknown observable(s): %s", unknown.Data())));
66  }
67  return observables;
68 }

◆ getToyCalculation()

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

Use a toy calculation to apply some operations to the input observables.

The output file keeps (as much as possible) the same structure as the input.

Todo:
once the meta info is available, use metainfo.Seed<67325879>(slice)
Parameters
inputinput ROOT file
outputoutput ROOT file
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
85  {1,0}
86  )
87 {
88  Teddy::verbose = (steering & DT::verbose) == DT::verbose;
89 
90  auto fIn = make_unique<TFile>(input.c_str(), "READ");
91 
92  auto fOut = DT::GetOutputFile(output);
93 
94  auto pt_obs = config.get_child("unfolding.observables");
95  cout << "Setting observables" << endl;
96  vector<Observable *> observables = GetObservables(pt_obs);
97 
98  const TString hName = config.get<string>("uncertainties.toy.distribution"),
99  covName = config.get<string>("uncertainties.toy.covariance"),
100  level = config.get<string>("uncertainties.toy.level");
101 
102  const bool genLevel = level.Contains("gen");
103  if (!genLevel && !level.Contains("rec"))
104  BOOST_THROW_EXCEPTION( invalid_argument("Level was unrecognized") );
105 
106  TUnfoldBinning * pre = new TUnfoldBinning("pre");
107  TUnfoldBinning * post = new TUnfoldBinning("post");
108  vector<unique_ptr<Transformer>> transformers;
109  for (const Observable * obs: observables) {
110  // pre
111  auto preBinning = genLevel ? obs->genBinning : obs->recBinning;
112  pre->AddBinning(preBinning);
113  // post
114  unique_ptr<Transformer> t = obs->getTransformer(preBinning);
115  TUnfoldBinning * postSubBinning = t->postBinning;
116  post->AddBinning(postSubBinning);
117  transformers.push_back(move(t));
118  }
119 
120  auto N = config.get<int>("uncertainties.toy.N");
121 
122  fIn->cd();
123  // loop over the entries of a TFile / TDirectory
124  for (const auto&& obj: *(fIn->GetListOfKeys())) {
125 
126  auto key = dynamic_cast<TKey*>(obj);
127  const TString classname = key->ReadObj()->ClassName();
128  cout << __LINE__ << '\t' << key->GetName() << endl;
129  if (!classname.Contains("TDirectory")) continue;
130  auto dirIn = dynamic_cast<TDirectory*>(key->ReadObj());
131  dirIn->cd();
132 
133  auto hIn = unique_ptr<TH1D>(dirIn->Get<TH1D>( hName));
134  auto covIn = unique_ptr<TH2D>(dirIn->Get<TH2D>(covName));
135  if (!hIn)
136  BOOST_THROW_EXCEPTION( DE::BadInput("Can't find" + hName, *dirIn) );
137  if (!covIn)
138  BOOST_THROW_EXCEPTION( DE::BadInput("Can't find" + covName, *dirIn) );
139 
140  cout << __LINE__ << "\t`hIn` has " << hIn->GetNbinsX() << " bins\n"
141  << __LINE__ << "\t`covIn` has " << covIn->GetNbinsX()
142  << 'x' << covIn->GetNbinsY() << " bins" << endl;
143 
144  {
145  static int ivar = 0; // used for parallelisation
146  ++ivar;
147  if (slice.second != ivar % slice.first) continue;
148  }
149 
150  cout << __LINE__ << "\tRemove bad input bins" << endl;
151  for (auto& t: transformers)
152  t->RemoveBadInputBins(hIn.get(), covIn.get());
153 
154  cout << __LINE__ << "\tChecking input" << endl;
155  CheckSymmetry(hIn, covIn);
156  CheckEntries(hIn, covIn);
157 
158  int nBinsOut = 0;
159  for (auto node = post->GetChildNode(); node != nullptr;
160  node = node->GetNextNode()) {
161  int nBinsHere = node->GetDistributionNumberOfBins();
162  cout << __LINE__ << '\t' << node->GetName() << " has " << nBinsHere << " bins\n";
163  nBinsOut += nBinsHere;
164  }
165  cout << __LINE__ << "\tnBinsOut = " << nBinsOut << endl;
166  function<VectorXd(const VectorXd &)> transformation =
167  [nBinsOut,&transformers](const VectorXd& x) {
168  Transformer::y = VectorXd::Zero(nBinsOut);
169  for (const auto& t: transformers)
170  t->Transform(x);
171  return Transformer::y;
172  };
173 
174  cout << __LINE__ << "\tDeclaring toy" << endl;
175  auto seed = 42;
176  Teddy toy(hIn, covIn, transformation, seed);
177 
178  cout << __LINE__ << "\tPlaying" << endl;
179  for (long i = 0; i < N; ++i)
180  toy.play();
181 
182  cout << __LINE__ << "\tGetting output" << endl;
183  auto [hOut, covOut] = toy.buy("trans", nBinsOut, 0.5, nBinsOut + 0.5);
184  hOut->SetTitle(hIn->GetTitle());
185  covOut->SetTitle(covIn->GetTitle());
186 
187  cout << __LINE__ << "\tChecking output" << endl;
188  CheckSymmetry(hOut, covOut);
189  CheckEntries(hOut, covOut);
190 
191  fOut->cd();
192 
193  cout << __LINE__ << "\tWriting output" << endl;
194  auto dirOut = fOut->mkdir(dirIn->GetName());
195  dirOut->cd();
196  hOut->SetDirectory(dirOut);
197  hOut->Write(hName);
198  covOut->SetDirectory(dirOut);
199  covOut->Write(covName);
200  TString corrName = covName;
201  corrName.ReplaceAll("cov", "corr");
202  corrName.ReplaceAll("trans_", "");
203  SaveCorr<TH2D>(covOut, corrName, dirOut);
204  }
205 
206  cout << __LINE__ << "\tDone" << endl;
207 }

◆ getUnfBinning()

void DAS::Unfolding::getUnfBinning ( const fs::path &  outputGen,
const fs::path &  outputRec,
const pt::ptree &  config,
const int  steering 
)

Get binning in XML format.

Todo:
get post-transformation unfolding
Parameters
outputGenoutput XML file
outputRecoutput XML file
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
38 {
39  cout << __func__ << " start" << endl;
40 
41  cout << "Setting observables" << endl;
42  auto pt_obs = config.get_child("unfolding.observables");
43  vector<Observable *> observables = GetObservables(pt_obs);
44 
46 
47  cout << "Adding binnings" << endl;
48  TUnfoldBinningXML * genBinning = new TUnfoldBinningXML("gen"),
49  * recBinning = new TUnfoldBinningXML("rec");
50  for (Observable * obs: observables) {
51  genBinning->AddBinning(obs->genBinning);
52  recBinning->AddBinning(obs->recBinning);
53  }
54 
55  cout << "Exporting" << endl;
56  if (outputGen != "/dev/null")
57  genBinning->ExportXML(outputGen.c_str());
58  if (outputRec != "/dev/null")
59  genBinning->ExportXML(outputRec.c_str());
60 
61  cout << __func__ << " stop" << endl;
62 }

◆ getUnfHist()

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

Get histograms for unfolding only for pt (in d=1), without using TUnfoldBinning.

Then the output can be unfolded using unfold.

Todo:
First-class support in Darwin
Parameters
inputsinput ROOT files (n-tuple)
outputoutput ROOT file
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);
45  auto tIn = flow.GetInputTree(inputs, slice);
46  auto [fOut, tOut] = flow.GetOutput(output);
47 
48  DT::MetaInfo metainfo(tOut);
49  bool isMC = metainfo.Get<bool>("flags", "isMC");
50  int R = metainfo.Get<int>("flags", "R");
51 
52  TTreeReader reader(tIn);
53 
54  cout << "Setting observables" << endl;
55  Observable::isMC = isMC;
56  Observable::maxDR = R / 20.;
57  auto pt_obs = config.get_child("unfolding.observables");
58  vector<Observable *> observables = GetObservables(pt_obs);
59  vector<unique_ptr<Filler>> fillers;
60  for (const auto obs: observables)
61  fillers.push_back(obs->getFiller(reader));
62 
63  cout << "Setting variations" << endl;
64  DistVariation::isMC = isMC;
65  DistVariation::genBinning = new TUnfoldBinning("gen");
66  DistVariation::recBinning = new TUnfoldBinning("rec");
67  for (Observable * obs: observables) {
68  DistVariation::genBinning->AddBinning(obs->genBinning);
69  DistVariation::recBinning->AddBinning(obs->recBinning);
70  }
71 
72  [[ maybe_unused ]]
73  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
74 
75  vector<DistVariation> variations = GetVariations(metainfo,
76  steering & DT::syst, cout);
77 
78  // See Darwin::Tools::Looper constructor
80  const auto entries = reader.GetEntries();
81  const long long start = entries * ((slice.second + 0.) / slice.first);
82  const long long stop = entries * ((slice.second + 1.) / slice.first);
83  reader.SetEntriesRange(start, stop);
84  while (reader.Next()) {
85  if (auto status = reader.GetEntryStatus(); status != TTreeReader::kEntryValid)
86  BOOST_THROW_EXCEPTION(
87  DE::BadInput(Form("Loading entry %lld failed: error %d",
88  reader.GetCurrentEntry(), status),
89  tIn) );
90 
91  for (DistVariation& v: variations) {
92 
93  list<int> binIDs; // collect filled bins to calculate covariance
94  for (auto& filler: fillers)
95  binIDs.merge(filler->fillRec(v));
96  fillCov(v, binIDs);
97 
98  if (!isMC) continue;
99 
100  // run matching (possibly different for each observable)
101  for (auto& filler: fillers)
102  filler->match();
103 
104  // fill gen, RM, fake, miss
105  for (auto& filler: fillers)
106  filler->fillMC(v);
107  }
108  }
109 
110  fOut->cd();
111  for (DistVariation& v: variations)
112  v.Write(fOut);
113 
114  cout << __func__ << ' ' << slice << " stop" << endl;
115 }

◆ getUnfPhysDist()

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

Extract distributions with physical binning from output of the various getUnf* commands.

Parameters
inputinput ROOT file (histograms)
outputoutput ROOT file (histograms)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
240  {1,0}
241  )
242 {
243  cout << __func__ << " start" << endl;
244 
245  cout << "Setting observables" << endl;
246  auto pt_obs = config.get_child("unfolding.observables");
247  vector<Observable *> observables = GetObservables(pt_obs);
248 
249  cout << "Setting up axes" << endl;
250  auto genBinning = new TUnfoldBinning("gen"),
251  recBinning = new TUnfoldBinning("rec");
252  for (Observable * obs: observables) {
253  genBinning->AddBinning(obs->genBinning);
254  recBinning->AddBinning(obs->recBinning);
255  }
256 
257  // use the axis name to guess which binning is used
258  auto getBinning = [genBinning,recBinning](TString name) {
259  if (name == "gen") return genBinning;
260  if (name == "rec") return recBinning;
261  BOOST_THROW_EXCEPTION( runtime_error("The binning could "
262  "not be determined from \"" + name + "\"") );
263  };
264 
265  unique_ptr<TFile> fIn(TFile::Open(input.c_str(), "READ"));
266  auto fOut = DT::GetOutputFile(output);
267 
268  cout << "Looping over input `TDirectory`s (which in principle correspond to variations)" << endl;
269  for (const auto& dIn: GetObjects<TDirectory>(fIn.get())) {
270  [[ maybe_unused ]] // trick to only activate `cout` in the loop if option `-v` has been given
271  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
272 
273  // parallelisation (option `-j`)
274  {
275  static int ivar = 0; // used for parallelisation
276  ++ivar;
277  if (slice.second != ivar % slice.first) continue;
278  }
279 
280  TString dirname = dIn->GetName();
281  if (!(steering & DT::syst) && dirname != "nominal") continue;
282  auto dOut = fOut->mkdir(dIn->GetName());
283  dIn->cd();
284  dOut->SetTitle(dIn->GetTitle());
285 
286  cout << "Looping over input `TH1`s" << endl;
287  for (TH1 * hist: GetObjects<TH1>(dIn)) {
288  if (TString(hist->ClassName()).Contains("TH2")) continue;
289 
290  TString name = hist->GetName(),
291  axis = hist->GetXaxis()->GetTitle();
292  name.ReplaceAll(dirname, ""); // for some weird reason, ROOT remembers the very first
293  // name of the object, and not the one used in the TFile
294  cout << name << ' ' << axis << endl;
295 
296  TUnfoldBinning * binning = getBinning(axis);
297 
298  auto ddOut = dOut->mkdir(name, hist->GetTitle());
299  ExtractHistogram::Dist dist{hist};
300  ExtractHistogram::nodeLoop(ddOut, binning, dist, cout);
301 
302  delete hist;
303  }
304 
305  cout << "Looping over input `TH2`s" << endl;
306  for (TH2 * hist: GetObjects<TH2>(dIn)) {
307  TString name = hist->GetName(),
308  xaxis = hist->GetXaxis()->GetTitle(),
309  yaxis = hist->GetYaxis()->GetTitle();
310  name.ReplaceAll(dirname, ""); // for some weird reason, ROOT remembers the very first
311  // name of the object, and not the one used in the TFile
312  cout << name << ' ' << xaxis << ' ' << yaxis << endl;
313 
314  TUnfoldBinning * binning1 = getBinning(xaxis),
315  * binning2 = getBinning(yaxis);
316 
317  auto ddOut = dOut->mkdir(name, dIn->GetTitle());
318  ExtractHistogram::Matrix mat(hist, binning2);
319  ExtractHistogram::nodeLoop(ddOut, binning1, mat, cout);
320 
321  delete hist;
322  }
323  }
324 
325  cout << __func__ << " stop" << endl;
326 }

◆ GetVariations()

vector< DistVariation > GetVariations ( Darwin::Tools::MetaInfo &  metainfo,
bool  applySyst = false,
std::ostream &  cout = std::cout 
)

Get all variations availables according to metainfo.

Parameters
metainfofull meta info including variation block
applySystapply systematics
201 {
202  using namespace DAS::Unfolding;
203 
204  vector<DistVariation> variations;
205  variations.reserve(100);
206  DistVariation nominal{"nominal","nominal"};
207  variations.push_back(move(nominal)); // TODO: emplace_back?
208  if (!applySyst) return variations;
209 
210  const TList * groupList = metainfo.List("variations");
211  for (const TObject * group: *groupList) {
212  const auto groupContents = dynamic_cast<const TList*>(group);
213  if (!groupContents)
214  BOOST_THROW_EXCEPTION( invalid_argument(group->GetName()) );
215  size_t i = 0;
216  for (const TObject * obj: *groupContents) {
217  const auto string = dynamic_cast<const TObjString*>(obj);
218  if (!string)
219  BOOST_THROW_EXCEPTION( invalid_argument(obj->GetName()) );
220  cout << group->GetName() << " " << string->GetString() << " " << (i+1) << endl;
221  DistVariation v(group->GetName(), string->GetString(), ++i);
222  variations.push_back(move(v)); // TODO: emplace_back?
223  }
224  }
225 
226  return variations;
227 }

◆ initOptionalBranch()

auto DAS::Unfolding::initOptionalBranch ( TTreeReader &  reader,
const char *  name 
)

Utility function to init an optional TTreeReaderArray/TTreeReaderValue. Example:

struct A
{
std::optional<TTreeReaderValue<int>> value;
A(TTreeReader& reader)
: value(initOptionalBranch<decltype(value)>(reader, "branchName"))
{}
};
Copyright
Originally written for the Shears framework.
40 {
41  return reader.GetTree()->GetBranch(name) != nullptr
42  ? std::make_optional<typename T::value_type>(reader, name)
43  : std::nullopt;
44 }

◆ SaveCorr()

void DAS::Unfolding::SaveCorr ( const std::unique_ptr< TH2X > &  cov,
TString  name,
TDirectory *  d 
)

Save correlation matrix.

15 {
16  using namespace std;
17  auto corr = unique_ptr<TH2>(dynamic_cast<TH2*>(cov->Clone(name)));
18  corr->Reset();
19 
20  TString title = cov->GetTitle();
21  title.ReplaceAll("covariance", "correlation");
22  corr->SetTitle(title);
23 
24  const int Nx = corr->GetNbinsX(),
25  Ny = corr->GetNbinsY();
26 
27  for (int x = 1; x <= Nx; ++x) {
28  double denx = sqrt(cov->GetBinContent(x, x));
29  for (int y = 1; y <= Ny; ++y) {
30  double deny = sqrt(cov->GetBinContent(y, y));
31 
32  if (denx > 0 && deny > 0) {
33  double content = cov->GetBinContent(x, y);
34  content /= denx*deny;
35  corr->SetBinContent(x, y, content);
36  if (x != y && abs(content) >= 1.) cout << x << ' ' << y << ' ' << content << endl;
37  }
38  else
39  corr->SetBinContent(x, y, 0);
40  }
41  }
42  corr->SetMinimum(-1);
43  corr->SetMaximum( 1);
44  corr->SetDirectory(d);
45  d->cd();
46  corr->Write(name);
47 }

◆ unfold()

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

Unfolding output from one of the getUnfHistograms*.C executables (indifferently)

Note that this code does not require TUnfoldBinning. Indeed, TUnfoldBinning is used to make the mapping from ND to 1D in other executables, but it is not necessary to make the unfolding itself. In fact, on purpose, we do not want to know the meaning of the distribution in this executable, to keep it as generic as possible.

Comment on the structure of the code:

  • nested loops, one over the variations in MC, the other over the variations in data
  • we vary either MC or Data, but never both at the same time

  • we keep the same structure of directories as in the inputs, and add additional directories corresponding to unfolding variations, and additional control distributions about the unfolding

    Todo:
    DAgostini (at least for nominal value)
Todo:
merge the MetaInfo from the two (or more??) inputs
  • compare flags, git, corrections
    • merge
    • send warning if relevant but proceed anyway
  • merge histories & variations
  • find a solution for the preseed (Q: is it necessary for the data?)
  • then add block on unfolding but preserve the cycle
Todo:
if (!RM) BOOST_THROW_EXCEPTION( DE::BadInput("Missing response matrix", dirMC) );
Todo:
Propagate metainfo
Todo:
Eps matrix of MyTUnfoldDensity != 1e-100?
Todo:
syst?
Todo:
run CT if the gen level is available (= pseudodata)
  • if inputData = inputMC -> expect perfect agreement
  • otherwise, interpretation may depend...
Todo:
remove?
Todo:
title of new histogram
Todo:
cov? (needed for smoothing)
Todo:
title of new histogram
Todo:
cov? (needed for smoothing)
Todo:
cov? (needed for smoothing)
Todo:
cov? (needed for smoothing)
Parameters
inputDatainput ROOT files (histograms)
inputMCinput ROOT files (histograms)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::Options`
steeringparameters obtained from explicit options
slicenumber and index of slice
442  {1,0}
443  )
444 {
445  cout << __func__ << ' ' << slice << " start" << endl;
446 
447  auto sysUncorr = config.get_optional<fs::path>("unfolding.sysUncorr"),
448  Lmatrix = config.get_optional<fs::path>("unfolding.TUnfold.regularisation.Lmatrix");
449  const bool applySyst = steering & DT::syst;
450 
451  cout << __LINE__ << "\tOpening input files" << endl;
452  auto fData = make_unique<TFile>(inputData.c_str(), "READ"),
453  fMC = make_unique<TFile>(inputMC .c_str(), "READ");
454 
462 
463  cout << __LINE__ << "\tPreparing output file" << endl;
464  auto fOut = DT::GetOutputFile(output);
465 
466  fMC->cd();
467  // loop over the entries of a TFile / TDirectory
468  for (const auto&& objMC: *(fMC->GetListOfKeys())) {
469 
470  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
471 
472  auto keyMC = dynamic_cast<TKey*>(objMC);
473  const TString classname = keyMC->ReadObj()->ClassName();
474  //if (classname == "TTree") {
475  // unique_ptr<TTree> tree = ReadObj<TTree>(keyMC);
476  // auto List = tree->GetUserInfo();
477  // fOut->cd();
478  // List->Write("MC");
479  // fMC->cd();
480  //}
481  if (!classname.Contains("TDirectory")) continue;
482 
483  const TString MCvar = objMC->GetName();
484  cout << __LINE__ << '\t' << MCvar << endl;
485  if (MCvar != "nominal" && !applySyst) continue;
486 
487  auto dirMC = dynamic_cast<TDirectory*>(keyMC->ReadObj());
488  dirMC->cd();
489 
490  cout << __LINE__ << "\tExtracting miss and fake counts" << endl;
491  map<TString, unique_ptr<TH1>> miss, // only to be added by hand at the end
492  fake; // seen as backgrounds from the point of view of TUnfold
493  // reminder: in principle, one source of miss / fake is enough for the unfolding,
494  // but it's good to be able to disentangle the different sources
495 
496  // loop over the entries of a TFile / TDirectory
497  // here, we want to retrieve only miss and fake
498  for (const auto&& objMC2: *(dirMC->GetListOfKeys())) {
499  auto keyMC = dynamic_cast<TKey*>(objMC2);
500  if (!TString(keyMC->ReadObj()->ClassName()).Contains("TH") ) continue;
501 
502  unique_ptr<TH1> h = ReadObj<TH1>(keyMC);
503  TString name = h->GetName();
504  if (!name.Contains("miss") && !name.Contains("fake")) continue;
505  if (h->GetEntries() == 0)
506  cerr << orange << "`" << name << "` is empty\n" << def;
507  name.ReplaceAll(MCvar, "");
508  h->SetDirectory(nullptr);
509  if (name.Contains("fake")) {
510  name.ReplaceAll("fake", "");
511  cout << __LINE__ << "\tfake\t" << name << endl;
512  fake[name] = move(h);
513  }
514  else {
515  name.ReplaceAll("miss", "");
516  cout << __LINE__ << "\tmiss\t" << name << endl;
517  miss[name] = move(h);
518  }
519  } // end of retrieving of miss/fake histograms from MC
520 
521  cout << __LINE__ << "\tExtracting RM from MC" << endl;
522  auto RM = unique_ptr<TH2>(dirMC->Get<TH2>("RM"));
524  RM->SetDirectory(nullptr);
525 
526  auto genMC = unique_ptr<TH1>(dirMC->Get<TH1>("gen")),
527  recMC = unique_ptr<TH1>(dirMC->Get<TH1>("rec"));
528  genMC->SetName("genMC");
529  recMC->SetName("recMC");
530  genMC->SetDirectory(nullptr);
531  recMC->SetDirectory(nullptr);
532 
533  unique_ptr<TH2> PM = Checks::GetPM(RM);
534  PM->SetDirectory(nullptr);
535  Checks::folding(genMC, recMC, PM, miss, fake);
536 
537  cout << __LINE__ << "\tCalculating miss and fake rates" << endl;
538 
539  // we want a fake/miss >> RATE << (i.e. a value between 0 and 1)
540  // -> we will use it later to estimate the fake/miss contribution
541  // in the data using a bin-by-bin method.
542  // currently: the elements of `fake`/`miss` are currently fake/miss *counts*
543  // (including model dependence)
544  unique_ptr<TH1> genMC_noMiss = Clone(genMC, "genMC_noMiss");
545  for (auto& f: fake) f.second ->Divide(f.second.get(), recMC.get(), 1, 1, "b");
546  for (auto& m: miss) genMC_noMiss->Add (m.second.get(), -1);
547  for (auto& m: miss) m.second ->Divide(m.second.get(), genMC_noMiss.get(), 1, 1, "b");
548  // now: the elements of `fake`/`miss` are currently fake/miss *rates*
549 
550  // now we are done with retrieving objects from fMC, and we move to fData
551 
552  fData->cd();
553 
554  // now we estimate the fake contribution in data with bin-by-bin method
555  {
556  auto rec = unique_ptr<TH1>(fData->Get<TH1>("nominal/rec"));
557  // reminder: the elements of `fake` are currently fake *rates*
558  for (auto& f: fake) f.second->Multiply(rec.get());
559  // now: the elements of `fake` are currently fake *counts* (corrected to data)
560  }
561  // (similar thing will be done with the miss contribution *after* the unfolding)
562 
563  // reminder: now we have absolute fake count BUT miss rate
564 
565  cout << __LINE__ << "\tExtracting data variations and running unfolding" << endl;
566 
567  // here we loop again on the entries of a file
568  for (const auto&& objData: *(fData->GetListOfKeys())) {
569 
570  auto keyData = dynamic_cast<TKey*>(objData);
571  const TString classname = keyData->ReadObj()->ClassName();
573  //if (classname == "TTree") {
574  // unique_ptr<TTree> tree = ReadObj<TTree>(keyData);
575  // auto List = tree->GetUserInfo();
576  // fOut->cd();
577  // List->Write("Data");
578  // fMC->cd();
579  //}
580  if (!classname.Contains("TDirectory")) continue;
581 
582  const TString DATAvar = objData->GetName();
583  if (DATAvar != "nominal") {
584  if (!applySyst) continue;
585  if (MCvar != "nominal") continue;
586  }
587  cout << __LINE__ << "\tvariation: " << DATAvar << endl;
588  TString var = MCvar == "nominal" ? DATAvar : MCvar;
589 
590  if (var != "nominal" && MCvar == DATAvar)
591  cerr << orange << "The same variation has been found "
592  "in both data and MC.\n" << def;
593 
594  {
595  static int ivar = 0; // used for parallelisation
596  ++ivar;
597  if (slice.second != ivar % slice.first) continue;
598  }
599 
600  auto dirData = dynamic_cast<TDirectory*>(keyData->ReadObj());
601  dirData->cd();
602 
603  auto recData = unique_ptr<TH1>(dirData->Get<TH1>("rec"));
604  recData->SetName("recData");
605  recData->SetTitle("detector level (" + var + ")");
606  recData->SetDirectory(nullptr);
607  auto covInData = unique_ptr<TH2>(dirData->Get<TH2>("cov"));
608  covInData->SetName("covInData");
609  covInData->SetTitle("covariance matrix at detector level (" + var + ")");
610  covInData->SetDirectory(nullptr);
611  auto covInMC = unique_ptr<TH2>(dirMC->Get<TH2>("cov"));
612  covInMC->SetName("covInMC");
613  covInMC->SetTitle("covariance matrix at detector level (" + var + ")");
614  covInMC->SetDirectory(nullptr);
615 
616  cout << __LINE__ << '\t' << "\e[1m" << DATAvar << "\e[0m" << endl;
617 
618  if (DATAvar.Contains("nominal")) { // if data variation is nominal,
619  // but the condition number will
620  // be calculated for each variation of MC
621  auto condition = Migrations::Condition(RM);
622  RM->SetTitle(Form("%.1f",condition)); // keep track of the condition number
623  }
624 
625  // underflow and overflow are expected to be empty (else, set to 0 with a warning)
626  Checks::CovMargin(recData, covInData, "underflow", 0);
627  Checks::CovMargin(recData, covInData, "overflow", recData->GetNbinsX()+1);
628 
629  // so far, still no unfolding done
630 
631  cout << __LINE__ << "\tSetting up density & input (checking "
632  "and comparing the bin contents of data and MC)" << endl;
633 
634  auto density = make_unique<MyTUnfoldDensity>(
635  RM.get(), // reminder: this is varied by the first loop
636  TUnfold::kHistMapOutputHoriz, // truth level on x-axis of the response matrix
637  TUnfold::kRegModeNone, // do not use *predefined* regularisation schemes
638  TUnfold::kEConstraintNone, // no constraint on area
639  TUnfoldDensity::kDensityModeBinWidthAndUser); // no scale factors
640 
641  // reminder: the epsilon is the smallest possible number
642  density->SetEpsMatrix(1e-100);
643 
644  int status = density->SetInput(
645  recData.get(), // reconstructed level
646  0, // bias factor (only relevant for Tikhonov regularisation)
647  0, // for bins with zero error, this number defines 1/error
648  covInData.get()); // covariance matrix
649 
650  if (status > 0)
651  cerr << __LINE__ << '\t' << orange << "unfold status = " << status << '\n' << def;
652 
653  // now, we should have (nearly) everything settled for the unfolding...
654 
655  cout << __LINE__ << "\tSubtract background" << endl;
656  auto fakeNormVar = config.get_optional<float>("unfolding.fakeNormVar").value_or(0.05);
657  cout << __LINE__ << "\tfakeNormVar = " << fakeNormVar << endl;
658  for (auto& f: fake) { // fake now contains esimations of the fake contributions in data
659  cout << __LINE__ << '\t' << f.first << endl;
660  density->SubtractBackground(
661  f.second.get(), // 1D histogram in rec-level binning (with its own stat use
662  // --> we will use later `GetEmatrixSysBackgroundUncorr()`
663  "fake" + f.first, // name
664  1.0, // scale factor
665  fakeNormVar); // factor to vary the normalisation
666  // --> we will use later `GetDeltaSysBackgroundScale()`
667  // the lines just above are equivalent to three things:
668  // - `h->Add(f.second, -1);`, i.e. subtract the fake contribution
669  // from the rec level
670  // - it propagates the contribution from the limited statistics
671  // to the covariance matrix used in the unfolding
672  // - TUnfold provides systematic variations
673  }
674 
675  if (sysUncorr) {
676  cout << __LINE__ << "\tAdding bin-to-bin uncertainty" << endl;
677 
678  if (!fs::exists(*sysUncorr))
679  BOOST_THROW_EXCEPTION(fs::filesystem_error("File with systematic "
680  "bin-to-bin uncorrelated uncertainties was not found!",
681  make_error_code(errc::no_such_file_or_directory)));
682 
683  auto fSys = make_unique<TFile>(sysUncorr->c_str(), "READ");
684  auto hSysUncorr = unique_ptr<TH1>( fSys->GetDirectory(DATAvar)
685  ->Get<TH1>("sysUncorr1D") );
686  density->SubtractBackground(hSysUncorr.get(), "sysUncorr", 1.0, 0.);
687  }
688 
689  fOut->cd();
690  if (fOut->GetDirectory(var)) continue;
691  auto dOut = fOut->mkdir(var);
692  dOut->cd();
693 
694  cout << __LINE__ << "\tUnfolding" << endl;
695  static auto tau = config.get_optional<float>("unfolding.TUnfold.regularisation.tau");
696  if (Lmatrix) {
697  cout << __LINE__ << "\tFetching L-matrix" << endl;
698 
699  if (!fs::exists(*Lmatrix))
700  BOOST_THROW_EXCEPTION(fs::filesystem_error("File with L-matrix not found!",
701  make_error_code(errc::no_such_file_or_directory)));
702 
703  Tikhonov::PrepareRegularisation("nominal", Lmatrix->c_str(), density);
704  if (!tau && var == "nominal") {
705 
706  tau = Tikhonov::ScanLcurve(density, dOut, config);
707 
708  if (applySyst) {
709  cout << __LINE__ << "\tSaving shifts for variations of tau parameter" << endl;
710  Tikhonov::RedoUnfold(density, *tau*2, "Regularisation" + SysUp,
711  "upper variation of regularisation", fOut, miss);
712  Tikhonov::RedoUnfold(density, *tau/2, "Regularisation" + SysDown,
713  "lower variation of regularisation", fOut, miss);
714  }
715  }
716  } // end of regularise
717  cout << __LINE__ << "\ttau = " << tau << endl;
718  density->DoUnfold(tau.value_or(0));
719 
720  // now that the unfolding is done, retrieve the result (use TUnfold)
721  // and apply the correction for miss entries (by hand)
722 
723  cout << __LINE__ << "\tGetting output" << endl;
724  auto unf = unique_ptr<TH1>(density->GetOutput("unfold", "particle level (" + var + ")"));
725  Checks::NegativeBins(unf);
726  CorrectInefficiencies(unf, miss);
727  dOut->cd();
728  unf->SetDirectory(dOut);
729  unf->GetXaxis()->SetTitle("gen");
730  unf->Write();
731 
735 
736  cout << __LINE__ << "\tCovariance matrix from measurement" << endl;
737  auto covOutData = unique_ptr<TH2>(density->GetEmatrixInput("covOutData", // name
738  "covariance matrix coming from input data distribution (" + var + ")")); // title
739  covOutData->SetDirectory(dOut);
740  covOutData->GetXaxis()->SetTitle("gen");
741  covOutData->GetYaxis()->SetTitle("gen");
742  covOutData->Write();
743 
744  cout << __LINE__ << "\tSaving correlation matrix from measurement" << endl;
745  SaveCorr(covInData , "corrInData" , dOut); // transform cov matrix into corr matrix
746  SaveCorr(covOutData, "corrOutData", dOut); // transform cov matrix into corr matrix
747 
748  // main reason to do it is to get the condition number (saved in the title)
749  // -> otherwise, it should be identical to the input `fMC` file
750  cout << __LINE__ << "\tSaving response and probability matrices" << endl;
751  RM->SetDirectory(dOut);
752  RM->Write("RM"); // title was already set earlier (corresponding to condition number)
753  RM->SetDirectory(nullptr);
754  PM->SetDirectory(dOut);
755  PM->Write();
756 
757  cout << __LINE__ << "\tGetting backfolding" << endl;
758  auto backfolding = unique_ptr<TH1>(density->GetFoldedOutput(
759  "backfold", // histogram name
760  "backfolded (detector level)", // histogram title
761  nullptr /* dist name */, nullptr /* axis steering */,
762  true /* axis binning */, true /* add bg */));
763  backfolding->SetDirectory(dOut);
764  backfolding->GetXaxis()->SetTitle("rec");
765  backfolding->Write();
766  recData->Write("recData");
767  covInData->Write();
768 
769  genMC->SetDirectory(dOut);
770  genMC->Write();
771  genMC->SetDirectory(nullptr);
772 
773  recMC->SetDirectory(dOut);
774  recMC->Write();
775  recMC->SetDirectory(nullptr);
776  covInMC->SetDirectory(dOut);
777  covInMC->Write();
778 
779  cout << __LINE__ << "\tGetting systematic uncertainties from MC and from unfolding"
780  << endl;
781 
782  // just initialising histogram to contain contribution from the limited statistics
783  // of the MC (not only from RM, but also from fake and miss histograms)
784  auto covOutUncor = Clone<TH2>(covOutData, "covOutUncor");
785  covOutUncor->SetTitle("covariance matrix at particle level coming from the MC "
786  "(RM + miss + fake)");
787  covOutUncor->Reset();
788 
789  auto dUncorr = dOut->mkdir("Uncorr");
790  dUncorr->cd();
791  cout << __LINE__ << "\tCovariance matrix from RM" << endl;
792  auto SysUncorr = unique_ptr<TH2>(density->GetEmatrixSysUncorr("RM", // name
793  "covariance matrix at particle level coming from the RM")); // title
794  SysUncorr->SetDirectory(dUncorr);
795  SysUncorr->GetXaxis()->SetTitle("gen");
796  SysUncorr->Write();
797  covOutUncor->Add(SysUncorr.get());
798  cout << __LINE__ << "\tCovariance matrices from backgrounds" << endl;
799  for (auto& f: fake) {
800  auto SysBackgroundUncorr = unique_ptr<TH2>(density->GetEmatrixSysBackgroundUncorr(
801  "fake" + f.first, // source
802  "fake" + f.first, // histogram name
803  Form("covariance matrix coming from %s", f.second->GetTitle())
804  ));
805  SysBackgroundUncorr->SetDirectory(dUncorr);
806  SysBackgroundUncorr->Write();
807  covOutUncor->Add(SysBackgroundUncorr.get());
808  }
809 
810  if (sysUncorr) {
811  cout << __LINE__ << "\tCovariance matrices from input bin-to-bin "
812  "uncorrelated systematic uncertainties." << endl;
813  auto SysBackgroundUncorr = unique_ptr<TH2>(
814  density->GetEmatrixSysBackgroundUncorr(
815  "sysUncorr", "sysUncorr",
816  "covariance matrix from trigger and prefiring uncertainties"));
817  SysBackgroundUncorr->SetDirectory(dUncorr);
818  SysBackgroundUncorr->Write();
819  covOutUncor->Add(SysBackgroundUncorr.get());
820  }
821  cout << __LINE__ << "\tCovariance matrices from inefficiencies" << endl;
822  for (const auto& m: miss) { // only using miss for the names
823  auto SysIneffUncorr = Clone<TH2>(SysUncorr, "miss" + m.first);
824  SysIneffUncorr->Reset();
825  for (int i = 1; i <= m.second->GetNbinsX(); ++i) {
826  double error = m.second->GetBinError(i) * ::abs(unf->GetBinContent(i));
827  SysIneffUncorr->SetBinContent(i, i, pow(error, 2));
828  }
829  SysIneffUncorr->SetDirectory(dUncorr);
830  SysIneffUncorr->Write();
831  covOutUncor->Add(SysIneffUncorr.get());
832  }
833  dOut->cd();
834  covOutUncor->SetDirectory(dOut);
835  covOutUncor->Write();
836 
837  auto covOutTotal = Clone<TH2>(covOutData, "covOutTotal");
838  covOutTotal->Add(covOutUncor.get());
839  covOutTotal->SetTitle("covariance matrix at particle level including "
840  "contributions from data and MC");
841  covOutTotal->SetDirectory(dOut);
842  covOutTotal->Write();
843  SaveCorr(covOutTotal, "corrOutTotal", dOut);
844 
845  cout << __LINE__ << "\tSaving estimates of fake counts" << endl;
846  for (auto& f: fake) {
847  TString name = "fake" + f.first,
848  title = f.second->GetTitle();
849  cout << __LINE__ << '\t' << name << endl;
850 
851  // count
852  f.second->SetDirectory(dOut);
853  f.second->SetTitle(title + " (corrected)");
854  f.second->Write(name);
855  f.second->SetDirectory(nullptr);
856  }
857 
858  if (var == "nominal" && applySyst) {
859  cout << __LINE__ << "\tSaving shifts from background variations "
860  "in parallel directories" << endl;
861  fOut->cd();
862  for (auto& f: fake) {
863  TString name = "fake" + f.first,
864  title = f.second->GetTitle();
865  cout << __LINE__ << '\t' << name << endl;
866 
867  // upper shift
868  {
869  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
870  name, // name of the source given in `SubtractBackground()`
871  name + SysUp, // name of new histogram
872  ""));
873  shift->Add(unf.get());
874  auto dOut = fOut->mkdir(name + SysUp, "upper shift of " + title);
875  shift->SetDirectory(dOut);
876  dOut->cd();
877  shift->GetXaxis()->SetTitle("gen");
878  shift->Write("unfold");
880  }
881 
882  // lower shift
883  {
884  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
885  name, // name of the source given in `SubtractBackground()`
886  name + SysDown, // name of new histogram
887  ""));
888  shift->Scale(-1);
889  shift->Add(unf.get());
890  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
891  shift->SetDirectory(dOut);
892  dOut->cd();
893  shift->GetXaxis()->SetTitle("gen");
894  shift->Write("unfold");
896  }
897  }
898  }
899 
900  cout << __LINE__ << "\tSaving estimates of miss counts" << endl;
901  for (auto& m: miss) {
902  TString name = "miss" + m.first,
903  title = m.second->GetTitle();
904  cout << __LINE__ << '\t' << name << endl;
905 
906  auto M = Clone<TH1>(m.second, name);
907 
908  // count
909  dOut->cd();
910  for (int i = 1; i <= unf->GetNbinsX(); ++i) {
911  double u = unf->GetBinContent(i),
912  content = m.second->GetBinContent(i) * u ,
913  error = m.second->GetBinError (i) * ::abs(u);
914  M->SetBinContent(i, content);
915  M->SetBinError (i, error );
916  }
917  M->SetDirectory(dOut);
918  M->Write();
919  M->SetDirectory(nullptr);
920 
921  if (var == "nominal" && applySyst) {
922  cout << __LINE__ << "\tSaving shifts from inefficiency variations "
923  "in parallel directories" << endl;
924  fOut->cd();
925 
926  auto missNormVar = config.get_optional<float>("unfolding.missNormVar").value_or(0.05);
927  cout << __LINE__ << "\tmissNormVar = " << missNormVar << endl;
928 
929  // upper shifts
930  {
931  auto shift = Clone<TH1>(M, name + SysUp);
932  shift->Scale(missNormVar);
933  shift->Add(unf.get());
934  auto dOut = fOut->mkdir(name + SysUp, "lower shift of " + title);
935  shift->SetDirectory(dOut);
936  dOut->cd();
937  shift->GetXaxis()->SetTitle("gen");
938  shift->Write("unfold");
940  }
941 
942  // lower shifts
943  {
944  auto shift = Clone<TH1>(M, name + SysDown);
945  shift->Scale(-1);
946  shift->Scale(missNormVar);
947  shift->Add(unf.get());
948  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
949  shift->SetDirectory(dOut);
950  dOut->cd();
951  shift->GetXaxis()->SetTitle("gen");
952  shift->Write("unfold");
954  }
955  }
956  }
957 
958  // basically, everything has now been retrieved
959  // -> now we just perform some last test
960 
961  cout << __LINE__ << "\tBottom line test" << endl;
962  {
963  int rebin = recMC->GetNbinsX()/genMC->GetNbinsX();
964 
965  cout << __LINE__ << "\tDetector level (`rebin` = " << rebin << "): ";
966  Checks::BLT(recData, covInData, recMC, rebin);
967 
968  cout << __LINE__ << "\tHadron level (with data unc. only): ";
969  Checks::BLT(unf, covOutData, genMC);
970 
971  cout << __LINE__ << "\tHadron level (with data + MC unc.): ";
972  Checks::BLT(unf, covOutTotal, genMC);
973  }
974  } // end of data variations
975  } // end of MC variations
976 
977  cout << __func__ << ' ' << slice << " stop" << endl;
978 } // end of `unfold`
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
DAS::Unfolding::Tikhonov::PrepareRegularisation
void PrepareRegularisation(TString MCvar, const char *fname, unique_ptr< MyTUnfoldDensity > &density)
Definition: unfold.cc:275
DAS::Unfolding::Checks::NegativeBins
void NegativeBins(const unique_ptr< TH1 > &unf)
Definition: unfold.cc:122
jmarExample.pt
pt
Definition: jmarExample.py:19
DAS::Unfolding::Tikhonov::RedoUnfold
void RedoUnfold(unique_ptr< MyTUnfoldDensity > &density, double newtau, TString name, TString subtitle, shared_ptr< TFile > fOut, const map< TString, unique_ptr< TH1 >> &miss)
Definition: unfold.cc:330
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
Step::def
static const char * def
Definition: Step.h:36
DAS::SysUp
const std::string SysUp
Suffix used for "up" uncertainties. Follows the Combine convention.
Definition: Format.h:8
Ntupliser_cfg.f
f
Definition: Ntupliser_cfg.py:256
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:29
DAS::Unfolding::GetVariations
std::vector< DistVariation > GetVariations(Darwin::Tools::MetaInfo &, bool=false, std::ostream &=std::cout)
Get all variations availables according to metainfo.
Definition: DistVariation.cc:199
jercExample.key
string key
Definition: jercExample.py:109
DAS::JetEnergy::getBinning
std::vector< double > getBinning(int nBins, float first, float last)
Definition: common.h:53
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
DAS::Unfolding::CorrectInefficiencies
void CorrectInefficiencies(unique_ptr< TH1 > &unf, const map< TString, unique_ptr< TH1 >> &miss)
Definition: unfold.cc:58
DAS::Unfolding
Definition: getToyCalculation.cc:37
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:99
DAS::Unfolding::CheckEntries
void CheckEntries(const unique_ptr< TH1X > &h, const unique_ptr< TH2X > &cov)
Definition: getToyCalculation.cc:58
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:264
OBS_TYPES
#define OBS_TYPES
DAS::Uncertainties::nominal
const Variation nominal
Definition: Variation.h:55
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
IF_OBS
#define IF_OBS(r, data, TYPE)
deps
static const auto deps
Definition: getToyCalculation.cc:35
DAS::Unfolding::Checks::folding
void folding(const unique_ptr< TH1 > &genIn, const unique_ptr< TH1 > &recIn, const unique_ptr< TH2 > &PMin, const map< TString, unique_ptr< TH1 >> &miss, const map< TString, unique_ptr< TH1 >> &fake)
Definition: unfold.cc:161
DAS::Unfolding::Checks::GetPM
unique_ptr< TH2 > GetPM(const unique_ptr< TH2 > &RM)
Calculate probability matrix.
Definition: unfold.cc:136
DAS::Unfolding::fillCov
void fillCov(const DistVariation &, const std::list< int > &)
Definition: DistVariation.cc:185
DAS::Unfolding::initOptionalBranch
auto initOptionalBranch(TTreeReader &reader, const char *name)
Definition: Observable.h:39
DAS::Unfolding::Checks::BLT
void BLT(const unique_ptr< TH1 > &dataDist, const unique_ptr< TH2 > &dataCov, const unique_ptr< TH1 > &MC, int rebin=1)
Definition: unfold.cc:207
DAS::SysDown
const std::string SysDown
Suffix used for "down" uncertainties. Follows the Combine convention.
Definition: Format.h:10
Darwin::Tools::UserInfo::List
TList * List(TList *mother, const char *key) const
Accest do daughter list from a mother list.
Definition: UserInfo.cc:33
Ntupliser_cfg.isMC
string isMC
Definition: Ntupliser_cfg.py:59
DAS::Unfolding::Tikhonov::ScanLcurve
double ScanLcurve(unique_ptr< MyTUnfoldDensity > &density, TDirectory *dOut, const pt::ptree &config)
Save L-curve and other curves, as well as solutions.
Definition: unfold.cc:348
jercExample.inputs
def inputs
Definition: jercExample.py:118
DAS::Unfolding::ExtractHistogram::nodeLoop
void nodeLoop(TDirectory *, TUnfoldBinning *, T, ostream &, TString="", int=1)
Loop on the nodes of a given binning.
Definition: getUnfPhysDist.cc:201
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
DAS::Unfolding::GetObservables
std::vector< Observable * > GetObservables(boost::property_tree::ptree)
Get the observables to unfold.
Definition: Observable.cc:43
DAS::Unfolding::CheckSymmetry
void CheckSymmetry(const unique_ptr< TH1X > &h, const unique_ptr< TH2X > &cov)
Definition: getToyCalculation.cc:40
DAS::Unfolding::Checks::CovMargin
void CovMargin(unique_ptr< TH1 > &h, unique_ptr< TH2 > &cov, TString name, int ibin)
Set underflow/overflow to 0.
Definition: unfold.cc:250
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83
DAS::Unfolding::Migrations::Condition
double Condition(const unique_ptr< TH2 > &RM)
Definition: unfold.cc:83
DAS::Unfolding::SaveCorr
void SaveCorr(const std::unique_ptr< TH2X > &cov, TString name, TDirectory *d)
Save correlation matrix.
Definition: toolbox.h:14