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