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)
 
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.

Todo:
make this function a static member of Observable
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::Flow(steering, {input}).GetInputHist<TH2>("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
45 {
46  vector<Observable *> observables;
47 #define OBS_TYPES (InclusiveJet::PtY)(DijetMass::MjjYmax)(DijetMass3D::MjjYbYs)(Rij::HTn)(Rij::Ratios)(DrellYan::ZPtY)(ZmmY::Dalitz)(ZmmY::BF)(ZmmY::DYJetsToMuMu)(ZmmY::DYJetsToTauTau)(ZmmY::TTTo2L2Nu)(ZmmY::QCD)(ZmmY::WW)(ZmmY::WZ)(ZmmY::ZZ)(MNjets::DEtaDPhi)
48 #define IF_OBS(r, data, TYPE) \
49  if (pt.find(BOOST_PP_STRINGIZE(TYPE)) != pt.not_found()) { \
50  observables.push_back(new TYPE); \
51  for (auto it = pt.begin(); it != pt.end(); ++it) { \
52  if (it->first != BOOST_PP_STRINGIZE(TYPE)) continue; \
53  pt.erase(it); \
54  break; \
55  } \
56  }
57  BOOST_PP_SEQ_FOR_EACH(IF_OBS, _, OBS_TYPES)
58 #undef OBS_TYPES
59 #undef IF_OBS
60  if (!pt.empty()) {
61  TString unknown;
62  for (auto it = pt.begin(); it != pt.end(); ++it) {
63  unknown += it->first;
64  unknown += ' ';
65  }
66  BOOST_THROW_EXCEPTION(invalid_argument(Form("Unknown observable(s): %s", unknown.Data())));
67  }
68  return observables;
69 }

◆ 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(std::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.

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, inputs);
45  auto tIn = flow.GetInputTree(slice);
46  auto [fOut, tOut] = flow.GetOutput(output);
47 
48  DT::MetaInfo metainfo(tOut);
49  metainfo.Check(config);
50  bool isMC = metainfo.Get<bool>("flags", "isMC");
51  int R = metainfo.Get<int>("flags", "R");
52 
53  cout << "Setting observables" << endl;
54  Observable::isMC = isMC;
55  Observable::maxDR = R / 20.;
56  auto pt_obs = config.get_child("unfolding.observables");
57  vector<Observable *> observables = GetObservables(pt_obs);
58  vector<unique_ptr<Filler>> fillers;
59  for (const auto obs: observables)
60  fillers.push_back(obs->getFiller(flow));
61  if (observables.empty())
62  BOOST_THROW_EXCEPTION( invalid_argument("You should provide at least one observable "
63  "(e.g. `InclusiveJet::PtY`)") );
64 
65  cout << "Setting variations" << endl;
66  DistVariation::isMC = isMC;
67  DistVariation::genBinning = new TUnfoldBinning("gen");
68  DistVariation::recBinning = new TUnfoldBinning("rec");
69  for (Observable * obs: observables) {
70  DistVariation::genBinning->AddBinning(obs->genBinning);
71  DistVariation::recBinning->AddBinning(obs->recBinning);
72  }
73 
74  [[ maybe_unused ]]
75  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
76 
77  vector<DistVariation> variations = GetVariations(metainfo,
78  steering & DT::syst, cout);
79 
80  for (DT::Looper looper(tIn); looper(); ++looper) {
81  for (DistVariation& v: variations) {
82 
83  list<int> binIDs; // collect filled bins to calculate covariance
84  for (auto& filler: fillers)
85  binIDs.merge(filler->fillRec(v));
86  fillCov(v, binIDs);
87 
88  if (!isMC) continue;
89 
90  // run matching (possibly different for each observable)
91  for (auto& filler: fillers)
92  filler->match();
93 
94  // fill gen, RM, fake, miss
95  for (auto& filler: fillers)
96  filler->fillMC(v);
97  }
98  }
99 
100  metainfo.Set<bool>("git", "complete", true);
101 
102  fOut->cd();
103  for (DistVariation& v: variations)
104  v.Write(fOut);
105 
106  cout << __func__ << ' ' << slice << " stop" << endl;
107 }

