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  unique_ptr<TFile> fOut(DT_GetOutput(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  unique_ptr<TFile> fOut(DT_GetOutput(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  unique_ptr<TChain> tIn = DT::GetChain(inputs);
45  unique_ptr<TFile> fOut(DT_GetOutput(output));
46  auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
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.get());
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.get());
113 
114  fOut->cd();
115  tOut->Write();
116 
117  cout << __func__ << ' ' << slice << " stop" << endl;
118 }

◆ 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  fOut(DT_GetOutput(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  unique_ptr<TFile> fOut(DT_GetOutput(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  auto dOut = fOut->mkdir(var);
691  dOut->cd();
692 
693  cout << __LINE__ << "\tUnfolding" << endl;
694  static auto tau = config.get_optional<float>("unfolding.TUnfold.regularisation.tau");
695  if (Lmatrix) {
696  cout << __LINE__ << "\tFetching L-matrix" << endl;
697 
698  if (!fs::exists(*Lmatrix))
699  BOOST_THROW_EXCEPTION(fs::filesystem_error("File with L-matrix not found!",
700  make_error_code(errc::no_such_file_or_directory)));
701 
702  Tikhonov::PrepareRegularisation("nominal", Lmatrix->c_str(), density);
703  if (!tau && var == "nominal") {
704 
705  tau = Tikhonov::ScanLcurve(density, dOut, config);
706 
707  if (applySyst) {
708  cout << __LINE__ << "\tSaving shifts for variations of tau parameter" << endl;
709  Tikhonov::RedoUnfold(density, *tau*2, "Regularisation" + SysUp,
710  "upper variation of regularisation", fOut, miss);
711  Tikhonov::RedoUnfold(density, *tau/2, "Regularisation" + SysDown,
712  "lower variation of regularisation", fOut, miss);
713  }
714  }
715  } // end of regularise
716  cout << __LINE__ << "\ttau = " << tau << endl;
717  density->DoUnfold(tau.value_or(0));
718 
719  // now that the unfolding is done, retrieve the result (use TUnfold)
720  // and apply the correction for miss entries (by hand)
721 
722  cout << __LINE__ << "\tGetting output" << endl;
723  auto unf = unique_ptr<TH1>(density->GetOutput("unfold", "particle level (" + var + ")"));
724  Checks::NegativeBins(unf);
725  CorrectInefficiencies(unf, miss);
726  dOut->cd();
727  unf->SetDirectory(dOut);
728  unf->GetXaxis()->SetTitle("gen");
729  unf->Write();
730 
734 
735  cout << __LINE__ << "\tCovariance matrix from measurement" << endl;
736  auto covOutData = unique_ptr<TH2>(density->GetEmatrixInput("covOutData", // name
737  "covariance matrix coming from input data distribution (" + var + ")")); // title
738  covOutData->SetDirectory(dOut);
739  covOutData->GetXaxis()->SetTitle("gen");
740  covOutData->GetYaxis()->SetTitle("gen");
741  covOutData->Write();
742 
743  cout << __LINE__ << "\tSaving correlation matrix from measurement" << endl;
744  SaveCorr(covInData , "corrInData" , dOut); // transform cov matrix into corr matrix
745  SaveCorr(covOutData, "corrOutData", dOut); // transform cov matrix into corr matrix
746 
747  // main reason to do it is to get the condition number (saved in the title)
748  // -> otherwise, it should be identical to the input `fMC` file
749  cout << __LINE__ << "\tSaving response and probability matrices" << endl;
750  RM->SetDirectory(dOut);
751  RM->Write("RM"); // title was already set earlier (corresponding to condition number)
752  RM->SetDirectory(nullptr);
753  PM->SetDirectory(dOut);
754  PM->Write();
755 
756  cout << __LINE__ << "\tGetting backfolding" << endl;
757  auto backfolding = unique_ptr<TH1>(density->GetFoldedOutput(
758  "backfold", // histogram name
759  "backfolded (detector level)", // histogram title
760  nullptr /* dist name */, nullptr /* axis steering */,
761  true /* axis binning */, true /* add bg */));
762  backfolding->SetDirectory(dOut);
763  backfolding->GetXaxis()->SetTitle("rec");
764  backfolding->Write();
765  recData->Write("recData");
766  covInData->Write();
767 
768  genMC->SetDirectory(dOut);
769  genMC->Write();
770  genMC->SetDirectory(nullptr);
771 
772  recMC->SetDirectory(dOut);
773  recMC->Write();
774  recMC->SetDirectory(nullptr);
775  covInMC->SetDirectory(dOut);
776  covInMC->Write();
777 
778  cout << __LINE__ << "\tGetting systematic uncertainties from MC and from unfolding"
779  << endl;
780 
781  // just initialising histogram to contain contribution from the limited statistics
782  // of the MC (not only from RM, but also from fake and miss histograms)
783  auto covOutUncor = Clone<TH2>(covOutData, "covOutUncor");
784  covOutUncor->SetTitle("covariance matrix at particle level coming from the MC "
785  "(RM + miss + fake)");
786  covOutUncor->Reset();
787 
788  auto dUncorr = dOut->mkdir("Uncorr");
789  dUncorr->cd();
790  cout << __LINE__ << "\tCovariance matrix from RM" << endl;
791  auto SysUncorr = unique_ptr<TH2>(density->GetEmatrixSysUncorr("RM", // name
792  "covariance matrix at particle level coming from the RM")); // title
793  SysUncorr->SetDirectory(dUncorr);
794  SysUncorr->GetXaxis()->SetTitle("gen");
795  SysUncorr->Write();
796  covOutUncor->Add(SysUncorr.get());
797  cout << __LINE__ << "\tCovariance matrices from backgrounds" << endl;
798  for (auto& f: fake) {
799  auto SysBackgroundUncorr = unique_ptr<TH2>(density->GetEmatrixSysBackgroundUncorr(
800  "fake" + f.first, // source
801  "fake" + f.first, // histogram name
802  Form("covariance matrix coming from %s", f.second->GetTitle())
803  ));
804  SysBackgroundUncorr->SetDirectory(dUncorr);
805  SysBackgroundUncorr->Write();
806  covOutUncor->Add(SysBackgroundUncorr.get());
807  }
808 
809  if (sysUncorr) {
810  cout << __LINE__ << "\tCovariance matrices from input bin-to-bin "
811  "uncorrelated systematic uncertainties." << endl;
812  auto SysBackgroundUncorr = unique_ptr<TH2>(
813  density->GetEmatrixSysBackgroundUncorr(
814  "sysUncorr", "sysUncorr",
815  "covariance matrix from trigger and prefiring uncertainties"));
816  SysBackgroundUncorr->SetDirectory(dUncorr);
817  SysBackgroundUncorr->Write();
818  covOutUncor->Add(SysBackgroundUncorr.get());
819  }
820  cout << __LINE__ << "\tCovariance matrices from inefficiencies" << endl;
821  for (const auto& m: miss) { // only using miss for the names
822  auto SysIneffUncorr = Clone<TH2>(SysUncorr, "miss" + m.first);
823  SysIneffUncorr->Reset();
824  for (int i = 1; i <= m.second->GetNbinsX(); ++i) {
825  double error = m.second->GetBinError(i) * ::abs(unf->GetBinContent(i));
826  SysIneffUncorr->SetBinContent(i, i, pow(error, 2));
827  }
828  SysIneffUncorr->SetDirectory(dUncorr);
829  SysIneffUncorr->Write();
830  covOutUncor->Add(SysIneffUncorr.get());
831  }
832  dOut->cd();
833  covOutUncor->SetDirectory(dOut);
834  covOutUncor->Write();
835 
836  auto covOutTotal = Clone<TH2>(covOutData, "covOutTotal");
837  covOutTotal->Add(covOutUncor.get());
838  covOutTotal->SetTitle("covariance matrix at particle level including "
839  "contributions from data and MC");
840  covOutTotal->SetDirectory(dOut);
841  covOutTotal->Write();
842  SaveCorr(covOutTotal, "corrOutTotal", dOut);
843 
844  cout << __LINE__ << "\tSaving estimates of fake counts" << endl;
845  for (auto& f: fake) {
846  TString name = "fake" + f.first,
847  title = f.second->GetTitle();
848  cout << __LINE__ << '\t' << name << endl;
849 
850  // count
851  f.second->SetDirectory(dOut);
852  f.second->SetTitle(title + " (corrected)");
853  f.second->Write(name);
854  f.second->SetDirectory(nullptr);
855  }
856 
857  if (var == "nominal" && applySyst) {
858  cout << __LINE__ << "\tSaving shifts from background variations "
859  "in parallel directories" << endl;
860  fOut->cd();
861  for (auto& f: fake) {
862  TString name = "fake" + f.first,
863  title = f.second->GetTitle();
864  cout << __LINE__ << '\t' << name << endl;
865 
866  // upper shift
867  {
868  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
869  name, // name of the source given in `SubtractBackground()`
870  name + SysUp, // name of new histogram
871  ""));
872  shift->Add(unf.get());
873  auto dOut = fOut->mkdir(name + SysUp, "upper shift of " + title);
874  shift->SetDirectory(dOut);
875  dOut->cd();
876  shift->GetXaxis()->SetTitle("gen");
877  shift->Write("unfold");
879  }
880 
881  // lower shift
882  {
883  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
884  name, // name of the source given in `SubtractBackground()`
885  name + SysDown, // name of new histogram
886  ""));
887  shift->Scale(-1);
888  shift->Add(unf.get());
889  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
890  shift->SetDirectory(dOut);
891  dOut->cd();
892  shift->GetXaxis()->SetTitle("gen");
893  shift->Write("unfold");
895  }
896  }
897  }
898 
899  cout << __LINE__ << "\tSaving estimates of miss counts" << endl;
900  for (auto& m: miss) {
901  TString name = "miss" + m.first,
902  title = m.second->GetTitle();
903  cout << __LINE__ << '\t' << name << endl;
904 
905  auto M = Clone<TH1>(m.second, name);
906 
907  // count
908  dOut->cd();
909  for (int i = 1; i <= unf->GetNbinsX(); ++i) {
910  double u = unf->GetBinContent(i),
911  content = m.second->GetBinContent(i) * u ,
912  error = m.second->GetBinError (i) * ::abs(u);
913  M->SetBinContent(i, content);
914  M->SetBinError (i, error );
915  }
916  M->SetDirectory(dOut);
917  M->Write();
918  M->SetDirectory(nullptr);
919 
920  if (var == "nominal" && applySyst) {
921  cout << __LINE__ << "\tSaving shifts from inefficiency variations "
922  "in parallel directories" << endl;
923  fOut->cd();
924 
925  auto missNormVar = config.get_optional<float>("unfolding.missNormVar").value_or(0.05);
926  cout << __LINE__ << "\tmissNormVar = " << missNormVar << endl;
927 
928  // upper shifts
929  {
930  auto shift = Clone<TH1>(M, name + SysUp);
931  shift->Scale(missNormVar);
932  shift->Add(unf.get());
933  auto dOut = fOut->mkdir(name + SysUp, "lower shift of " + title);
934  shift->SetDirectory(dOut);
935  dOut->cd();
936  shift->GetXaxis()->SetTitle("gen");
937  shift->Write("unfold");
939  }
940 
941  // lower shifts
942  {
943  auto shift = Clone<TH1>(M, name + SysDown);
944  shift->Scale(-1);
945  shift->Scale(missNormVar);
946  shift->Add(unf.get());
947  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
948  shift->SetDirectory(dOut);
949  dOut->cd();
950  shift->GetXaxis()->SetTitle("gen");
951  shift->Write("unfold");
953  }
954  }
955  }
956 
957  // basically, everything has now been retrieved
958  // -> now we just perform some last test
959 
960  cout << __LINE__ << "\tBottom line test" << endl;
961  {
962  int rebin = recMC->GetNbinsX()/genMC->GetNbinsX();
963 
964  cout << __LINE__ << "\tDetector level (`rebin` = " << rebin << "): ";
965  Checks::BLT(recData, covInData, recMC, rebin);
966 
967  cout << __LINE__ << "\tHadron level (with data unc. only): ";
968  Checks::BLT(unf, covOutData, genMC);
969 
970  cout << __LINE__ << "\tHadron level (with data + MC unc.): ";
971  Checks::BLT(unf, covOutTotal, genMC);
972  }
973  } // end of data variations
974  } // end of MC variations
975 
976  cout << __func__ << ' ' << slice << " stop" << endl;
977 } // 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
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:30
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
Darwin::Tools::GetChain
std::unique_ptr< TChain > GetChain(std::vector< std::filesystem::path > inputs, const char *name="events")
Load chain from a list of files.
Definition: FileUtils.cc:67
Ntupliser_cfg.f
f
Definition: Ntupliser_cfg.py:255
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:49
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:65
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
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:263
OBS_TYPES
#define OBS_TYPES
DAS::Uncertainties::nominal
const Variation nominal
Definition: Variation.h:109
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::RedoUnfold
void RedoUnfold(unique_ptr< MyTUnfoldDensity > &density, double newtau, TString name, TString subtitle, unique_ptr< TFile > &fOut, const map< TString, unique_ptr< TH1 >> &miss)
Definition: unfold.cc:330
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
DT_GetOutput
#define DT_GetOutput(output)
Definition: FileUtils.h:96
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:68
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