ATLAS Offline Software
Loading...
Searching...
No Matches
Functions
compareFlatTrees.cxx File Reference
#include <ROOT/RDataFrame.hxx>
#include <ROOT/RLogger.hxx>
#include <TCanvas.h>
#include <TF1.h>
#include <TH1D.h>
#include <TRatioPlot.h>
#include <TROOT.h>
#include <TStyle.h>
#include <CxxUtils/checker_macros.h>
#include <boost/program_options.hpp>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <memory>
#include <string>
#include <unordered_set>

Go to the source code of this file.

Functions

std::vector< std::string > intersection (std::vector< std::string > &v1, std::vector< std::string > &v2)
std::vector< std::string > remainder (const std::vector< std::string > &v1, const std::vector< std::string > &v2)
int main ATLAS_NOT_THREAD_SAFE (int argc, char *argv[])

Function Documentation

◆ ATLAS_NOT_THREAD_SAFE()

int main ATLAS_NOT_THREAD_SAFE ( int argc,
char * argv[] )

Definition at line 58 of file compareFlatTrees.cxx.

59{
60 namespace po = boost::program_options;
61 po::options_description poDescription("Common options");
62 poDescription.add_options()
63 ("help", "produce help message")
64 ("require-same-branches", "require both trees to have the same branches")
65 ("tree-name", po::value<std::string>()->required(), "tree name")
66 ("reference-file", po::value<std::string>()->required(), "reference file(s), wildcards supported")
67 ("test-file", po::value<std::string>()->required(), "test file(s), wildcards supported")
68 ("branch-name", po::value<std::string>(), "base branch name (optional)");
69
70 po::options_description poDescriptionAdvanced("Advanced options");
71 poDescriptionAdvanced.add_options()
72 ("scale", "scale histograms that both have the same event count")
73 ("rebin", "do smart rebinning")
74 ("benchmark", "benchmark the code")
75 ("verbose", "verbose logging");
76
77 po::options_description poDescriptionAll;
78 poDescriptionAll.add(poDescription).add(poDescriptionAdvanced);
79
80 po::positional_options_description poPositionalOptions;
81 poPositionalOptions.add("tree-name", 1);
82 poPositionalOptions.add("reference-file", 1);
83 poPositionalOptions.add("test-file", 1);
84 poPositionalOptions.add("branch-name", 1);
85
86 po::variables_map poVariablesMap;
87 po::store(po::command_line_parser(argc, argv)
88 .options(poDescriptionAll)
89 .positional(poPositionalOptions).run(),
90 poVariablesMap);
91
92 if (poVariablesMap.count("help"))
93 {
94 std::cout << "Usage: compareFlatTrees [OPTION] tree-name reference-file test-file [branch-name]" << std::endl;
95 std::cout << poDescriptionAll << std::endl;
96 return 0;
97 }
98
99 po::notify(poVariablesMap);
100
101 // Base name of branches to read
102 std::string treeName = poVariablesMap["tree-name"].as<std::string>();
103 std::string treeNameOut = treeName;
104 std::replace( treeNameOut.begin(), treeNameOut.end(), '/', '_');
105 std::string referenceInput = poVariablesMap["reference-file"].as<std::string>();
106 std::string testInput = poVariablesMap["test-file"].as<std::string>();
107 std::string baseBranchName;
108 std::string outputPDF;
109 if (poVariablesMap.count("branch-name") > 0)
110 {
111 baseBranchName = poVariablesMap["branch-name"].as<std::string>();
112 outputPDF = "comparison_" + treeNameOut + "_" + baseBranchName + ".pdf";
113 }
114 else
115 {
116 outputPDF = "comparison_" + treeNameOut + ".pdf";
117 }
118
119 bool scale = poVariablesMap.count("scale") > 0;
120 bool rebin = poVariablesMap.count("rebin") > 0;
121 bool benchmark = poVariablesMap.count("benchmark") > 0;
122 bool verbose = poVariablesMap.count("verbose") > 0;
123
124 // Verbose logging
125 // auto verbosity = ROOT::Experimental::RLogScopedVerbosity(ROOT::Detail::RDF::RDFLogChannel(), ROOT::Experimental::ELogLevel::kInfo);
126
127 // Run in batch mode - output is pdf
128 gROOT->SetBatch(kTRUE);
129 // Suppress logging
130 if (!verbose)
131 {
132 gErrorIgnoreLevel = kWarning;
133 }
134 // No stats box
135 gStyle->SetOptStat(0);
136 // No error bar ticks
137 gStyle->SetEndErrorSize(0);
138 // Parallel processing where possible
139 ROOT::EnableImplicitMT(4);
140
141 // Create RDataFrame
142 ROOT::RDataFrame dataFrameRefr(treeName, referenceInput);
143 ROOT::RDataFrame dataFrameTest(treeName, testInput);
144
145 // Event count and ratio
146 auto eventsRefr{dataFrameRefr.Count()};
147 auto eventsTest{dataFrameTest.Count()};
148 double eventRatio{static_cast<float>(eventsRefr.GetValue()) / static_cast<float>(eventsTest.GetValue())};
149
150 // Get column names for each file and then the intersection
151 auto colNamesRefr = dataFrameRefr.GetColumnNames();
152 auto colNamesTest = dataFrameTest.GetColumnNames();
153 auto colNames = intersection(colNamesRefr, colNamesTest);
154
155 bool checkMissColumns = false;
156 std::vector<std::string> missColNamesRefr, missColNamesTest;
157 if (poVariablesMap.count("require-same-branches"))
158 {
159 checkMissColumns = true;
160 missColNamesRefr = remainder (colNamesRefr, colNames);
161 missColNamesTest = remainder (colNamesTest, colNames);
162 }
163
164 // Loop over column names and get a list of the required columns
165 std::vector<std::string> requiredColumns;
166 std::cout << "Will attempt to plot the following columns:" << std::endl;
167 for (auto &&colName : colNames)
168 {
169 if ((baseBranchName.empty() || colName.find(baseBranchName) != std::string::npos) && // include
170 (colName.find("Link") == std::string::npos) && // exclude, elementlinks
171 (colName.find("m_persIndex") == std::string::npos) && // exclude, elementlinks
172 (colName.find("m_persKey") == std::string::npos) && // exclude, elementlinks
173 (colName.find("Parent") == std::string::npos) && // exclude, elementlinks
174 (colName.find("original") == std::string::npos) && // exclude, elementlinks
175 (colName.find("EventInfoAuxDyn.detDescrTags") == std::string::npos) && // exclude, std::pair
176 (dataFrameRefr.GetColumnType(colName).find("xAOD") == std::string::npos) && // exclude, needs ATLAS s/w
177 (dataFrameRefr.GetColumnType(colName) != "ROOT::VecOps::RVec<string>") && // exclude, needs ATLAS s/w
178 (!dataFrameRefr.GetColumnType(colName).starts_with("ROOT::VecOps::RVec<pair")) && // exclude, std::pair
179 (dataFrameRefr.GetColumnType(colName).find("vector") == std::string::npos))
180 { // exclude, needs unwrapping
181 requiredColumns.push_back(colName);
182 std::cout << " " << colName << " " << dataFrameRefr.GetColumnType(colName) << std::endl;
183 }
184 }
185
186 // Set binning
187 const int nBins{128};
188
189 // Loop over the required columns and plot them for each sample along with the ratio
190 // Write resulting plots to a pdf file
191 bool fileOpen{};
192 size_t counter{};
193 size_t failedCount{};
194 std::chrono::seconds totalDuration{};
195 std::unordered_map<std::string, ROOT::RDF::RResultPtr<double>> mapMinValuesRefr;
196 std::unordered_map<std::string, ROOT::RDF::RResultPtr<double>> mapMinValuesTest;
197 std::unordered_map<std::string, ROOT::RDF::RResultPtr<double>> mapMaxValuesRefr;
198 std::unordered_map<std::string, ROOT::RDF::RResultPtr<double>> mapMaxValuesTest;
199 std::unordered_map<std::string, ROOT::RDF::RResultPtr<TH1D>> mapHistRefr;
200 std::unordered_map<std::string, ROOT::RDF::RResultPtr<TH1D>> mapHistTest;
201
202 std::cout << "Preparing ranges..." << std::endl;
203 for (const std::string &colName : requiredColumns)
204 {
205 mapMinValuesRefr.emplace(colName, dataFrameRefr.Min(colName));
206 mapMinValuesTest.emplace(colName, dataFrameTest.Min(colName));
207 mapMaxValuesRefr.emplace(colName, dataFrameRefr.Max(colName));
208 mapMaxValuesTest.emplace(colName, dataFrameTest.Max(colName));
209 }
210
211 std::cout << "Preparing histograms..." << std::endl;
212 auto start = std::chrono::high_resolution_clock::now();
213 for (auto it = requiredColumns.begin(); it != requiredColumns.end();)
214 {
215 const std::string &colName = *it;
216 const char* colNameCh = colName.c_str();
217
218 // Initial histogram range
219 int nBinsForColumn = nBins;
220 double min = std::min(mapMinValuesRefr[colName].GetValue(), mapMinValuesTest[colName].GetValue());
221 double max = std::max(mapMaxValuesRefr[colName].GetValue(), mapMaxValuesTest[colName].GetValue());
222 if (min == max) {
223 min -= 0.5;
224 max += 0.5;
225 nBinsForColumn = 1;
226 } else {
227 max *= 1.02;
228 }
229 if (std::isinf(min) || std::isinf(max))
230 {
231 std::cout << " skipping " << colName << " ..." << std::endl;
232 it = requiredColumns.erase(it);
233 continue;
234 } else {
235 ++it;
236 }
237
238 if (max > 250e3 && min > 0.0)
239 {
240 min = 0.0;
241 }
242
243 if (verbose)
244 {
245 std::cout << " " << colName << " type: " << dataFrameRefr.GetColumnType(colName) << " min: " << min << " max: " << max << " nbins: " << nBinsForColumn << std::endl;
246 }
247
248 // Initial histograms
249 mapHistRefr.emplace(colName, dataFrameRefr.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh));
250 mapHistTest.emplace(colName, dataFrameTest.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh));
251 }
252 auto stop = std::chrono::high_resolution_clock::now();
253 auto duration = std::chrono::duration_cast<std::chrono::seconds>(stop - start);
254 totalDuration += duration;
255 if (benchmark)
256 {
257 std::cout << " Time for this step: " << duration.count() << " seconds " << std::endl;
258 std::cout << " Elapsed time: " << totalDuration.count() << " seconds (" << std::chrono::duration_cast<std::chrono::minutes>(totalDuration).count() << " minutes)" << std::endl;
259 }
260
261 if (rebin)
262 {
263 std::cout << "Rebinning histograms..." << std::endl;
264 auto start = std::chrono::high_resolution_clock::now();
265 for (const std::string &colName : requiredColumns)
266 {
267 const char* colNameCh = colName.c_str();
268 auto &histRefr = mapHistRefr[colName];
269 auto &histTest = mapHistTest[colName];
270
271 // Initial histogram range
272 int nBinsForColumn = nBins;
273 double min = std::min(mapMinValuesRefr[colName].GetValue(), mapMinValuesTest[colName].GetValue());
274 double max = std::max(mapMaxValuesRefr[colName].GetValue(), mapMaxValuesTest[colName].GetValue());
275 if (min == max) {
276 min -= 0.5;
277 max += 0.5;
278 nBinsForColumn = 1;
279 } else {
280 max *= 1.02;
281 }
282 if (max > 250e3 && min > 0.0)
283 {
284 min = 0.0;
285 }
286
287 // Check range - make sure that bins other than the first contain at least one per mille events
288 // Avoids case where max is determined by a single outlier leading to most events being in the 1st bin
289 bool rangeSatisfactory{};
290 size_t rangeItrCntr{};
291 while (!rangeSatisfactory && rangeItrCntr < 10)
292 {
293 ++rangeItrCntr;
294 if (verbose)
295 {
296 std::cout << std::endl
297 << " Range tuning... iteration number " << rangeItrCntr << std::endl;
298 }
299 double entriesFirstBin = histRefr.GetPtr()->GetBinContent(1);
300 double entriesLastBin = histRefr.GetPtr()->GetBinContent(nBinsForColumn);
301 double entriesOtherBins{};
302 for (size_t i{2}; i < static_cast<size_t>(nBinsForColumn); ++i)
303 {
304 entriesOtherBins += histRefr.GetPtr()->GetBinContent(i);
305 }
306 bool firstBinOK{((entriesOtherBins + entriesLastBin) / entriesFirstBin > 0.001f)};
307 bool lastBinOK{((entriesOtherBins + entriesFirstBin) / entriesLastBin > 0.001f)};
308 rangeSatisfactory = ((firstBinOK && lastBinOK) || entriesOtherBins == 0.0f);
309 if (!rangeSatisfactory)
310 {
311 if (verbose)
312 {
313 std::cout << "Min " << min << std::endl;
314 std::cout << "Max " << max << std::endl;
315 std::cout << "1st " << entriesFirstBin << std::endl;
316 std::cout << "Mid " << entriesOtherBins << std::endl;
317 std::cout << "End " << entriesLastBin << std::endl;
318 std::cout << "R/F " << (entriesOtherBins + entriesLastBin) / entriesFirstBin << std::endl;
319 std::cout << "R/L " << (entriesOtherBins + entriesFirstBin) / entriesLastBin << std::endl;
320 }
321 if (!firstBinOK)
322 {
323 max = (max - min) / static_cast<double>(nBinsForColumn);
324 if (verbose)
325 {
326 std::cout << " " << colName << " min: " << min << " max: " << max << " nbins: " << nBinsForColumn << std::endl;
327 }
328 histRefr = dataFrameRefr.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh);
329 histTest = dataFrameTest.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh);
330 }
331 if (!lastBinOK)
332 {
333 min = max * (1.0f - (1.0f / static_cast<double>(nBinsForColumn)));
334 if (verbose)
335 {
336 std::cout << " " << colName << " min: " << min << " max: " << max << " nbins: " << nBinsForColumn << std::endl;
337 }
338 histRefr = dataFrameRefr.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh);
339 histTest = dataFrameTest.Histo1D({colNameCh, colNameCh, nBinsForColumn, min, max}, colNameCh);
340 }
341 }
342 }
343 }
344
345 auto stop = std::chrono::high_resolution_clock::now();
346 auto duration = std::chrono::duration_cast<std::chrono::seconds>(stop - start);
347 totalDuration += duration;
348 if (benchmark)
349 {
350 std::cout << " Time for this step: " << duration.count() << " seconds " << std::endl;
351 std::cout << " Elapsed time: " << totalDuration.count() << " seconds (" << std::chrono::duration_cast<std::chrono::minutes>(totalDuration).count() << " minutes)" << std::endl;
352 }
353 }
354
355 std::cout << "Running comparisons..." << std::endl;
356 start = std::chrono::high_resolution_clock::now();
357
358 // Store only last canvas
359 std::unique_ptr<TCanvas> lastCanvas;
360
361 for (const std::string &colName : requiredColumns)
362 {
363 ++counter;
364
365 std::cout << "Processing column " << counter << " of " << requiredColumns.size() << " : " << colName << " ... ";
366
367 auto h1 = mapHistRefr[colName].GetPtr();
368 auto h2 = mapHistTest[colName].GetPtr();
369
370 if (scale)
371 {
372 h2->Scale(eventRatio);
373 }
374 h2->SetMarkerStyle(20);
375 h2->SetMarkerSize(0.8);
376
377 if (!verbose)
378 {
379 gErrorIgnoreLevel = kError; // this is spammy due to empty bins
380 }
381 auto c1 = std::make_unique<TCanvas>();
382 auto rp = std::make_unique<TRatioPlot>(h2, h1);
383 if (!verbose)
384 {
385 gErrorIgnoreLevel = kWarning;
386 }
387
388 rp->SetH1DrawOpt("PE");
389 rp->SetH2DrawOpt("hist");
390 rp->SetGraphDrawOpt("PE");
391 rp->Draw();
392 rp->GetUpperRefXaxis()->SetTitle(colName.c_str());
393 rp->GetUpperRefYaxis()->SetTitle("Count");
394 rp->GetLowerRefYaxis()->SetTitle("Test / Ref.");
395 rp->GetLowerRefGraph()->SetMarkerStyle(20);
396 rp->GetLowerRefGraph()->SetMarkerSize(0.8);
397
398 bool valid{true};
399 for (int i{}; i < rp->GetLowerRefGraph()->GetN(); i++)
400 {
401 if (rp->GetLowerRefGraph()->GetY()[i] != 1.0)
402 {
403 valid = false;
404 break;
405 }
406 }
407 if (valid)
408 {
409 std::cout << "PASS" << std::endl;
410 continue;
411 }
412 else
413 {
414 std::cout << "FAILED" << std::endl;
415 ++failedCount;
416 }
417
418 rp->GetLowerRefGraph()->SetMinimum(0.5);
419 rp->GetLowerRefGraph()->SetMaximum(1.5);
420 rp->GetLowYaxis()->SetNdivisions(505);
421
422 if (!fileOpen)
423 {
424 // Open file
425 c1->Print((outputPDF + "[").c_str());
426 fileOpen = true;
427 }
428 // Actual plot
429 c1->Print(outputPDF.c_str());
430 c1->Clear();
431 lastCanvas = std::move(c1);
432 }
433
434 if (fileOpen)
435 {
436 // Close file
437 lastCanvas->Print((outputPDF + "]").c_str());
438 lastCanvas.reset();
439 }
440
441 stop = std::chrono::high_resolution_clock::now();
442 duration = std::chrono::duration_cast<std::chrono::seconds>(stop - start);
443 totalDuration += duration;
444 if (benchmark)
445 {
446 std::cout << " Time for this step: " << duration.count() << " seconds " << std::endl;
447 std::cout << " Elapsed time: " << totalDuration.count() << " seconds (" << std::chrono::duration_cast<std::chrono::minutes>(totalDuration).count() << " minutes)" << std::endl;
448 }
449
450 std::cout << "========================" << std::endl;
451 std::cout << "Reference events: " << eventsRefr.GetValue() << std::endl;
452 std::cout << "Test events: " << eventsTest.GetValue() << std::endl;
453 std::cout << "Ratio: " << eventRatio << std::endl;
454 std::cout << "========================" << std::endl;
455 std::cout << "Tested columns: " << requiredColumns.size() << std::endl;
456 std::cout << "Passed: " << requiredColumns.size() - failedCount << std::endl;
457 std::cout << "Failed: " << failedCount << std::endl;
458 std::cout << "========================" << std::endl;
459 if (checkMissColumns)
460 {
461 std::cout << "Columns only in reference: " << missColNamesRefr.size();
462 for (const auto& column : missColNamesRefr)
463 std::cout << " " << column;
464 std::cout << std::endl;
465 std::cout << "Columns only in test: " << missColNamesTest.size();
466 for (const auto& column : missColNamesTest)
467 std::cout << " " << column;
468 std::cout << std::endl;
469 failedCount += missColNamesRefr.size() + missColNamesTest.size();
470 std::cout << "========================" << std::endl;
471 }
472
473 if (failedCount)
474 {
475 return 1;
476 }
477
478 return 0;
479}
rp
ReadCards * rp
Definition IReadCards.cxx:26
min
#define min(a, b)
Definition cfImp.cxx:40
max
#define max(a, b)
Definition cfImp.cxx:41
intersection
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Definition compareFlatTrees.cxx:25
remainder
std::vector< std::string > remainder(const std::vector< std::string > &v1, const std::vector< std::string > &v2)
Definition compareFlatTrees.cxx:44
verbose
bool verbose
Definition hcg.cxx:75
PixelModuleFeMask_create_db.stop
int stop
Definition PixelModuleFeMask_create_db.py:76
calibdata.valid
list valid
Definition calibdata.py:44
makeegammaturnon.rebin
rebin(binning, data)
Definition makeegammaturnon.py:15
read_hist_ntuple.h1
h1
Definition read_hist_ntuple.py:24
test_FEX_BSEncoding.h2
h2
Definition test_FEX_BSEncoding.py:71
test_pyathena.counter
counter
Definition test_pyathena.py:15
yodamerge_tmp.scale
scale
Definition yodamerge_tmp.py:138
run
int run(int argc, char *argv[])
Definition ttree2hdf5.cxx:28