◆ 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  variations.emplace_back("nominal","nominal");
207  if (!applySyst) return variations;
208 
209  const TList * groupList = metainfo.List("variations");
210  for (const TObject * group: *groupList) {
211  const auto groupContents = dynamic_cast<const TList*>(group);
212  if (!groupContents)
213  BOOST_THROW_EXCEPTION( invalid_argument(group->GetName()) );
214  size_t i = 0;
215  for (const TObject * obj: *groupContents) {
216  const auto string = dynamic_cast<const TObjString*>(obj);
217  if (!string)
218  BOOST_THROW_EXCEPTION( invalid_argument(obj->GetName()) );
219  cout << group->GetName() << " " << string->GetString() << " " << (i+1) << endl;
220  variations.emplace_back(group->GetName(), string->GetString(), ++i);
221  }
222  }
223 
224  return variations;
225 }

◆ 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:
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] = std::move(h);
513  }
514  else {
515  name.ReplaceAll("miss", "");
516  cout << __LINE__ << "\tmiss\t" << name << endl;
517  miss[name] = std::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"));
523  if (!RM)
524  BOOST_THROW_EXCEPTION( DE::BadInput("Missing response "
525  "matrix", *dirMC) );
526  RM->SetDirectory(nullptr);
527 
528  auto genMC = unique_ptr<TH1>(dirMC->Get<TH1>("gen")),
529  recMC = unique_ptr<TH1>(dirMC->Get<TH1>("rec"));
530  if (!genMC)
531  BOOST_THROW_EXCEPTION( DE::BadInput("Missing gen "
532  "histogram", *dirMC) );
533  genMC->SetName("genMC");
534  recMC->SetName("recMC");
535  genMC->SetDirectory(nullptr);
536  recMC->SetDirectory(nullptr);
537 
538  unique_ptr<TH2> PM = Checks::GetPM(RM);
539  PM->SetDirectory(nullptr);
540  Checks::folding(genMC, recMC, PM, miss, fake);
541 
542  cout << __LINE__ << "\tCalculating miss and fake rates" << endl;
543 
544  // we want a fake/miss >> RATE << (i.e. a value between 0 and 1)
545  // -> we will use it later to estimate the fake/miss contribution
546  // in the data using a bin-by-bin method.
547  // currently: the elements of `fake`/`miss` are currently fake/miss *counts*
548  // (including model dependence)
549  unique_ptr<TH1> genMC_noMiss = Clone(genMC, "genMC_noMiss");
550  for (auto& f: fake) f.second ->Divide(f.second.get(), recMC.get(), 1, 1, "b");
551  for (auto& m: miss) genMC_noMiss->Add (m.second.get(), -1);
552  for (auto& m: miss) m.second ->Divide(m.second.get(), genMC_noMiss.get(), 1, 1, "b");
553  // now: the elements of `fake`/`miss` are currently fake/miss *rates*
554 
555  // now we are done with retrieving objects from fMC, and we move to fData
556 
557  fData->cd();
558 
559  // now we estimate the fake contribution in data with bin-by-bin method
560  {
561  auto rec = unique_ptr<TH1>(fData->Get<TH1>("nominal/rec"));
562  // reminder: the elements of `fake` are currently fake *rates*
563  for (auto& f: fake) f.second->Multiply(rec.get());
564  // now: the elements of `fake` are currently fake *counts* (corrected to data)
565  }
566  // (similar thing will be done with the miss contribution *after* the unfolding)
567 
568  // reminder: now we have absolute fake count BUT miss rate
569 
570  cout << __LINE__ << "\tExtracting data variations and running unfolding" << endl;
571 
572  // here we loop again on the entries of a file
573  for (const auto&& objData: *(fData->GetListOfKeys())) {
574 
575  auto keyData = dynamic_cast<TKey*>(objData);
576  const TString classname = keyData->ReadObj()->ClassName();
578  //if (classname == "TTree") {
579  // unique_ptr<TTree> tree = ReadObj<TTree>(keyData);
580  // auto List = tree->GetUserInfo();
581  // fOut->cd();
582  // List->Write("Data");
583  // fMC->cd();
584  //}
585  if (!classname.Contains("TDirectory")) continue;
586 
587  const TString DATAvar = objData->GetName();
588  if (DATAvar != "nominal") {
589  if (!applySyst) continue;
590  if (MCvar != "nominal") continue;
591  }
592  cout << __LINE__ << "\tvariation: " << DATAvar << endl;
593  TString var = MCvar == "nominal" ? DATAvar : MCvar;
594 
595  if (var != "nominal" && MCvar == DATAvar)
596  cerr << orange << "The same variation has been found "
597  "in both data and MC.\n" << def;
598 
599  {
600  static int ivar = 0; // used for parallelisation
601  ++ivar;
602  if (slice.second != ivar % slice.first) continue;
603  }
604 
605  auto dirData = dynamic_cast<TDirectory*>(keyData->ReadObj());
606  dirData->cd();
607 
608  auto recData = unique_ptr<TH1>(dirData->Get<TH1>("rec"));
609  recData->SetName("recData");
610  recData->SetTitle("detector level (" + var + ")");
611  recData->GetXaxis()->SetTitle("rec");
612  recData->SetDirectory(nullptr);
613  auto covInData = unique_ptr<TH2>(dirData->Get<TH2>("cov"));
614  covInData->SetName("covInData");
615  covInData->SetTitle("covariance matrix at detector level (" + var + ")");
616  covInData->GetXaxis()->SetTitle("rec");
617  covInData->GetYaxis()->SetTitle("rec");
618  covInData->SetDirectory(nullptr);
619  auto covInMC = unique_ptr<TH2>(dirMC->Get<TH2>("cov"));
620  covInMC->SetName("covInMC");
621  covInMC->SetTitle("covariance matrix at detector level (" + var + ")");
622  covInMC->GetXaxis()->SetTitle("rec");
623  covInMC->GetYaxis()->SetTitle("rec");
624  covInMC->SetDirectory(nullptr);
625 
626  cout << __LINE__ << '\t' << "\e[1m" << DATAvar << "\e[0m" << endl;
627 
628  if (DATAvar.Contains("nominal")) { // if data variation is nominal,
629  // but the condition number will
630  // be calculated for each variation of MC
631  auto condition = Migrations::Condition(RM);
632  RM->SetTitle(Form("%.1f",condition)); // keep track of the condition number
633  }
634 
635  // underflow and overflow are expected to be empty (else, set to 0 with a warning)
636  Checks::CovMargin(recData, covInData, "underflow", 0);
637  Checks::CovMargin(recData, covInData, "overflow", recData->GetNbinsX()+1);
638 
639  // so far, still no unfolding done
640 
641  cout << __LINE__ << "\tSetting up density & input (checking "
642  "and comparing the bin contents of data and MC)" << endl;
643 
644  auto density = make_unique<MyTUnfoldDensity>(
645  RM.get(), // reminder: this is varied by the first loop
646  TUnfold::kHistMapOutputHoriz, // truth level on x-axis of the response matrix
647  TUnfold::kRegModeNone, // do not use *predefined* regularisation schemes
648  TUnfold::kEConstraintNone, // no constraint on area
649  TUnfoldDensity::kDensityModeBinWidthAndUser); // no scale factors
650 
651  // reminder: the epsilon is the smallest possible number
652  density->SetEpsMatrix(1e-100);
653 
654  int status = density->SetInput(
655  recData.get(), // reconstructed level
656  0, // bias factor (only relevant for Tikhonov regularisation)
657  0, // for bins with zero error, this number defines 1/error
658  covInData.get()); // covariance matrix
659 
660  if (status > 0)
661  cerr << __LINE__ << '\t' << orange << "unfold status = " << status << '\n' << def;
662 
663  // now, we should have (nearly) everything settled for the unfolding...
664 
665  cout << __LINE__ << "\tSubtract background" << endl;
666  auto fakeNormVar = config.get_optional<float>("unfolding.fakeNormVar").value_or(0.05);
667  cout << __LINE__ << "\tfakeNormVar = " << fakeNormVar << endl;
668  for (auto& f: fake) { // fake now contains esimations of the fake contributions in data
669  cout << __LINE__ << '\t' << f.first << endl;
670  density->SubtractBackground(
671  f.second.get(), // 1D histogram in rec-level binning (with its own stat use
672  // --> we will use later `GetEmatrixSysBackgroundUncorr()`
673  "fake" + f.first, // name
674  1.0, // scale factor
675  fakeNormVar); // factor to vary the normalisation
676  // --> we will use later `GetDeltaSysBackgroundScale()`
677  // the lines just above are equivalent to three things:
678  // - `h->Add(f.second, -1);`, i.e. subtract the fake contribution
679  // from the rec level
680  // - it propagates the contribution from the limited statistics
681  // to the covariance matrix used in the unfolding
682  // - TUnfold provides systematic variations
683  }
684 
685  if (sysUncorr) {
686  cout << __LINE__ << "\tAdding bin-to-bin uncertainty" << endl;
687 
688  if (!fs::exists(*sysUncorr))
689  BOOST_THROW_EXCEPTION(fs::filesystem_error("File with systematic "
690  "bin-to-bin uncorrelated uncertainties was not found!",
691  make_error_code(errc::no_such_file_or_directory)));
692 
693  auto fSys = make_unique<TFile>(sysUncorr->c_str(), "READ");
694  auto hSysUncorr = unique_ptr<TH1>( fSys->GetDirectory(DATAvar)
695  ->Get<TH1>("sysUncorr1D") );
696  density->SubtractBackground(hSysUncorr.get(), "sysUncorr", 1.0, 0.);
697  }
698 
699  fOut->cd();
700  if (fOut->GetDirectory(var)) continue;
701  auto dOut = fOut->mkdir(var);
702  dOut->cd();
703 
704  cout << __LINE__ << "\tUnfolding" << endl;
705  static auto tau = config.get_optional<float>("unfolding.TUnfold.regularisation.tau");
706  if (Lmatrix) {
707  cout << __LINE__ << "\tFetching L-matrix" << endl;
708 
709  if (!fs::exists(*Lmatrix))
710  BOOST_THROW_EXCEPTION(fs::filesystem_error("File with L-matrix not found!",
711  make_error_code(errc::no_such_file_or_directory)));
712 
713  Tikhonov::PrepareRegularisation("nominal", Lmatrix->c_str(), density);
714  if (!tau && var == "nominal") {
715 
716  tau = Tikhonov::ScanLcurve(density, dOut, config);
717 
718  if (applySyst) {
719  cout << __LINE__ << "\tSaving shifts for variations of tau parameter" << endl;
720  Tikhonov::RedoUnfold(density, *tau*2, "Regularisation" + SysUp,
721  "upper variation of regularisation", fOut, miss);
722  Tikhonov::RedoUnfold(density, *tau/2, "Regularisation" + SysDown,
723  "lower variation of regularisation", fOut, miss);
724  }
725  }
726  } // end of regularise
727  cout << __LINE__ << "\ttau = " << tau << endl;
728  density->DoUnfold(tau.value_or(0));
729 
730  // now that the unfolding is done, retrieve the result (use TUnfold)
731  // and apply the correction for miss entries (by hand)
732 
733  cout << __LINE__ << "\tGetting output" << endl;
734  auto unf = unique_ptr<TH1>(density->GetOutput("unfold", "particle level (" + var + ")"));
735  Checks::NegativeBins(unf);
736  CorrectInefficiencies(unf, miss);
737  dOut->cd();
738  unf->SetDirectory(dOut);
739  unf->GetXaxis()->SetTitle("gen");
740  unf->Write();
741 
745 
746  cout << __LINE__ << "\tCovariance matrix from measurement" << endl;
747  auto covOutData = unique_ptr<TH2>(density->GetEmatrixInput("covOutData", // name
748  "covariance matrix coming from input data distribution (" + var + ")")); // title
749  covOutData->SetDirectory(dOut);
750  covOutData->GetXaxis()->SetTitle("gen");
751  covOutData->GetYaxis()->SetTitle("gen");
752  covOutData->Write();
753 
754  cout << __LINE__ << "\tSaving correlation matrix from measurement" << endl;
755  SaveCorr(covInData , "corrInData" , dOut); // transform cov matrix into corr matrix
756  SaveCorr(covOutData, "corrOutData", dOut); // transform cov matrix into corr matrix
757 
758  // main reason to do it is to get the condition number (saved in the title)
759  // -> otherwise, it should be identical to the input `fMC` file
760  cout << __LINE__ << "\tSaving response and probability matrices" << endl;
761  RM->SetDirectory(dOut);
762  RM->Write("RM"); // title was already set earlier (corresponding to condition number)
763  RM->SetDirectory(nullptr);
764  PM->SetDirectory(dOut);
765  PM->Write();
766 
767  cout << __LINE__ << "\tGetting backfolding" << endl;
768  auto backfolding = unique_ptr<TH1>(density->GetFoldedOutput(
769  "backfold", // histogram name
770  "backfolded (detector level)", // histogram title
771  nullptr /* dist name */, nullptr /* axis steering */,
772  true /* axis binning */, true /* add bg */));
773  backfolding->SetDirectory(dOut);
774  backfolding->GetXaxis()->SetTitle("rec");
775  backfolding->Write();
776  recData->Write("recData");
777  covInData->Write();
778 
779  genMC->SetDirectory(dOut);
780  genMC->Write();
781  genMC->SetDirectory(nullptr);
782 
783  recMC->SetDirectory(dOut);
784  recMC->Write();
785  recMC->SetDirectory(nullptr);
786  covInMC->SetDirectory(dOut);
787  covInMC->Write();
788 
789  cout << __LINE__ << "\tGetting systematic uncertainties from MC and from unfolding"
790  << endl;
791 
792  // just initialising histogram to contain contribution from the limited statistics
793  // of the MC (not only from RM, but also from fake and miss histograms)
794  auto covOutUncor = Clone<TH2>(covOutData, "covOutUncor");
795  covOutUncor->SetTitle("covariance matrix at particle level coming from the MC "
796  "(RM + miss + fake)");
797  covOutUncor->Reset();
798 
799  auto dUncorr = dOut->mkdir("Uncorr");
800  dUncorr->cd();
801  cout << __LINE__ << "\tCovariance matrix from RM" << endl;
802  auto SysUncorr = unique_ptr<TH2>(density->GetEmatrixSysUncorr("RM", // name
803  "covariance matrix at particle level coming from the RM")); // title
804  SysUncorr->SetDirectory(dUncorr);
805  SysUncorr->GetXaxis()->SetTitle("gen");
806  SysUncorr->Write();
807  covOutUncor->Add(SysUncorr.get());
808  cout << __LINE__ << "\tCovariance matrices from backgrounds" << endl;
809  for (auto& f: fake) {
810  auto SysBackgroundUncorr = unique_ptr<TH2>(density->GetEmatrixSysBackgroundUncorr(
811  "fake" + f.first, // source
812  "fake" + f.first, // histogram name
813  Form("covariance matrix coming from %s", f.second->GetTitle())
814  ));
815  SysBackgroundUncorr->SetDirectory(dUncorr);
816  SysBackgroundUncorr->Write();
817  covOutUncor->Add(SysBackgroundUncorr.get());
818  }
819 
820  if (sysUncorr) {
821  cout << __LINE__ << "\tCovariance matrices from input bin-to-bin "
822  "uncorrelated systematic uncertainties." << endl;
823  auto SysBackgroundUncorr = unique_ptr<TH2>(
824  density->GetEmatrixSysBackgroundUncorr(
825  "sysUncorr", "sysUncorr",
826  "covariance matrix from trigger and prefiring uncertainties"));
827  SysBackgroundUncorr->SetDirectory(dUncorr);
828  SysBackgroundUncorr->Write();
829  covOutUncor->Add(SysBackgroundUncorr.get());
830  }
831  cout << __LINE__ << "\tCovariance matrices from inefficiencies" << endl;
832  for (const auto& m: miss) { // only using miss for the names
833  auto SysIneffUncorr = Clone<TH2>(SysUncorr, "miss" + m.first);
834  SysIneffUncorr->Reset();
835  for (int i = 1; i <= m.second->GetNbinsX(); ++i) {
836  double error = m.second->GetBinError(i) * ::abs(unf->GetBinContent(i));
837  SysIneffUncorr->SetBinContent(i, i, pow(error, 2));
838  }
839  SysIneffUncorr->SetDirectory(dUncorr);
840  SysIneffUncorr->Write();
841  covOutUncor->Add(SysIneffUncorr.get());
842  }
843  dOut->cd();
844  covOutUncor->SetDirectory(dOut);
845  covOutUncor->Write();
846 
847  auto covOutTotal = Clone<TH2>(covOutData, "covOutTotal");
848  covOutTotal->Add(covOutUncor.get());
849  covOutTotal->SetTitle("covariance matrix at particle level including "
850  "contributions from data and MC");
851  covOutTotal->SetDirectory(dOut);
852  covOutTotal->Write();
853  SaveCorr(covOutTotal, "corrOutTotal", dOut);
854 
855  cout << __LINE__ << "\tSaving estimates of fake counts" << endl;
856  for (auto& f: fake) {
857  TString name = "fake" + f.first,
858  title = f.second->GetTitle();
859  cout << __LINE__ << '\t' << name << endl;
860 
861  // count
862  f.second->SetDirectory(dOut);
863  f.second->SetTitle(title + " (corrected)");
864  f.second->Write(name);
865  f.second->SetDirectory(nullptr);
866  }
867 
868  if (var == "nominal" && applySyst) {
869  cout << __LINE__ << "\tSaving shifts from background variations "
870  "in parallel directories" << endl;
871  fOut->cd();
872  for (auto& f: fake) {
873  TString name = "fake" + f.first,
874  title = f.second->GetTitle();
875  cout << __LINE__ << '\t' << name << endl;
876 
877  // upper shift
878  {
879  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
880  name, // name of the source given in `SubtractBackground()`
881  name + SysUp, // name of new histogram
882  ""));
883  shift->Add(unf.get());
884  auto dOut = fOut->mkdir(name + SysUp, "upper shift of " + title);
885  shift->SetDirectory(dOut);
886  dOut->cd();
887  shift->GetXaxis()->SetTitle("gen");
888  shift->Write("unfold");
890  }
891 
892  // lower shift
893  {
894  auto shift = unique_ptr<TH1>(density->GetDeltaSysBackgroundScale(
895  name, // name of the source given in `SubtractBackground()`
896  name + SysDown, // name of new histogram
897  ""));
898  shift->Scale(-1);
899  shift->Add(unf.get());
900  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
901  shift->SetDirectory(dOut);
902  dOut->cd();
903  shift->GetXaxis()->SetTitle("gen");
904  shift->Write("unfold");
906  }
907  }
908  }
909 
910  cout << __LINE__ << "\tSaving estimates of miss counts" << endl;
911  for (auto& m: miss) {
912  TString name = "miss" + m.first,
913  title = m.second->GetTitle();
914  cout << __LINE__ << '\t' << name << endl;
915 
916  auto M = Clone<TH1>(m.second, name);
917 
918  // count
919  dOut->cd();
920  for (int i = 1; i <= unf->GetNbinsX(); ++i) {
921  double u = unf->GetBinContent(i),
922  content = m.second->GetBinContent(i) * u ,
923  error = m.second->GetBinError (i) * ::abs(u);
924  M->SetBinContent(i, content);
925  M->SetBinError (i, error );
926  }
927  M->SetDirectory(dOut);
928  M->Write();
929  M->SetDirectory(nullptr);
930 
931  if (var == "nominal" && applySyst) {
932  cout << __LINE__ << "\tSaving shifts from inefficiency variations "
933  "in parallel directories" << endl;
934  fOut->cd();
935 
936  auto missNormVar = config.get_optional<float>("unfolding.missNormVar").value_or(0.05);
937  cout << __LINE__ << "\tmissNormVar = " << missNormVar << endl;
938 
939  // upper shifts
940  {
941  auto shift = Clone<TH1>(M, name + SysUp);
942  shift->Scale(missNormVar);
943  shift->Add(unf.get());
944  auto dOut = fOut->mkdir(name + SysUp, "lower shift of " + title);
945  shift->SetDirectory(dOut);
946  dOut->cd();
947  shift->GetXaxis()->SetTitle("gen");
948  shift->Write("unfold");
950  }
951 
952  // lower shifts
953  {
954  auto shift = Clone<TH1>(M, name + SysDown);
955  shift->Scale(-1);
956  shift->Scale(missNormVar);
957  shift->Add(unf.get());
958  auto dOut = fOut->mkdir(name + SysDown, "lower shift of " + title);
959  shift->SetDirectory(dOut);
960  dOut->cd();
961  shift->GetXaxis()->SetTitle("gen");
962  shift->Write("unfold");
964  }
965  }
966  }
967 
968  // basically, everything has now been retrieved
969  // -> now we just perform some last test
970 
971  cout << __LINE__ << "\tBottom line test" << endl;
972  {
973  int rebin = recMC->GetNbinsX()/genMC->GetNbinsX();
974 
975  cout << __LINE__ << "\tDetector level (`rebin` = " << rebin << "): ";
976  Checks::BLT(recData, covInData, recMC, rebin);
977 
978  cout << __LINE__ << "\tHadron level (with data unc. only): ";
979  Checks::BLT(unf, covOutData, genMC);
980 
981  cout << __LINE__ << "\tHadron level (with data + MC unc.): ";
982  Checks::BLT(unf, covOutTotal, genMC);
983  }
984  } // end of data variations
985  } // end of MC variations
986 
987  cout << __func__ << ' ' << slice << " stop" << endl;
988 } // 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:101
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:78
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::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:22
Ntupliser_cfg.f
f
Definition: Ntupliser_cfg.py:309
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:55
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:141
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:317
OBS_TYPES
#define OBS_TYPES
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::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
dictionary isMC
Definition: Ntupliser_cfg.py:61
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:44
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