Plugins/Root/RootHistogramFitTests.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 "Acts/Utilities/Histogram.hpp"
0012 #include "ActsPlugins/Root/HistogramConverter.hpp"
0013 #include "ActsPlugins/Root/RootHistogramFit.hpp"
0014
0015 #include <cstdint>
0016 #include <random>
0017 #include <string>
0018 #include <vector>
0019
0020 #include <TH1.h>
0021 #include <TH2.h>
0022
0023 using namespace Acts::Experimental;
0024 using ActsPlugins::RootHistogramFit;
0025 using ActsPlugins::toRoot;
0026
0027 namespace {
0028
0029 Histogram1 makeHistogram(const std::string& name, int nBins, double xMin,
0030                          double xMax) {
0031   auto axis = AxisVariant(BoostRegularAxis(nBins, xMin, xMax, "x"));
0032   return Histogram1(name, name, {axis});
0033 }
0034
0035 /// Sample `count` entries from N(mean, sigma) into a fresh histogram
0036 Histogram1 sampled(const std::string& name, std::size_t count, double mean,
0037                    double sigma, int nBins, double xMin, double xMax,
0038                    std::uint32_t seed) {
0039   Histogram1 hist = makeHistogram(name, nBins, xMin, xMax);
0040
0041   std::mt19937 generator(seed);
0042   std::normal_distribution<double> distribution(mean, sigma);
0043   for (std::size_t i = 0; i < count; ++i) {
0044     hist.fill({distribution(generator)});
0045   }
0046
0047   return hist;
0048 }
0049
0050 }  // namespace
0051
0052 BOOST_AUTO_TEST_SUITE(RootHistogramFitSuite)
0053
0054 BOOST_AUTO_TEST_CASE(SingleFit_RecoversKnownParameters) {
0055   const RootHistogramFit fitter;
0056   const Histogram1 hist =
0057       sampled("single", 100000, 0.0, 1.0, 100, -8.0, 8.0, 101);
0058
0059   const auto result = fitter(hist);
0060   BOOST_REQUIRE(result.has_value());
0061
0062   const auto& [mean, sigma, meanError, sigmaError] = *result;
0063   BOOST_CHECK_SMALL(mean, 0.05);
0064   BOOST_CHECK_CLOSE(sigma, 1.0, 5.0);
0065   BOOST_CHECK_GT(meanError, 0.0);
0066   BOOST_CHECK_GT(sigmaError, 0.0);
0067 }
0068
0069 BOOST_AUTO_TEST_CASE(RestrictedRange_Works) {
0070   const RootHistogramFit fitter;
0071   const Histogram1 hist =
0072       sampled("restricted", 30000, 0.4, 1.2, 80, -8.0, 8.0, 202);
0073
0074   const auto result = fitter(hist, RootHistogramFit::Range{-2.6, 3.4});
0075   BOOST_REQUIRE(result.has_value());
0076
0077   const auto& [mean, sigma, meanError, sigmaError] = *result;
0078   BOOST_CHECK_CLOSE(mean, 0.4, 10.0);
0079   BOOST_CHECK_CLOSE(sigma, 1.2, 10.0);
0080 }
0081
0082 BOOST_AUTO_TEST_CASE(DegenerateInputs_Fail) {
0083   const RootHistogramFit fitter;
0084   auto axis = AxisVariant(BoostRegularAxis(20, -5.0, 5.0, "x"));
0085
0086   Histogram1 empty("empty", "empty", {axis});
0087   BOOST_CHECK(!fitter(empty).has_value());
0088
0089   Histogram1 spike("spike", "spike", {axis});
0090   spike.setBinContent({10}, 500.0);
0091   BOOST_CHECK(!fitter(spike).has_value());
0092 }
0093
0094 BOOST_AUTO_TEST_CASE(Histogram1D_ConvertsWithErrors) {
0095   std::vector<double> edges = {0.0, 1.0, 3.0, 7.0};
0096   auto axis = AxisVariant(BoostVariableAxis(edges, "eta"));
0097   Histogram1 hist("resmean_d0_vs_eta", "Mean", {axis});
0098
0099   hist.setBin({0}, 1.5, 0.25);
0100   hist.setBin({1}, -2.5, 0.5);
0101   // Bin 2 deliberately left untouched
0102
0103   const auto rootHist = toRoot(hist);
0104   BOOST_REQUIRE(rootHist != nullptr);
0105   BOOST_CHECK_EQUAL(rootHist->GetName(), "resmean_d0_vs_eta");
0106   BOOST_CHECK_EQUAL(rootHist->GetTitle(), "Mean");
0107   BOOST_CHECK_EQUAL(rootHist->GetNbinsX(), 3);
0108   BOOST_CHECK_EQUAL(std::string(rootHist->GetXaxis()->GetTitle()), "eta");
0109
0110   BOOST_CHECK_CLOSE(rootHist->GetBinContent(1), 1.5, 1e-4);
0111   BOOST_CHECK_CLOSE(rootHist->GetBinError(1), 0.25, 1e-4);
0112   BOOST_CHECK_CLOSE(rootHist->GetBinContent(2), -2.5, 1e-4);
0113   BOOST_CHECK_CLOSE(rootHist->GetBinError(2), 0.5, 1e-4);
0114   BOOST_CHECK_EQUAL(rootHist->GetBinContent(3), 0.0);
0115   BOOST_CHECK_EQUAL(rootHist->GetBinError(3), 0.0);
0116
0117   // Variable binning must be carried over
0118   BOOST_CHECK_CLOSE(rootHist->GetXaxis()->GetBinLowEdge(1), 0.0, 1e-6);
0119   BOOST_CHECK_CLOSE(rootHist->GetXaxis()->GetBinUpEdge(3), 7.0, 1e-6);
0120 }
0121
0122 BOOST_AUTO_TEST_CASE(Histogram2D_ConvertsWithErrors) {
0123   auto xAxis = AxisVariant(BoostRegularAxis(2, 0.0, 2.0, "eta"));
0124   auto yAxis = AxisVariant(BoostRegularAxis(3, 0.0, 3.0, "pt"));
0125   Histogram2 hist("reswidth_d0_vs_eta_pt", "Width", {xAxis, yAxis});
0126
0127   hist.setBin({1, 2}, 0.75, 0.1);
0128
0129   const auto rootHist = toRoot(hist);
0130   BOOST_REQUIRE(rootHist != nullptr);
0131   BOOST_CHECK_EQUAL(rootHist->GetNbinsX(), 2);
0132   BOOST_CHECK_EQUAL(rootHist->GetNbinsY(), 3);
0133   BOOST_CHECK_EQUAL(std::string(rootHist->GetXaxis()->GetTitle()), "eta");
0134   BOOST_CHECK_EQUAL(std::string(rootHist->GetYaxis()->GetTitle()), "pt");
0135
0136   BOOST_CHECK_CLOSE(rootHist->GetBinContent(2, 3), 0.75, 1e-4);
0137   BOOST_CHECK_CLOSE(rootHist->GetBinError(2, 3), 0.1, 1e-4);
0138   BOOST_CHECK_EQUAL(rootHist->GetBinContent(1, 1), 0.0);
0139 }
0140
0141 BOOST_AUTO_TEST_SUITE_END()