ATLAS Offline Software
TElectronEfficiencyCorrectionTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3  */
4 
15 // This class header
17 // STL includes
18 #include <cmath>
19 #include <iostream>
20 #include <memory>
21 #include <cstring>
22 // ROOT includes
23 #include "TClass.h"
24 #include "TFile.h"
25 #include "TKey.h"
26 #include "TMD5.h"
27 #include "TObjString.h"
28 #include "TROOT.h"
29 #include "TString.h"
30 #include "TSystem.h"
31 
32 namespace mapkey {
33 // The "keys" below are indices for the vectors
34 //(of size mapkey::end)
35 // we use to look up the numbers we need
36 enum key
37 {
38  sf = 0,
39  stat = 1,
40  eig = 2,
41  uncorr = 3,
42  sys = 4,
43  end = 5
44 };
45 const char*
47 {
48  switch (input) {
49  case (sf):
50  return "sf";
51  case (stat):
52  return "stat";
53  case (eig):
54  return "eig";
55  case (uncorr):
56  return "uncorr";
57  case (sys):
58  return "sys";
59  }
60  return "";
61 }
62 }
63 
64 namespace {
65 char const * const LowPt_string = "LowPt" ;
66 const std::vector<int> s_keys = { mapkey::sf,
70 }
71 
73  const char* name)
74  : asg::AsgMessaging(std::string(name))
75  , m_doToyMC(false)
76  , m_doCombToyMC(false)
77  , m_nToyMC(0)
78  , m_seed(0)
79  , m_nSysMax(0)
80  , m_histList{ mapkey::end }
81  , m_fastHistList{ mapkey::end }
82  , m_Rndm()
83 {
84 }
85 
86 int
88 {
89  // use an int as a StatusCode
90  int sc(1);
91  // Check if files are present
92  if (m_corrFileNameList.empty()) {
93  ATH_MSG_ERROR(" No file added!");
94  return 0;
95  }
96  ATH_MSG_DEBUG("Initializing tool with " << m_corrFileNameList.size()
97  << " configuration file(s)");
98 
99  if (m_doToyMC && m_doCombToyMC) {
100  ATH_MSG_ERROR(" Both regular and combined toy MCs booked!"
101  << " Only use one!");
102  return 0;
103  }
104  /*
105  * initialize the random number generator if toyMC propagation booked
106  * Use the 1st 4 bytes of the CheckSum of the reccomendation file as seed
107  */
108  if (m_doToyMC || m_doCombToyMC) {
109  if (m_seed == 0) {
110  // Use the name of the correction for auto-setting of the seed based on
111  // the md5-sum of the file
112  const std::unique_ptr<char[]> fname(
113  gSystem->ExpandPathName(m_corrFileNameList[0].c_str()));
114  std::unique_ptr<TMD5> tmd = std::make_unique<TMD5>();
115  const char* tmd_as_string = tmd->FileChecksum(fname.get())->AsString();
116  m_seed = *(reinterpret_cast<const unsigned long int*>(tmd_as_string));
117  ATH_MSG_DEBUG("Seed (automatically) set to " << m_seed);
118  } else {
119  ATH_MSG_DEBUG("Seed set to " << m_seed);
120  }
121  m_Rndm = TRandom3(m_seed);
122  }
123  /*
124  * Load the needed histograms
125  */
126  if (0 == getHistograms()) {
127  ATH_MSG_ERROR(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
128  << "! Problem when calling getHistograms()");
129  return 0;
130  }
131  const unsigned int nRunNumbersFull = m_begRunNumberList.size();
132  const unsigned int nRunNumbersFast = m_begRunNumberListFastSim.size();
133 
134  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
135  << "Found " << nRunNumbersFast
136  << " run number ranges for fast sim with a total of "
137  << m_fastHistList[mapkey::sf].size()
138  << " scale factor histograms.");
139 
140  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
141  << "Found " << nRunNumbersFull
142  << " run number ranges for full sim with a total of "
143  << m_histList[mapkey::sf].size()
144  << " scale factor histograms.");
145 
146  ATH_MSG_DEBUG("Tool succesfully initialized!");
147 
148  return sc;
149 }
150 
152 {
154  const auto& currentmap = (isFastSim) ? m_fastHistList : m_histList;
155 
156  return currentmap.at(mapkey::stat).empty() && currentmap.at(mapkey::uncorr).empty();
157 }
158 
159 int
162  const unsigned int runnumber,
163  const double cluster_eta,
164  const double et, /* in MeV */
165  Result& result,
166  const bool onlyTotal) const
167 {
168  // Set up the non-0 defaults
169  result.SF = -999;
170  result.Total = 1;
171  /*
172  * Determine Simulation flavour and find the run period
173  */
175  int runnumberIndex = -1;
176  if (isFastSim) {
177  for (unsigned int i = 0; i < m_begRunNumberListFastSim.size(); ++i) {
178  if (m_begRunNumberListFastSim[i] <= runnumber &&
179  runnumber <= m_endRunNumberListFastSim[i]) {
180  runnumberIndex = i;
181  break;
182  }
183  }
184  } else {
185  for (unsigned int i = 0; i < m_begRunNumberList.size(); ++i) {
186  if (m_begRunNumberList[i] <= runnumber &&
187  runnumber <= m_endRunNumberList[i]) {
188  runnumberIndex = i;
189  break;
190  }
191  }
192  }
193  if (runnumberIndex < 0) {
194  return 0;
195  }
196  /* What we have is a vector<vector<HistArray>> acting as a map:
197  * Key/Index is SF,Stat,Eigen,UnCorr
198  * The Entry in this index is a vector<HistArray>
199  * Each vector<HistArray> has as many entries as supported Run periods.
200  * Each HistArray has 2D histos (could be standard, low-et, or forward
201  * electrons) The 2D Histo then has the number we want. What follows is the
202  * logic to get to this number.
203  */
204  const auto& currentmap = (isFastSim) ? m_fastHistList : m_histList;
205  const std::vector<HistArray>& sfVector = currentmap.at(mapkey::sf);
206  /*
207  * See if we can find a vector for key SF in the map
208  * and then if we can get the corresponding HisArray
209  * for the run period.
210  */
211  if (sfVector.empty() || runnumberIndex >= static_cast<int>(sfVector.size())) {
212  return 0;
213  }
214  /*
215  * At this stage we have found the relevant period
216  * So we need to locate the right histogram.
217  */
218  const HistArray& sfObjectArray = sfVector[runnumberIndex];
219  const auto& edges = (isFastSim) ? m_fastHistEdges[runnumberIndex]
220  : m_histEdges[runnumberIndex];
221  const int entries = edges.size();
222  /*
223  * Now the logic of finding the histogram
224  * Some parts of the code can be perhaps improved ...
225  */
226  double xValue(et);
227  double yValue(cluster_eta);
228  int smallEt(0);
229  int etaCov(0);
230  int nSF(0);
231  bool invalid = false;
232  bool changedEt = false;
233  int index = -1;
234 
235  for (int i = 0; i < entries; ++i) {
236  invalid = false;
237  const HistEdge& histEdge = edges[i];
238  // invalid if we are below minimum et
239  if (et < histEdge.etMin) {
240  smallEt++;
241  invalid = true;
242  }
243  // invalid if we are above max eta
244  if (std::abs(yValue) >= histEdge.etaMax) {
245  etaCov++;
246  invalid = true;
247  }
248  // invalid if we are less than minimum eta (forward electrons)
249  if (std::abs(yValue) < histEdge.etaMin) {
250  etaCov++;
251  invalid = true;
252  }
253  /*
254  * Invalid if above max et and is a low Et histogram.
255  * If not low Et histogram then change the xValue to the maximum
256  * availabe Et of ths histogram. As we assume that the SF stays the same
257  * for very high Et
258  */
259  if (et > histEdge.etMax) {
260  if (histEdge.isLowPt) {
261  invalid = true;
262  } else {
263  xValue = histEdge.etMax - 1000 ; //1000 is 1 GeV
264  changedEt = true;
265  }
266  }
267  /*
268  * Get the histogram index in the TObjArray
269  * Also mark how many times we found something
270  * as SF should be unique
271  */
272  if (!invalid) {
273  index = i;
274  if (!changedEt) {
275  nSF++;
276  }
277  }
278  }
279  if (smallEt == entries) {
280  return 0;
281  }
282  if (etaCov == entries) {
283  return 0;
284  }
285  if (nSF > 1) {
286  ATH_MSG_WARNING("More than 1 SF found for eta="
287  << yValue << " , et = " << et << " , run number = "
288  << runnumber << ". Please check your input files!");
289  }
290  if (index < 0) {
291  return 0;
292  }
293 
294  /*
295  * Now we have the index of the histogram
296  */
297  const HistEdge& currentEdge = edges[index];
298  /*
299  * If SF is only given in Abs(eta) convert eta input to std::abs()
300  */
301  constexpr double epsilon = 1e-6;
302  if (currentEdge.etaMin >= (0 - epsilon)) {
303  yValue = std::abs(yValue);
304  }
305 
306  const TH2* currentHist = static_cast<TH2*>(sfObjectArray[index].get()) ;
307  const int globalBinNumber = currentHist->FindFixBin(xValue, yValue);
308  const double scaleFactor = currentHist->GetBinContent(globalBinNumber);
309  const double scaleFactorErr = currentHist->GetBinError(globalBinNumber);
310  /*
311  * Write the retrieved values to the output
312  */
313  result.SF= scaleFactor;
314  result.Total = scaleFactorErr;
315  result.histIndex = index;
316  result.histBinNum = globalBinNumber;
317  /*
318  * if we only wanted the Total we can exit here
319  */
320  if (onlyTotal) {
321  return 1;
322  }
323  /*
324  * Do the stat error using the available info from the above (SF)
325  */
326  double statErr = -999;
327  const std::vector<HistArray>& statVector = currentmap.at(mapkey::stat);
328  if (!statVector.empty()) {
329  if (!statVector[runnumberIndex].empty()) {
330  statErr = static_cast<TH1*>(statVector[runnumberIndex][index].get())
331  ->GetBinContent(globalBinNumber);
332  result.Stat = statErr;
333  }
334  }
335  /*
336  * Do the Uncorr uncertainty
337  */
338  double val = statErr;
339  const std::vector<HistArray>& uncorrVector = currentmap.at(mapkey::uncorr);
340  if (!uncorrVector.empty()) {
341  if (!uncorrVector.at(runnumberIndex).empty()) {
342  const double valAdd =
343  static_cast<TH1*>(uncorrVector[runnumberIndex][index].get())
344  ->GetBinContent(globalBinNumber);
345  val = sqrt(val * val + valAdd * valAdd);
346  }
347  }
348  result.UnCorr = val;
349  /*
350  * Do the correlated part
351  * For the N~16 systematic variations
352  * we keep them in a vector of vector of HistArray
353  * The first vector index being the runnumber
354  * The second the systematic
355  * And them the HistArray for high low etc.
356  */
357  result.Corr.resize(m_nSysMax);
358  const std::vector<std::vector<HistArray>>& sysList =
359  (isFastSim) ? m_fastSysList : m_sysList;
360  if (sysList.size() > static_cast<unsigned int>(index)) {
361  if (sysList.at(index).size() > static_cast<unsigned int>(runnumberIndex)) {
362  const int sys_entries = sysList.at(index).at(runnumberIndex).size();
363  for (int sys = 0; sys < sys_entries; ++sys) {
364  double sysVal =
365  static_cast<TH2*>(sysList[index][runnumberIndex][sys].get())
366  ->GetBinContent(globalBinNumber);
367  result.Corr[sys] = sysVal;
368  }
369  }
370  }
371  /*
372  * Do the toys if requested
373  */
374  if (m_doToyMC || m_doCombToyMC) {
375  result.toys.resize(static_cast<size_t>(m_nToyMC));
376  const std::vector<std::vector<HistArray>>& toyMCList =
377  ((isFastSim) ? m_uncorrToyMCSystFast : m_uncorrToyMCSystFull);
378  if (toyMCList.size() > (unsigned int)runnumberIndex) {
379  for (int toy = 0; toy < m_nToyMC; ++toy) {
380  if (toyMCList[runnumberIndex][toy].size() >
381  static_cast<unsigned int>(index)) {
382  result.toys[toy] =
383  static_cast<TH2*>(toyMCList[runnumberIndex][toy][index].get())
384  ->GetBinContent(globalBinNumber);
385  }
386  }
387  }
388  }
389  return 1;
390 }
391 /*
392  * Build the toyMC tables from inputs
393  * Ownership should be tranfered to the map of the tables
394  * and the proper delete happens in the dtor
395  */
396 std::vector<TH2*>
398  const TH2* sf,
399  const TH2* stat,
400  const TH2* uncorr,
401  const std::vector<TH1*>& corr,
402  int& randomCounter)
403 {
404  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")! "
405  << "Entering function buildSingleToyMC");
406  std::vector<TH2*> tmpHists;
407  int nBins = (stat->GetNbinsX() + 2) * (stat->GetNbinsY() + 2);
408  tmpHists.reserve(m_nToyMC);
409  for (int toy = 0; toy < m_nToyMC; toy++) {
410  tmpHists.push_back((TH2*)corr.at(0)->Clone());
411  }
412  // Loop over all bins
413  for (int bin = 0; bin < nBins; bin++) {
414  double val = stat->GetBinContent(bin);
415 
416  // Add uncorrelated systematics
417  if (uncorr != nullptr) {
418  double valAdd = uncorr->GetBinContent(bin);
419  val = sqrt(val * val + valAdd * valAdd);
420  }
421  for (int toy = 0; toy < m_nToyMC; toy++) {
422  tmpHists.at(toy)->SetBinContent(
423  bin, (val * m_Rndm.Gaus(0, 1)) + sf->GetBinContent(bin));
424  randomCounter++;
425  tmpHists.at(toy)->SetDirectory(nullptr);
426  }
427  }
428  return tmpHists;
429 }
430 /*
431  * Build the combined toyMC tables from inputs
432  * Ownership should be tranfered to the unique_ptr
433  * in the lookup tables and the proper delete
434  * happens in the dtor
435  */
436 TH2*
438  const TH2* sf,
439  const TH2* stat,
440  const TH2* uncorr,
441  const std::vector<TH1*>& corr,
442  const int nSys,
443  int& randomCounter)
444 {
445 
446  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
447  << "Entering function buildSingleCombToyMC");
448 
449  TH2* tmpHist;
450  const int nBins = (stat->GetNbinsX() + 2) * (stat->GetNbinsY() + 2);
451  tmpHist = (TH2*)corr.at(0)->Clone();
452  // Create random numbers for the corr. uncertainties
453  std::vector<double> rnd(nSys, 0);
454  for (int s = 0; s < nSys; ++s) {
455  rnd[s] = m_Rndm.Gaus(0, 1);
456  randomCounter++;
457  }
458  // Loop over all bins
459  for (int bin = 0; bin < nBins; ++bin) {
460  double val = stat->GetBinContent(bin);
461 
462  // Add uncorrelated systematics
463  if (uncorr != nullptr) {
464  double valAdd = uncorr->GetBinContent(bin);
465  val = sqrt(val * val + valAdd * valAdd);
466  }
467  val = val * m_Rndm.Gaus(0, 1);
468  randomCounter++;
469  // Add larger correlated systematics
470  for (int s = 0; s < nSys; ++s) {
471  if (corr.at(s) != nullptr) {
472  val += static_cast<TH2*>(corr.at(s))->GetBinContent(bin) * rnd[s];
473  }
474  }
475  tmpHist->SetBinContent(bin, val + sf->GetBinContent(bin));
476  }
477  tmpHist->SetDirectory(nullptr);
478  return tmpHist;
479 }
480 /*
481  * Build the toyMC tables from inputs
482  */
483 std::vector<Root::TElectronEfficiencyCorrectionTool::HistArray>
485  const std::vector<TH1*>& sf,
486  const std::vector<TH1*>& eig,
487  const std::vector<TH1*>& stat,
488  const std::vector<TH1*>& uncorr,
489  const std::vector<std::vector<TH1*>>& corr)
490 {
491 
492  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
493  << "Entering function buildToyMCTable");
494 
495  int nSys{};
496  int randomCounter(0);
497  std::vector<HistArray> tmpVec;
498  const int stat_entries = stat.size();
499  if (m_doCombToyMC) {
500  for (int toyMC = 0; toyMC < m_nToyMC; toyMC++) {
501  HistArray tmpArray;
502  for (int i = 0; i < stat_entries; ++i) {
503  if (!eig.empty() && !uncorr.empty()) {
504  nSys = ((TH1*)eig.at(i))->GetNbinsX() - 1;
505  tmpArray.emplace_back(buildSingleCombToyMC((TH2*)sf.at(i),
506  (TH2*)stat.at(i),
507  (TH2*)uncorr.at(i),
508  corr.at(i),
509  nSys,
510  randomCounter));
511  } else {
512  tmpArray.emplace_back(buildSingleCombToyMC((TH2*)sf.at(i),
513  (TH2*)stat.at(i),
514  nullptr,
515  corr.at(i),
516  nSys,
517  randomCounter));
518  }
519  }
520  tmpVec.emplace_back(std::move(tmpArray));
521  }
522  } else {
523  std::vector<std::vector<TH2*>> tmpVec2;
524  for (int i = 0; i < stat_entries; ++i) {
525  nSys = ((TH1*)eig.at(i))->GetNbinsX() - 1;
526  tmpVec2.push_back(buildSingleToyMC((TH2*)sf.at(i),
527  (TH2*)stat.at(i),
528  (TH2*)uncorr.at(i),
529  corr.at(i),
530  randomCounter));
531  }
532  for (int toy = 0; toy < m_nToyMC; toy++) {
533  HistArray tmpArray;
534  for (auto& i : tmpVec2) {
535  tmpArray.emplace_back(i.at(toy));
536  }
537  tmpVec.emplace_back(std::move(tmpArray));
538  }
539  }
540  return tmpVec;
541 }
542 /*
543  * Helper function to retrieve number of uncorrelated bins
544  */
545 int
547  std::map<float, std::vector<float>>& pt_eta1) const
548 {
549  // Get sf histograms
550  const std::vector<HistArray>& tmpVec = (!m_histList[mapkey::sf].empty())
551  ? m_histList[mapkey::sf]
552  : m_fastHistList[mapkey::sf];
553 
554  int nbinsTotal = 0;
555  pt_eta1.clear();
556  std::vector<float> eta1;
557  eta1.clear();
558 
559  // Loop over the different Run range (one TObjeArray for each)
560  for (const auto& ikey : tmpVec) {
561  // Loop over the histograms for a given run numbers
562  for (const auto& entries : ikey) {
563  eta1.clear();
564  // Get number of bins
565  TH2* h_tmp = ((TH2*)entries.get());
566  int nbinsX = h_tmp->GetNbinsX();
567  int nbinsY = h_tmp->GetNbinsY();
568  // fill in the eta pushing back
569  for (int biny = 1; biny <= nbinsY; ++biny) {
570  eta1.push_back(h_tmp->GetYaxis()->GetBinLowEdge(biny));
571  }
572  // associate each pt (bin) with the corresponding/available eta ones
573  for (int binx = 1; binx <= nbinsX; ++binx) {
574  pt_eta1[h_tmp->GetXaxis()->GetBinLowEdge(binx)] = eta1;
575  }
576  }
577  }
578  for (auto& i : pt_eta1) {
579  nbinsTotal += i.second.size();
580  }
581  return nbinsTotal;
582 }
583 /*
584  * Get the histograms from the input files
585  */
586 int
588 {
589 
590  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
591  << "Entering function getHistograms");
592  // Cache the current directory in the root file
593  TDirectory* origDir = gDirectory;
594  /*
595  * Get the name of the first input ROOT file and
596  * interpret from that what we have:
597  * efficiency vs. efficiencySF; offline vs. trigger; medium, loose,...
598  */
599  if (!m_corrFileNameList.empty()) {
600  TString firstFileNameAndPath = m_corrFileNameList[0].c_str();
601  std::unique_ptr<TObjArray> myStringList(firstFileNameAndPath.Tokenize("/"));
602  int lastIdx = myStringList->GetLast();
603  TString fileName = ((TObjString*)myStringList->At(lastIdx))->GetString();
604  std::unique_ptr<TObjArray> myFileNameTokensList(fileName.Tokenize("."));
605 
606  if (myFileNameTokensList->GetLast() < 3) {
607  ATH_MSG_ERROR("input file name has wrong format!");
608  return 0;
609  }
610  }
611  /*
612  * Get all ROOT files and histograms
613  */
614  for (auto& ifile : m_corrFileNameList) {
615  // Load the ROOT file
616  const std::unique_ptr<char[]> fname(gSystem->ExpandPathName(ifile.c_str()));
617  std::unique_ptr<TFile> rootFile(TFile::Open(fname.get(), "READ"));
618  if (!rootFile) {
619  ATH_MSG_ERROR("No ROOT file found here: " << ifile);
620  return 0;
621  }
622  // Loop over all directories inside the root file (correspond to the run
623  // number ranges
624  TIter nextdir(rootFile->GetListOfKeys());
625  TKey* dir = nullptr;
626  TObject* obj = nullptr;
627  while ((dir = (TKey*)nextdir())) {
628  obj = dir->ReadObj();
629  if (obj->IsA()->InheritsFrom("TDirectory")) {
630  // splits string by delimiter --> e.g RunNumber1_RunNumber2
631  std::unique_ptr<TObjArray> dirNameArray(
632  TString(obj->GetName()).Tokenize("_"));
633  // returns index of last string --> if one, the directory name does not
634  // contain any run numbers
635  int lastIdx = dirNameArray->GetLast();
636  if (lastIdx != 1) {
638  "The folder name seems to have the wrong format! Directory name:"
639  << obj->GetName());
640  return 0;
641  }
642  rootFile->cd(obj->GetName());
643  if (0 == this->setupHistogramsInFolder(*dirNameArray, lastIdx)) {
644  ATH_MSG_ERROR("Unable to setup the histograms in directory "
645  << dir->GetName() << "in file " << ifile);
646  return 0;
647  }
648  } else {
649  ATH_MSG_ERROR("Wrong file content! Expected only Directories "
650  << gDirectory->cd());
651  return 0;
652  }
653  // Return to the original ROOT directory
655  } // End: directory loop
656  } // End: file loop
657  return 1;
658 }
659 /*
660  * Get the input histograms from a given folder/run number range
661  */
662 int
664  const TObjArray& dirNameArray,
665  int lastIdx)
666 {
667 
668  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
669  << "Entering funtion setupHistogramsInFolder");
670 
671  int runNumBegin(-1);
672  TString myBegRunNumString =
673  ((TObjString*)dirNameArray.At(lastIdx - 1))->GetString();
674  if (myBegRunNumString.IsDigit()) {
675  runNumBegin = myBegRunNumString.Atoi();
676  }
677  int runNumEnd(-1);
678  TString myEndRunNumString =
679  ((TObjString*)dirNameArray.At(lastIdx))->GetString();
680  if (myEndRunNumString.IsDigit()) {
681  runNumEnd = myEndRunNumString.Atoi();
682  }
683  if (runNumBegin < 0 || runNumEnd < 0 || runNumEnd < runNumBegin) {
684  ATH_MSG_ERROR(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
685  << "Could NOT interpret the run number range: "
686  << runNumBegin << " - " << runNumEnd);
687  return 0;
688  }
690  //--> e.g. "sf" index 0 ,
691  std::vector<std::vector<TH1*>> objsFull(mapkey::end);
692  std::vector<std::vector<TH1*>> objsFast(mapkey::end);
693  // Vector to hold the N~16 systematic variations
694  std::vector<std::vector<TH1*>> sysObjsFull;
695  std::vector<std::vector<TH1*>> sysObjsFast;
696  TIter nextkey(gDirectory->GetListOfKeys());
697  TKey* key = nullptr;
698  TObject* obj = nullptr;
699  int seenSystematics = 0;
700  // Loop of the keys
701  while ((key = (TKey*)nextkey())) {
702  obj = key->ReadObj();
703  if (obj->IsA()->InheritsFrom("TH1")) {
704  // The histogram containing the scale factors need to end with _sf and
705  // need to contain either the string "FullSim" or "AtlFast"!
706  if (std::strstr(obj->GetName(), "FullSim") != nullptr) {
707  setupTempMapsHelper(
708  static_cast<TH1*>(obj), objsFull, sysObjsFull, seenSystematics);
709  } else if (std::strstr(obj->GetName(), "AtlFast") != nullptr) {
710  setupTempMapsHelper(
711  static_cast<TH1*>(obj), objsFast, sysObjsFast, seenSystematics);
712  } else {
713  ATH_MSG_ERROR("Could NOT interpret if the histogram: "
714  << obj->GetName() << " is full or fast simulation!");
715  return 0;
716  }
717  }
718  }
719  ATH_MSG_DEBUG(" (file: " << __FILE__ << ", line: " << __LINE__ << ")\n"
720  << "Setting up histograms for Run range "
721  << runNumEnd);
722  // Copy from the temporaries to the actual member variables
723  // via the setup function
724  for (int key : s_keys) {
725  if (!objsFull.at(key).empty()) {
726  if (0 == setup(objsFull.at(key),
727  m_histList[key],
728  m_begRunNumberList,
729  m_endRunNumberList,
730  runNumBegin,
731  runNumEnd)) {
732  ATH_MSG_ERROR("! Could NOT setup histogram " << key
733  << " for full sim!");
734  return 0;
735  }
736  }
737  if (!objsFast.at(key).empty()) {
738  if (0 == setup(objsFast.at(key),
739  m_fastHistList[key],
740  m_begRunNumberListFastSim,
741  m_endRunNumberListFastSim,
742  runNumBegin,
743  runNumEnd)) {
744  ATH_MSG_ERROR("! Could NOT setup histogram " << key << " for fast sim");
745  return 0;
746  }
747  }
748  }
749  m_fastSysList.resize(sysObjsFast.size());
750  for (unsigned int sys = 0; sys < sysObjsFast.size(); sys++) {
751  if (0 == setup(sysObjsFast.at(sys),
752  m_fastSysList[sys],
753  m_begRunNumberListFastSim,
754  m_endRunNumberListFastSim,
755  runNumBegin,
756  runNumEnd)) {
757  ATH_MSG_ERROR("! Could NOT setup systematic histograms for fast sim");
758  return 0;
759  }
760  }
761  m_sysList.resize(sysObjsFull.size());
762  for (unsigned int sys = 0; sys < sysObjsFull.size(); sys++) {
763  if (0 == setup(sysObjsFull.at(sys),
764  m_sysList[sys],
765  m_begRunNumberList,
766  m_endRunNumberList,
767  runNumBegin,
768  runNumEnd)) {
769  ATH_MSG_ERROR("! Could NOT setup systematic histograms for fast sim");
770  return 0;
771  }
772  }
773  //Histogram edges
774  fillHistEdges(m_histList.at(mapkey::sf), m_histEdges);
775  fillHistEdges(m_fastHistList.at(mapkey::sf), m_fastHistEdges);
776 
777  // Toys
778  if (m_doToyMC || m_doCombToyMC) {
779  bool fullToysBooked =
780  setupUncorrToySyst(objsFull, sysObjsFull, m_uncorrToyMCSystFull);
781  bool fastToysBooked =
782  setupUncorrToySyst(objsFast, sysObjsFast, m_uncorrToyMCSystFast);
783  if (fullToysBooked || fastToysBooked) {
784  if (m_doToyMC) {
785  ATH_MSG_DEBUG("Created tables for " << m_nToyMC
786  << " ToyMC systematics ");
787  }
788  if (m_doCombToyMC) {
789  ATH_MSG_DEBUG("Created tables for " << m_nToyMC
790  << " combined ToyMC systematics ");
791  }
792  }
793  }
794  return 1;
795 }
796 /*
797  * Helper for Setting up the temporary/intermediate maps
798  * from the histos
799  */
800 void
802  TH1* obj,
803  std::vector<std::vector<TH1*>>& objs,
804  std::vector<std::vector<TH1*>>& sysObjs,
805  int& seenSystematics)
806 {
807  // Add all except the correlated
808  for (int key : s_keys) {
809  if (TString(obj->GetName())
810  .EndsWith("_" + TString(mapkey::keytostring(key)))) {
811  objs.at(key).emplace_back(obj);
812  }
813  }
814 
815  const TString tmpName(obj->GetName());
816  // Special treatment , this is only for photons
817  if (tmpName.EndsWith("_sys")) {
818  objs.at(mapkey::sys).emplace_back(obj);
819  std::vector<TH1*> tmpArray;
820  // clone
821  tmpArray.emplace_back(static_cast<TH1*>(obj->Clone()));
822  sysObjs.emplace_back(tmpArray);
823  seenSystematics++;
824  }
825  // See if we are dealing with correlated
826  if (tmpName.Contains("_corr")) {
827  /*
828  * corr0 in the name triggers a few things
829  * We assume that 0 is the 1st
830  * histogram in a series of corr(i) that
831  * we see for each of the vector entries that
832  * can be one for LowPt,Standard,Forward etc
833  */
834  if (tmpName.EndsWith("corr0")) {
835  // 1st create a TObjectArray
836  std::vector<TH1*> tmpArray;
837  // Register it to the vector
838  sysObjs.emplace_back(tmpArray);
839  // Reset the counter here
840  seenSystematics = 0;
841  }
842  /*
843  * Now we can add to the TObjeArray
844  * This can be Low Pt or high Pt
845  */
846  sysObjs.back().emplace_back(obj);
847  //Increase the counter
848  seenSystematics++;
849  }
850 
851  if (seenSystematics > m_nSysMax) {
852  m_nSysMax = seenSystematics;
853  }
854 }
855 /*
856  * Helper for Setting up the uncorrelated syst for the toys
857  */
858 bool
860  std::vector<std::vector<TH1*>>& objs,
861  std::vector<std::vector<TH1*>>& sysObjs,
862  std::vector<std::vector<HistArray>>& uncorrToyMCSyst)
863 {
864  bool toysBooked = false;
865  if (!m_histList[mapkey::sf].empty()) {
866  if (objs.at(mapkey::eig).empty() || objs.at(mapkey::stat).empty() ||
867  objs.at(mapkey::uncorr).empty()) {
868 
869  if (objs.at(mapkey::stat).size() > 1 || objs.at(mapkey::sys).size() > 1) {
870 
871  uncorrToyMCSyst.push_back(buildToyMCTable(
872  objs.at(mapkey::sf), {}, objs.at(mapkey::stat), {}, sysObjs));
873  toysBooked = true;
874  }
875  } else {
876  uncorrToyMCSyst.push_back(buildToyMCTable(objs.at(mapkey::sf),
877  objs.at(mapkey::eig),
878  objs.at(mapkey::stat),
879  objs.at(mapkey::uncorr),
880  sysObjs));
881  toysBooked = true;
882  }
883  }
884  return toysBooked;
885 }
886 /*
887  * Fill and interpret the setup, depending
888  * on which histograms are found in the input file(s)
889  */
890 int
892  const std::vector<TH1*>& hists,
893  std::vector<HistArray>& histList,
894  std::vector<unsigned int>& beginRunNumberList,
895  std::vector<unsigned int>& endRunNumberList,
896  const int runNumBegin,
897  const int runNumEnd) const
898 {
899  if (hists.empty()) {
900  ATH_MSG_ERROR("! Could NOT find histogram with name *_sf in folder");
901  return 0;
902  }
903  TH1* tmpHist(nullptr);
904  HistArray tmpArray;
905  for (const auto& hist : hists) {
906  tmpHist = static_cast<TH1*>(hist);
907  tmpHist->SetDirectory(nullptr);
908  tmpArray.emplace_back(tmpHist);
909  }
910  histList.emplace_back(std::move(tmpArray));
911  // Now, we have all the needed info. Fill the vectors accordingly
912  if (!beginRunNumberList.empty()) {
913  if (runNumBegin != (int)beginRunNumberList.back()) {
914  beginRunNumberList.push_back(runNumBegin);
915  }
916  } else {
917  beginRunNumberList.push_back(runNumBegin);
918  }
919  if (!endRunNumberList.empty()) {
920  if (runNumEnd != (int)endRunNumberList.back()) {
921  endRunNumberList.push_back(runNumEnd);
922  }
923  } else {
924  endRunNumberList.push_back(runNumEnd);
925  }
926  return 1;
927 }
928 
929 
931  const std::vector<HistArray>& sfPerPeriodHist,
932  std::vector<std::vector<HistEdge>>& sfPerPeriodEdges) {
933 
934  for (const auto& vec : sfPerPeriodHist) {
935  const size_t vecSize = vec.size();
936  std::vector<HistEdge> periodVec;
937  periodVec.reserve(vecSize);
938  for (size_t i = 0; i < vecSize; ++i) {
939  const auto* tmpHist = static_cast<TH2*>(vec[i].get());
940  const auto* const xAxis = tmpHist->GetXaxis();
941  const auto* yAxis = tmpHist->GetYaxis();
942  HistEdge histEdge;
943  histEdge.etaMax = yAxis->GetXmax();
944  histEdge.etaMin = yAxis->GetXmin();
945  histEdge.etMax = xAxis->GetXmax();
946  histEdge.etMin = xAxis->GetXmin();
947  histEdge.isLowPt =
948  (std::strstr(tmpHist->GetName(), LowPt_string) != nullptr);
949 
950  periodVec.emplace_back(histEdge);
951  }
952  sfPerPeriodEdges.emplace_back(std::move(periodVec));
953  }
954 }
955 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Root::TElectronEfficiencyCorrectionTool::HistEdge::etaMin
double etaMin
Definition: TElectronEfficiencyCorrectionTool.h:162
et
Extra patterns decribing particle interation process.
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
get_generator_info.result
result
Definition: get_generator_info.py:21
mapkey::keytostring
const char * keytostring(int input)
Definition: TElectronEfficiencyCorrectionTool.cxx:46
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Root::TElectronEfficiencyCorrectionTool::setupUncorrToySyst
bool setupUncorrToySyst(std::vector< std::vector< TH1 * >> &objs, std::vector< std::vector< TH1 * >> &sysObjs, std::vector< std::vector< HistArray >> &uncorrToyMCSyst)
Definition: TElectronEfficiencyCorrectionTool.cxx:859
mapkey::uncorr
@ uncorr
Definition: TElectronEfficiencyCorrectionTool.cxx:41
index
Definition: index.py:1
plotmaker.hist
hist
Definition: plotmaker.py:148
PATCore::ParticleDataType::Fast
@ Fast
Definition: PATCoreEnums.h:22
asg
Definition: DataHandleTestTool.h:28
bin
Definition: BinsDiffFromStripMedian.h:43
downloadSingle.dataType
string dataType
Definition: downloadSingle.py:18
xAOD::eta1
setEt setPhi setE277 setWeta2 eta1
Definition: TrigEMCluster_v1.cxx:41
Root::TElectronEfficiencyCorrectionTool::uncorrEmpty
bool uncorrEmpty(const PATCore::ParticleDataType::DataType dataType)
Check if stat+uncorr has input.
Definition: TElectronEfficiencyCorrectionTool.cxx:151
vec
std::vector< size_t > vec
Definition: CombinationsGeneratorTest.cxx:12
empty
bool empty(TH1 *h)
Definition: computils.cxx:295
mapkey::sys
@ sys
Definition: TElectronEfficiencyCorrectionTool.cxx:42
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
PATCore::ParticleDataType::DataType
DataType
Definition: PATCoreEnums.h:22
Root::TElectronEfficiencyCorrectionTool::HistEdge::isLowPt
bool isLowPt
Definition: TElectronEfficiencyCorrectionTool.h:165
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
Root::TElectronEfficiencyCorrectionTool::Result
Definition: TElectronEfficiencyCorrectionTool.h:42
fitman.tmpHist
def tmpHist(what, wmin=-1e10, wmax=+1e10)
Definition: fitman.py:146
FortranAlgorithmOptions.fileName
fileName
Definition: FortranAlgorithmOptions.py:13
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
lumiFormat.i
int i
Definition: lumiFormat.py:85
Root::TElectronEfficiencyCorrectionTool::HistEdge::etaMax
double etaMax
Definition: TElectronEfficiencyCorrectionTool.h:161
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
mapkey
Definition: TElectronEfficiencyCorrectionTool.cxx:32
mergePhysValFiles.origDir
origDir
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:24
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
checkFileSG.objs
list objs
Definition: checkFileSG.py:93
Root::TElectronEfficiencyCorrectionTool::buildSingleToyMC
std::vector< TH2 * > buildSingleToyMC(const TH2 *sf, const TH2 *stat, const TH2 *uncorr, const std::vector< TH1 * > &corr, int &randomCounter)
Definition: TElectronEfficiencyCorrectionTool.cxx:397
Root::TElectronEfficiencyCorrectionTool::getHistograms
int getHistograms()
Load all histograms from the input file(s)
Definition: TElectronEfficiencyCorrectionTool.cxx:587
mapkey::end
@ end
Definition: TElectronEfficiencyCorrectionTool.cxx:43
MakeTH3DFromTH2Ds.hists
hists
Definition: MakeTH3DFromTH2Ds.py:72
makeComparison.rootFile
rootFile
Definition: makeComparison.py:27
jobOptions.isFastSim
string isFastSim
Definition: jobOptions.py:54
mapkey::eig
@ eig
Definition: TElectronEfficiencyCorrectionTool.cxx:40
DeMoScan.runnumber
runnumber
Definition: DeMoScan.py:266
beamspotman.dir
string dir
Definition: beamspotman.py:623
Root::TElectronEfficiencyCorrectionTool::initialize
int initialize()
Initialize this class.
Definition: TElectronEfficiencyCorrectionTool.cxx:87
python.InDetPriVxFinderConfig.setup
setup
Definition: InDetPriVxFinderConfig.py:124
dumpTgcDigiJitter.nBins
list nBins
Definition: dumpTgcDigiJitter.py:29
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:225
Root::TElectronEfficiencyCorrectionTool::setup
int setup(const std::vector< TH1 * > &hists, std::vector< HistArray > &histList, std::vector< unsigned int > &beginRunNumberList, std::vector< unsigned int > &endRunNumberList, const int runNumBegin, const int runNumEnd) const
Fill and interpret the setup, depending on which histograms are found in the input file(s)
Definition: TElectronEfficiencyCorrectionTool.cxx:891
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
Root::TElectronEfficiencyCorrectionTool::buildToyMCTable
std::vector< HistArray > buildToyMCTable(const std::vector< TH1 * > &sf, const std::vector< TH1 * > &eig, const std::vector< TH1 * > &stat, const std::vector< TH1 * > &uncorr, const std::vector< std::vector< TH1 * >> &corr)
Definition: TElectronEfficiencyCorrectionTool.cxx:484
Root::TElectronEfficiencyCorrectionTool::HistArray
std::vector< std::unique_ptr< TH1 > > HistArray
Definition: TElectronEfficiencyCorrectionTool.h:39
Root::TElectronEfficiencyCorrectionTool::fillHistEdges
static void fillHistEdges(const std::vector< HistArray > &sfPerPeriodHist, std::vector< std::vector< HistEdge >> &sfPerPeriodEdges)
Definition: TElectronEfficiencyCorrectionTool.cxx:930
mapkey::stat
@ stat
Definition: TElectronEfficiencyCorrectionTool.cxx:39
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:67
Root::TElectronEfficiencyCorrectionTool::HistEdge::etMin
double etMin
Definition: TElectronEfficiencyCorrectionTool.h:164
DeMoScan.index
string index
Definition: DeMoScan.py:364
Root::TElectronEfficiencyCorrectionTool::setupTempMapsHelper
void setupTempMapsHelper(TH1 *obj, std::vector< std::vector< TH1 * >> &objs, std::vector< std::vector< TH1 * >> &sysObjs, int &seenSystematics)
Definition: TElectronEfficiencyCorrectionTool.cxx:801
Root::TElectronEfficiencyCorrectionTool::TElectronEfficiencyCorrectionTool
TElectronEfficiencyCorrectionTool(const char *name="TElectronEfficiencyCorrectionTool")
Standard constructor.
Definition: TElectronEfficiencyCorrectionTool.cxx:72
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Root::TElectronEfficiencyCorrectionTool::HistEdge
Definition: TElectronEfficiencyCorrectionTool.h:160
Root::TElectronEfficiencyCorrectionTool::getNbins
int getNbins(std::map< float, std::vector< float >> &ptEta) const
Helpers to get the binning of the uncertainties in a std::map (pt, eta)
Definition: TElectronEfficiencyCorrectionTool.cxx:546
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
entries
double entries
Definition: listroot.cxx:49
Root::TElectronEfficiencyCorrectionTool::HistEdge::etMax
double etMax
Definition: TElectronEfficiencyCorrectionTool.h:163
LArCellNtuple.ifile
string ifile
Definition: LArCellNtuple.py:133
TElectronEfficiencyCorrectionTool.h
Root::TElectronEfficiencyCorrectionTool::setupHistogramsInFolder
int setupHistogramsInFolder(const TObjArray &dirNameArray, int lastIdx)
Definition: TElectronEfficiencyCorrectionTool.cxx:663
python.PyAthena.obj
obj
Definition: PyAthena.py:132
Root::TElectronEfficiencyCorrectionTool::calculate
int calculate(const PATCore::ParticleDataType::DataType dataType, const unsigned int runnumber, const double cluster_eta, const double et, Result &result, const bool onlyTotal=false) const
The main calculate method: dataType PATCore::ParticleDataType::DataType (e.g DATA,...
Definition: TElectronEfficiencyCorrectionTool.cxx:160
Root::TElectronEfficiencyCorrectionTool::buildSingleCombToyMC
TH2 * buildSingleCombToyMC(const TH2 *sf, const TH2 *stat, const TH2 *uncorr, const std::vector< TH1 * > &corr, const int nSys, int &randomCounter)
Definition: TElectronEfficiencyCorrectionTool.cxx:437
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37