Plugins/Root/BoostHistogramConversionTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "ActsPlugins/Root/HistogramConverter.hpp"
0012
0013 #include <cmath>
0014 #include <random>
0015 #include <string>
0016 #include <vector>
0017
0018 #include <TEfficiency.h>
0019 #include <TH1F.h>
0020 #include <TH2F.h>
0021 #include <TProfile.h>
0022
0023 using namespace Acts;
0024 using namespace Acts::Experimental;
0025 using namespace ActsPlugins;
0026
0027 BOOST_AUTO_TEST_SUITE(HistogramConversionSuite)
0028
0029 namespace {
0030 /// Helper function to test 1D histogram conversion
0031 void testHist1D(const Histogram1& boostHist) {
0032   TH1F* rootHist = toRoot(boostHist);
0033
0034   // Verify metadata
0035   BOOST_CHECK_EQUAL(std::string(rootHist->GetName()), boostHist.name());
0036   BOOST_CHECK_EQUAL(std::string(rootHist->GetTitle()), boostHist.title());
0037
0038   // Verify number of bins
0039   const auto& bh = boostHist.histogram();
0040   BOOST_CHECK_EQUAL(rootHist->GetNbinsX(), static_cast<int>(bh.axis(0).size()));
0041
0042   // Verify bin contents match
0043   for (int i = 1; i <= rootHist->GetNbinsX(); ++i) {
0044     BOOST_CHECK_CLOSE(rootHist->GetBinContent(i),
0045                       static_cast<double>(bh.at(i - 1)), 1e-10);
0046   }
0047
0048   delete rootHist;
0049 }
0050
0051 /// Helper function to test 2D histogram conversion
0052 void testHist2D(const Histogram2& boostHist) {
0053   TH2F* rootHist = toRoot(boostHist);
0054
0055   // Verify metadata
0056   BOOST_CHECK_EQUAL(std::string(rootHist->GetName()), boostHist.name());
0057   BOOST_CHECK_EQUAL(std::string(rootHist->GetTitle()), boostHist.title());
0058
0059   // Verify number of bins
0060   const auto& bh = boostHist.histogram();
0061   BOOST_CHECK_EQUAL(rootHist->GetNbinsX(), static_cast<int>(bh.axis(0).size()));
0062   BOOST_CHECK_EQUAL(rootHist->GetNbinsY(), static_cast<int>(bh.axis(1).size()));
0063
0064   // Verify bin contents match
0065   for (int i = 1; i <= rootHist->GetNbinsX(); ++i) {
0066     for (int j = 1; j <= rootHist->GetNbinsY(); ++j) {
0067       BOOST_CHECK_CLOSE(rootHist->GetBinContent(i, j),
0068                         static_cast<double>(bh.at(i - 1, j - 1)), 1e-10);
0069     }
0070   }
0071
0072   delete rootHist;
0073 }
0074 }  // namespace
0075
0076 BOOST_AUTO_TEST_CASE(Conversion_EmptyHistogram) {
0077   auto axis = AxisVariant(BoostRegularAxis(10, -10.0, 10.0, "x"));
0078   Histogram1 boostHist("empty", "Empty Histogram", {axis});
0079
0080   TH1F* rootHist = toRoot(boostHist);
0081
0082   BOOST_CHECK_EQUAL(rootHist->GetNbinsX(), 10);
0083   for (int i = 1; i <= rootHist->GetNbinsX(); ++i) {
0084     BOOST_CHECK_EQUAL(rootHist->GetBinContent(i), 0.0);
0085   }
0086
0087   delete rootHist;
0088 }
0089
0090 BOOST_AUTO_TEST_CASE(Conversion_Boost1D_to_ROOT_UniformBinning) {
0091   auto axis = AxisVariant(BoostRegularAxis(10, 0.0, 10.0, "x [cm]"));
0092   Histogram1 boostHist("test_hist", "Test Histogram", {axis});
0093
0094   // Fill with 100 random values
0095   std::mt19937 rng(42);
0096   std::normal_distribution<double> dist(5.0, 1.0);
0097   for (int i = 0; i < 100; ++i) {
0098     boostHist.fill({dist(rng)});
0099   }
0100
0101   testHist1D(boostHist);
0102 }
0103
0104 BOOST_AUTO_TEST_CASE(Conversion_Boost1D_to_ROOT_VariableBinning) {
0105   std::vector<double> edges = {0.0, 0.5, 1.0, 2.0, 5.0, 10.0};
0106   auto axis = BoostVariableAxis(edges, "eta");
0107   Histogram1 boostHist("res_eta", "Residual vs Eta", {axis});
0108
0109   // Fill bins with different counts
0110   boostHist.fill({0.3});
0111   boostHist.fill({0.3});
0112   boostHist.fill({1.5});
0113   boostHist.fill({3.0});
0114
0115   testHist1D(boostHist);
0116 }
0117
0118 BOOST_AUTO_TEST_CASE(Conversion_Boost2D_to_ROOT) {
0119   auto xAxis = BoostRegularAxis(10, -2.5, 2.5, "eta");
0120   auto yAxis = BoostRegularAxis(10, -5.0, 5.0, "residual");
0121   Histogram2 boostHist("res_vs_eta", "Residual vs Eta", {xAxis, yAxis});
0122
0123   // Fill with 100 2D Gaussian values
0124   std::mt19937 rng(456);
0125   std::normal_distribution<double> etaDist(0.0, 1.0);
0126   std::normal_distribution<double> resDist(0.0, 2.0);
0127   for (int i = 0; i < 100; ++i) {
0128     boostHist.fill({etaDist(rng), resDist(rng)});
0129   }
0130
0131   testHist2D(boostHist);
0132 }
0133
0134 BOOST_AUTO_TEST_CASE(Conversion_Boost2D_to_ROOT_VariableBinning) {
0135   std::vector<double> etaEdges = {-2.5, -1.5, -0.5, 0.5, 1.5, 2.5};
0136   std::vector<double> ptEdges = {0.5, 1.0, 2.0, 5.0, 10.0};
0137
0138   auto xAxis = BoostVariableAxis(etaEdges, "eta");
0139   auto yAxis = BoostVariableAxis(ptEdges, "pT [GeV]");
0140   Histogram2 boostHist("eff_vs_eta_pt", "Efficiency vs Eta and pT",
0141                        {xAxis, yAxis});
0142
0143   boostHist.fill({0.0, 3.0});
0144   boostHist.fill({-1.0, 7.0});
0145   boostHist.fill({2.0, 1.5});
0146
0147   testHist2D(boostHist);
0148 }
0149
0150 BOOST_AUTO_TEST_CASE(Conversion_BoostProfile_to_TProfile) {
0151   auto xAxis = BoostRegularAxis(10, -2.5, 2.5, "eta");
0152   ProfileHistogram1 profile("res_mean_vs_eta", "Mean Residual vs Eta", {xAxis},
0153                             "residual [mm]");
0154
0155   std::map<double, double> expectedMeans;
0156
0157   // Fill with known values to check mean calculation
0158   profile.fill({-2.0}, 1.0);
0159   profile.fill({-2.0}, 3.0);
0160   expectedMeans[-2.0] = 2.0;
0161
0162   profile.fill({0.0}, 5.0);
0163   profile.fill({0.0}, 7.0);
0164   expectedMeans[0.0] = 6.0;
0165
0166   profile.fill({2.0}, 9.0);
0167   profile.fill({2.0}, 11.0);
0168   profile.fill({2.0}, 13.0);
0169   expectedMeans[2.0] = 11.0;
0170
0171   TProfile* rootProfile = toRoot(profile);
0172
0173   // Verify metadata
0174   BOOST_CHECK_EQUAL(std::string(rootProfile->GetName()), "res_mean_vs_eta");
0175   BOOST_CHECK_EQUAL(std::string(rootProfile->GetTitle()),
0176                     "Mean Residual vs Eta");
0177   BOOST_CHECK_EQUAL(std::string(rootProfile->GetXaxis()->GetTitle()), "eta");
0178   BOOST_CHECK_EQUAL(std::string(rootProfile->GetYaxis()->GetTitle()),
0179                     "residual [mm]");
0180
0181   // Verify binning
0182   BOOST_CHECK_EQUAL(rootProfile->GetNbinsX(), 10);
0183
0184   // Verify mean values in filled bins
0185   const auto& bh = profile.histogram();
0186
0187   for (auto x : {-2.0, 0.0, 2.0}) {
0188     auto binIdx = bh.axis(0).index(x);
0189     BOOST_CHECK_CLOSE(rootProfile->GetBinContent(binIdx + 1),
0190                       bh.at(binIdx).value(), 1e-6);
0191     BOOST_CHECK_CLOSE(rootProfile->GetBinContent(binIdx + 1),
0192                       expectedMeans.at(x), 1e-6);
0193   }
0194
0195   delete rootProfile;
0196 }
0197
0198 BOOST_AUTO_TEST_CASE(Conversion_BoostProfile_to_TProfile_WithErrors) {
0199   auto xAxis = BoostRegularAxis(5, -2.5, 2.5, "eta");
0200
0201   // Create ACTS ProfileHistogram
0202   ProfileHistogram1 actsProfile("profile_test", "Profile Test", {xAxis},
0203                                 "value [units]");
0204
0205   // Create ROOT TProfile with identical binning for comparison
0206   const auto& bh = actsProfile.histogram();
0207   const auto& axis = bh.axis(0);
0208   std::vector<double> edges(axis.size() + 1);
0209   for (int i = 0; i < axis.size(); ++i) {
0210     edges[i] = axis.bin(i).lower();
0211   }
0212   edges.back() = axis.bin(axis.size() - 1).upper();
0213   TProfile rootProfile("root_profile", "ROOT Profile", 5, edges.data());
0214   rootProfile.Sumw2();
0215
0216   // Bin 1: mean=12, n=3
0217   for (auto y : {10.0, 12.0, 14.0}) {
0218     rootProfile.Fill(-2.0, y);
0219     actsProfile.fill({-2.0}, y);
0220   }
0221
0222   // Bin 2: mean=6, n=2
0223   for (auto y : {5.0, 7.0}) {
0224     rootProfile.Fill(0.0, y);
0225     actsProfile.fill({0.0}, y);
0226   }
0227
0228   // Bin 3: mean=23, n=4
0229   for (auto y : {20.0, 22.0, 24.0, 26.0}) {
0230     rootProfile.Fill(1.5, y);
0231     actsProfile.fill({1.5}, y);
0232   }
0233
0234   // Convert ACTS profile to ROOT
0235   TProfile* convertedProfile = toRoot(actsProfile);
0236
0237   // Compare binning
0238   BOOST_CHECK_EQUAL(convertedProfile->GetNbinsX(), rootProfile.GetNbinsX());
0239
0240   // Compare each filled bin
0241   for (int i = 0; i < axis.size(); ++i) {
0242     int rootBin = i + 1;
0243     const auto& acc = bh.at(i);
0244
0245     // Check count
0246     BOOST_CHECK_EQUAL(convertedProfile->GetBinEntries(rootBin),
0247                       rootProfile.GetBinEntries(rootBin));
0248     BOOST_CHECK_CLOSE(convertedProfile->GetBinEntries(rootBin), acc.count(),
0249                       1e-10);
0250
0251     // Check mean value
0252     BOOST_CHECK_CLOSE(convertedProfile->GetBinContent(rootBin),
0253                       rootProfile.GetBinContent(rootBin), 1e-6);
0254     BOOST_CHECK_CLOSE(convertedProfile->GetBinContent(rootBin), acc.value(),
0255                       1e-6);
0256
0257     // Check errors
0258     double rootError = rootProfile.GetBinError(rootBin);
0259     double convertedError = convertedProfile->GetBinError(rootBin);
0260     BOOST_CHECK_CLOSE(convertedError, rootError, 1e-6);
0261   }
0262
0263   delete convertedProfile;
0264 }
0265
0266 BOOST_AUTO_TEST_CASE(Conversion_Efficiency1D_to_TEfficiency) {
0267   auto axis = BoostRegularAxis(10, -3.0, 3.0, "eta");
0268   Efficiency1 eff("eff_vs_eta", "Efficiency vs Eta", {axis});
0269
0270   // Fill with known pass/fail patterns
0271   // Bin at eta=0.5: 3 accepted, 1 failed (75% efficiency)
0272   for (auto v : {true, true, true, false}) {
0273     eff.fill({0.5}, v);
0274   }
0275
0276   // Bin at eta=-1.5: 3 accepted, 3 failed (50% efficiency)
0277   for (auto v : {true, true, true, false, false, false}) {
0278     eff.fill({-1.5}, v);
0279   }
0280
0281   TEfficiency* rootEff = toRoot(eff);
0282
0283   // Verify metadata
0284   BOOST_CHECK_EQUAL(std::string(rootEff->GetName()), "eff_vs_eta");
0285   BOOST_CHECK_EQUAL(std::string(rootEff->GetTitle()), "Efficiency vs Eta");
0286
0287   // Verify binning
0288   BOOST_CHECK_EQUAL(rootEff->GetTotalHistogram()->GetNbinsX(), 10);
0289
0290   // Verify efficiency values
0291   const auto& accepted = eff.acceptedHistogram();
0292   const auto& total = eff.totalHistogram();
0293
0294   for (auto x : {0.5, -1.5}) {
0295     auto binIdx = accepted.axis(0).index(x);
0296     double expectedEff = static_cast<double>(accepted.at(binIdx)) /
0297                          static_cast<double>(total.at(binIdx));
0298     BOOST_CHECK_CLOSE(rootEff->GetEfficiency(binIdx + 1), expectedEff, 1e-6);
0299   }
0300
0301   delete rootEff;
0302 }
0303
0304 BOOST_AUTO_TEST_CASE(Conversion_Efficiency2D_to_TEfficiency) {
0305   auto xAxis = BoostRegularAxis(5, -2.5, 2.5, "eta");
0306   auto yAxis = BoostRegularAxis(5, 0.0, 5.0, "pt");
0307   Efficiency2 eff("eff_vs_eta_pt", "Efficiency vs Eta and pT", {xAxis, yAxis});
0308
0309   // Fill bin (0.0, 2.5): 3 accepted, 1 failed (75% efficiency)
0310   for (auto v : {true, true, true, false}) {
0311     eff.fill({0.0, 2.5}, v);
0312   }
0313
0314   // Fill bin (-1.5, 1.5): 3 accepted, 3 failed (50% efficiency)
0315   for (auto v : {true, true, true, false, false, false}) {
0316     eff.fill({-1.5, 1.5}, v);
0317   }
0318
0319   TEfficiency* rootEff = toRoot(eff);
0320
0321   // Verify metadata
0322   BOOST_CHECK_EQUAL(std::string(rootEff->GetName()), "eff_vs_eta_pt");
0323   BOOST_CHECK_EQUAL(std::string(rootEff->GetTitle()),
0324                     "Efficiency vs Eta and pT");
0325
0326   // Verify 2D binning
0327   BOOST_CHECK_EQUAL(rootEff->GetTotalHistogram()->GetNbinsX(), 5);
0328   BOOST_CHECK_EQUAL(rootEff->GetTotalHistogram()->GetNbinsY(), 5);
0329
0330   // Verify efficiency values
0331   const auto& accepted = eff.acceptedHistogram();
0332   const auto& total = eff.totalHistogram();
0333
0334   using P = std::pair<double, double>;
0335   for (auto coords : {P{0.0, 2.5}, P{-1.5, 1.5}}) {
0336     auto xIdx = accepted.axis(0).index(coords.first);
0337     auto yIdx = accepted.axis(1).index(coords.second);
0338     double expectedEff = static_cast<double>(accepted.at(xIdx, yIdx)) /
0339                          static_cast<double>(total.at(xIdx, yIdx));
0340     int globalBin = rootEff->GetGlobalBin(static_cast<int>(xIdx + 1),
0341                                           static_cast<int>(yIdx + 1));
0342     BOOST_CHECK_CLOSE(rootEff->GetEfficiency(globalBin), expectedEff, 1e-6);
0343   }
0344
0345   delete rootEff;
0346 }
0347
0348 BOOST_AUTO_TEST_SUITE_END()