◆ intersection()

std::vector< std::string > intersection ( std::vector< std::string > & v1,
std::vector< std::string > & v2 )

Definition at line 25 of file compareFlatTrees.cxx.

27{
28 std::vector<std::string> v3;
29
30 std::sort(v1.begin(), v1.end());
31 std::sort(v2.begin(), v2.end());
32
33 std::set_intersection(v1.begin(), v1.end(),
34 v2.begin(), v2.end(),
35 back_inserter(v3));
36 // Alphabetical order while we're at it
37 std::sort(v3.begin(), v3.end(), [](const std::string &a, const std::string &b) -> bool
38 { return a < b; });
39
40 return v3;
41}
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554

◆ remainder()

std::vector< std::string > remainder ( const std::vector< std::string > & v1,
const std::vector< std::string > & v2 )

Definition at line 44 of file compareFlatTrees.cxx.

46{
47 std::vector<std::string> result;
48
49 std::unordered_set<std::string> ignore {v2.begin(), v2.end()};
50 for (const auto& value : v1)
51 {
52 if (ignore.find (value) == ignore.end())
53 result.push_back (value);
54 }
55 return result;
56}
get_generator_info.result
result
Definition get_generator_info.py:21
Generated on for ATLAS Offline Software by doxygen 1.17.0