Root/src/HistogramConverter.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 "ActsPlugins/Root/HistogramConverter.hpp"
0010
0011 #include <cstddef>
0012 #include <format>
0013 #include <vector>
0014
0015 #include <TEfficiency.h>
0016 #include <TFitResult.h>
0017 #include <TGraphAsymmErrors.h>
0018 #include <TH1F.h>
0019 #include <TH2F.h>
0020 #include <TH3F.h>
0021 #include <TProfile.h>
0022 #include <boost/histogram/accumulators/weighted_sum.hpp>
0023
0024 using namespace Acts::Experimental;
0025
0026 std::unique_ptr<TH1F> ActsPlugins::toRoot(const Histogram1& boostHist) {
0027   const auto& bh = boostHist.histogram();
0028   const auto& axis = bh.axis(0);
0029
0030   // Extract bin edges from boost histogram axis
0031   std::vector<double> edges = extractBinEdges(axis);
0032
0033   // Create ROOT histogram with variable binning
0034   auto rootHist = std::make_unique<TH1F>(boostHist.name().c_str(),
0035                                          boostHist.title().c_str(), axis.size(),
0036                                          edges.data());
0037
0038   // Copy bin contents from boost to ROOT
0039   for (auto&& x : boost::histogram::indexed(bh)) {
0040     // Dereference to get bin content
0041     double content = *x;
0042
0043     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0044     int rootBinIndex = x.index(0) + 1;
0045     rootHist->SetBinContent(rootBinIndex, content);
0046   }
0047
0048   // Set axis titles from axis metadata
0049   rootHist->GetXaxis()->SetTitle(axis.metadata().c_str());
0050
0051   return rootHist;
0052 }
0053
0054 std::unique_ptr<TH2F> ActsPlugins::toRoot(const Histogram2& boostHist) {
0055   const auto& bh = boostHist.histogram();
0056   const auto& xAxis = bh.axis(0);
0057   const auto& yAxis = bh.axis(1);
0058
0059   // Extract bin edges from X axis
0060   std::vector<double> xEdges = extractBinEdges(xAxis);
0061
0062   // Extract bin edges from Y axis
0063   std::vector<double> yEdges = extractBinEdges(yAxis);
0064
0065   // Create ROOT histogram with 2D variable binning
0066   auto rootHist = std::make_unique<TH2F>(
0067       boostHist.name().c_str(), boostHist.title().c_str(), xAxis.size(),
0068       xEdges.data(), yAxis.size(), yEdges.data());
0069
0070   // Copy bin contents from boost to ROOT
0071   for (auto&& x : boost::histogram::indexed(bh)) {
0072     // Dereference to get bin content
0073     double content = *x;
0074
0075     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0076     // indexed() gives us 0-based bin indices for each axis
0077     int rootXBin = x.index(0) + 1;
0078     int rootYBin = x.index(1) + 1;
0079     rootHist->SetBinContent(rootXBin, rootYBin, content);
0080   }
0081
0082   // Set axis titles from axis metadata
0083   rootHist->GetXaxis()->SetTitle(xAxis.metadata().c_str());
0084   rootHist->GetYaxis()->SetTitle(yAxis.metadata().c_str());
0085
0086   return rootHist;
0087 }
0088
0089 std::unique_ptr<TH3F> ActsPlugins::toRoot(const Histogram3& boostHist) {
0090   const auto& bh = boostHist.histogram();
0091   const auto& xAxis = bh.axis(0);
0092   const auto& yAxis = bh.axis(1);
0093   const auto& zAxis = bh.axis(2);
0094
0095   // Extract bin edges from X axis
0096   std::vector<double> xEdges = extractBinEdges(xAxis);
0097
0098   // Extract bin edges from Y axis
0099   std::vector<double> yEdges = extractBinEdges(yAxis);
0100
0101   // Extract bin edges from Z axis
0102   std::vector<double> zEdges = extractBinEdges(zAxis);
0103
0104   // Create ROOT histogram with 3D variable binning
0105   auto rootHist = std::make_unique<TH3F>(
0106       boostHist.name().c_str(), boostHist.title().c_str(), xAxis.size(),
0107       xEdges.data(), yAxis.size(), yEdges.data(), zAxis.size(), zEdges.data());
0108
0109   // Copy bin contents from boost to ROOT
0110   for (auto&& x : boost::histogram::indexed(bh)) {
0111     // Dereference to get bin content
0112     double content = *x;
0113
0114     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0115     // indexed() gives us 0-based bin indices for each axis
0116     int rootXBin = x.index(0) + 1;
0117     int rootYBin = x.index(1) + 1;
0118     int rootZBin = x.index(2) + 1;
0119     rootHist->SetBinContent(rootXBin, rootYBin, rootZBin, content);
0120   }
0121
0122   // Set axis titles from axis metadata
0123   rootHist->GetXaxis()->SetTitle(xAxis.metadata().c_str());
0124   rootHist->GetYaxis()->SetTitle(yAxis.metadata().c_str());
0125   rootHist->GetZaxis()->SetTitle(zAxis.metadata().c_str());
0126
0127   return rootHist;
0128 }
0129
0130 std::unique_ptr<TProfile> ActsPlugins::toRoot(
0131     const ProfileHistogram1& boostProfile) {
0132   const auto& bh = boostProfile.histogram();
0133   const auto& axis = bh.axis(0);
0134
0135   // Extract bin edges from boost histogram axis
0136   std::vector<double> edges = extractBinEdges(axis);
0137
0138   // Create ROOT TProfile with variable binning
0139   auto rootProfile = std::make_unique<TProfile>(boostProfile.name().c_str(),
0140                                                 boostProfile.title().c_str(),
0141                                                 axis.size(), edges.data());
0142
0143   // Enable sum of weights squared storage for proper error calculation
0144   rootProfile->Sumw2();
0145
0146   // Copy data from boost profile to ROOT profile
0147   // - count(): number of fills
0148   // - value(): mean of y-values
0149   // - variance(): sample variance = sum((y - mean)^2) / (n - 1)
0150
0151   using Accumulator = boost::histogram::accumulators::mean<double>;
0152
0153   for (auto&& x : boost::histogram::indexed(bh)) {
0154     const Accumulator& acc = *x;
0155
0156     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0157     int rootBinIndex = x.index(0) + 1;
0158
0159     double count = acc.count();
0160
0161     if (count == 0) {
0162       continue;
0163     }
0164     double mean = acc.value();
0165     double variance = acc.variance();  // Sample variance from boost
0166
0167     // ROOT TProfile internally stores the following data (all weights=1):
0168     // - fArray[bin]      = sum of (w*x)    = count * mean
0169     // - fBinEntries[bin] = sum of w        = count
0170     // - fSumw2[bin]      = sum of (w*x)^2  = sum of x^2
0171     // - fBinSumw2[bin]   = sum of w^2      = count
0172
0173     double sum = mean * count;
0174
0175     // To compute sum of y^2 from sample variance s^2:
0176     // s^2 = Sum((x - μ)^2) / (n-1) = (Sum(x^2) - n*µ^2) / (n-1)
0177     // Sum(x^2) = (n-1) * s^2 + n*µ^2
0178     double sumOfSquares = (count - 1.0) * variance + count * mean * mean;
0179
0180     // Set bin content (sum) and entries via public interface
0181     rootProfile->SetBinContent(rootBinIndex, sum);
0182     rootProfile->SetBinEntries(rootBinIndex, count);
0183
0184     // Access internal arrays to set sum of squares for error calculation
0185     // This is necessary because ROOT has no public API to set these arrays
0186     TArrayD* sumw2 = rootProfile->GetSumw2();        // fSumw2 array
0187     TArrayD* binSumw2 = rootProfile->GetBinSumw2();  // fBinSumw2 array
0188
0189     assert(sumw2 && "Sumw2 is null");
0190     assert(binSumw2 && "BinSumw2 is null");
0191
0192     // Set sum of (weight * value)^2 = sum of y^2 for unweighted data
0193     sumw2->fArray[rootBinIndex] = sumOfSquares;
0194     // Set sum of weights^2 = count for unweighted data (all weights = 1)
0195     binSumw2->fArray[rootBinIndex] = count;
0196   }
0197
0198   // Set X axis title from axis metadata
0199   rootProfile->GetXaxis()->SetTitle(axis.metadata().c_str());
0200   // Set Y axis title from sampleAxisTitle member
0201   rootProfile->GetYaxis()->SetTitle(boostProfile.sampleAxisTitle().c_str());
0202
0203   return rootProfile;
0204 }
0205
0206 std::unique_ptr<TEfficiency> ActsPlugins::toRoot(const Efficiency1& boostEff) {
0207   const auto& accepted = boostEff.acceptedHistogram();
0208   const auto& total = boostEff.totalHistogram();
0209   const auto& axis = accepted.axis(0);
0210
0211   // Extract bin edges
0212   std::vector<double> edges = extractBinEdges(axis);
0213
0214   // Create accepted and total TH1F histograms
0215   auto acceptedHist =
0216       std::make_unique<TH1F>((boostEff.name() + "_accepted").c_str(), "Passed",
0217                              axis.size(), edges.data());
0218
0219   auto totalHist = std::make_unique<TH1F>((boostEff.name() + "_total").c_str(),
0220                                           "Total", axis.size(), edges.data());
0221
0222   // Fill histograms with counts
0223   for (int i = 0; i < axis.size(); ++i) {
0224     auto acceptedCount = static_cast<double>(accepted.at(i));
0225     auto totalCount = static_cast<double>(total.at(i));
0226
0227     acceptedHist->SetBinContent(i + 1, acceptedCount);
0228     totalHist->SetBinContent(i + 1, totalCount);
0229   }
0230
0231   // Set X axis title from axis metadata
0232   acceptedHist->GetXaxis()->SetTitle(axis.metadata().c_str());
0233   totalHist->GetXaxis()->SetTitle(axis.metadata().c_str());
0234
0235   // Set Y axis titles; use "Efficiency" for total histogram since TEfficiency
0236   // will inherit this title for the efficiency graph
0237   totalHist->GetYaxis()->SetTitle("Efficiency");
0238
0239   // Create TEfficiency from the two histograms
0240   // TEfficiency takes ownership of the histograms
0241   auto rootEff = std::make_unique<TEfficiency>(*acceptedHist, *totalHist);
0242   rootEff->SetName(boostEff.name().c_str());
0243   rootEff->SetTitle(boostEff.title().c_str());
0244
0245   return rootEff;
0246 }
0247
0248 std::unique_ptr<TEfficiency> ActsPlugins::toRoot(const Efficiency2& boostEff) {
0249   const auto& accepted = boostEff.acceptedHistogram();
0250   const auto& total = boostEff.totalHistogram();
0251   const auto& xAxis = accepted.axis(0);
0252   const auto& yAxis = accepted.axis(1);
0253
0254   // Extract X bin edges
0255   std::vector<double> xEdges = extractBinEdges(xAxis);
0256
0257   // Extract Y bin edges
0258   std::vector<double> yEdges = extractBinEdges(yAxis);
0259
0260   // Create accepted and total TH2F histograms
0261   auto acceptedHist = std::make_unique<TH2F>(
0262       (boostEff.name() + "_accepted").c_str(), "Accepted", xAxis.size(),
0263       xEdges.data(), yAxis.size(), yEdges.data());
0264
0265   auto totalHist = std::make_unique<TH2F>((boostEff.name() + "_total").c_str(),
0266                                           "Total", xAxis.size(), xEdges.data(),
0267                                           yAxis.size(), yEdges.data());
0268
0269   // Fill histograms with counts
0270   for (int i = 0; i < xAxis.size(); ++i) {
0271     for (int j = 0; j < yAxis.size(); ++j) {
0272       auto acceptedCount = static_cast<double>(accepted.at(i, j));
0273       auto totalCount = total.at(i, j);
0274
0275       acceptedHist->SetBinContent(i + 1, j + 1, acceptedCount);
0276       totalHist->SetBinContent(i + 1, j + 1, totalCount);
0277     }
0278   }
0279
0280   // Set X axis title from axis metadata
0281   acceptedHist->GetXaxis()->SetTitle(xAxis.metadata().c_str());
0282   totalHist->GetXaxis()->SetTitle(xAxis.metadata().c_str());
0283
0284   // Set Y axis title from axis metadata
0285   acceptedHist->GetYaxis()->SetTitle(yAxis.metadata().c_str());
0286   totalHist->GetYaxis()->SetTitle(yAxis.metadata().c_str());
0287
0288   // Set Z axis titles; use "Efficiency" for total histogram since TEfficiency
0289   // will inherit this title for the efficiency graph
0290   totalHist->GetZaxis()->SetTitle("Efficiency");
0291
0292   // Create TEfficiency from the two histograms
0293   // TEfficiency takes ownership of the histograms
0294   auto rootEff = std::make_unique<TEfficiency>(*acceptedHist, *totalHist);
0295   rootEff->SetName(boostEff.name().c_str());
0296   rootEff->SetTitle(boostEff.title().c_str());
0297
0298   return rootEff;
0299 }
0300
0301 std::tuple<std::unique_ptr<TH1F>, std::unique_ptr<TH1F>, double>
0302 ActsPlugins::extractMeanWidthProfiles(const TH2F& hist2d,
0303                                       const std::string& meanName,
0304                                       const std::string& widthName,
0305                                       const int minEntriesForFit,
0306                                       const std::string& fitOption,
0307                                       const Acts::Logger& logger) {
0308   const int nBinsX = hist2d.GetNbinsX();
0309
0310   // Create mean and width histograms with same X binning as the 2D histogram
0311   auto meanHist = std::make_unique<TH1F>(
0312       meanName.c_str(), (std::string(hist2d.GetTitle()) + " mean").c_str(),
0313       nBinsX, hist2d.GetXaxis()->GetXmin(), hist2d.GetXaxis()->GetXmax());
0314   auto widthHist = std::make_unique<TH1F>(
0315       widthName.c_str(), (std::string(hist2d.GetTitle()) + " width").c_str(),
0316       nBinsX, hist2d.GetXaxis()->GetXmin(), hist2d.GetXaxis()->GetXmax());
0317
0318   // Copy X axis bin edges for variable binning
0319   if (hist2d.GetXaxis()->GetXbins()->GetSize() > 0) {
0320     meanHist->SetBins(nBinsX, hist2d.GetXaxis()->GetXbins()->GetArray());
0321     widthHist->SetBins(nBinsX, hist2d.GetXaxis()->GetXbins()->GetArray());
0322   }
0323
0324   // Project each X bin and extract mean/width via Gaussian fit
0325   int fitFailures = 0;
0326   for (int i = 1; i <= nBinsX; ++i) {
0327     const auto proj = std::unique_ptr<TH1D>(hist2d.ProjectionY(
0328         std::format("{}_projy_bin_{}", hist2d.GetName(), i).c_str(), i, i));
0329
0330     if (proj->GetEntries() < minEntriesForFit) {
0331       continue;
0332     }
0333
0334     const TFitResultPtr r = proj->Fit("gaus", fitOption.c_str());
0335     if ((r.Get() == nullptr) || ((r->Status() % 1000) != 0)) {
0336       ++fitFailures;
0337       ACTS_DEBUG("Failed to fit Gaussian for bin "
0338                  << i << ": status "
0339                  << (r.Get() != nullptr ? r->Status() : -1));
0340       continue;
0341     }
0342
0343     // Fill mean
0344     meanHist->SetBinContent(i, r->Parameter(1));
0345     meanHist->SetBinError(i, r->ParError(1));
0346
0347     // Fill width (sigma)
0348     widthHist->SetBinContent(i, r->Parameter(2));
0349     widthHist->SetBinError(i, r->ParError(2));
0350   }
0351   const double fitFailureFraction =
0352       (nBinsX > 0) ? static_cast<double>(fitFailures) / nBinsX : 0;
0353
0354   meanHist->GetXaxis()->SetTitle(hist2d.GetXaxis()->GetTitle());
0355   meanHist->GetYaxis()->SetTitle(hist2d.GetYaxis()->GetTitle());
0356
0357   widthHist->GetXaxis()->SetTitle(hist2d.GetXaxis()->GetTitle());
0358   widthHist->GetYaxis()->SetTitle(hist2d.GetYaxis()->GetTitle());
0359
0360   return {std::move(meanHist), std::move(widthHist), fitFailureFraction};
0361 }
0362
0363 std::tuple<std::unique_ptr<TH2F>, std::unique_ptr<TH2F>, double>
0364 ActsPlugins::extractMeanWidthProfiles(const TH3F& hist3d,
0365                                       const std::string& meanName,
0366                                       const std::string& widthName,
0367                                       const int minEntriesForFit,
0368                                       const std::string& fitOption,
0369                                       const Acts::Logger& logger) {
0370   const int nBinsX = hist3d.GetNbinsX();
0371   const int nBinsY = hist3d.GetNbinsY();
0372
0373   // Create output histograms with same XY binning as input
0374   auto meanHist = std::make_unique<TH2F>(
0375       meanName.c_str(), (std::string(hist3d.GetTitle()) + " mean").c_str(),
0376       nBinsX, hist3d.GetXaxis()->GetXmin(), hist3d.GetXaxis()->GetXmax(),
0377       nBinsY, hist3d.GetYaxis()->GetXmin(), hist3d.GetYaxis()->GetXmax());
0378
0379   auto widthHist = std::make_unique<TH2F>(
0380       widthName.c_str(), (std::string(hist3d.GetTitle()) + " width").c_str(),
0381       nBinsX, hist3d.GetXaxis()->GetXmin(), hist3d.GetXaxis()->GetXmax(),
0382       nBinsY, hist3d.GetYaxis()->GetXmin(), hist3d.GetYaxis()->GetXmax());
0383
0384   // Copy X and Y axis bin edges for variable binning
0385   if (hist3d.GetXaxis()->GetXbins()->GetSize() > 0 ||
0386       hist3d.GetYaxis()->GetXbins()->GetSize() > 0) {
0387     meanHist->SetBins(nBinsX, hist3d.GetXaxis()->GetXbins()->GetArray(), nBinsY,
0388                       hist3d.GetYaxis()->GetXbins()->GetArray());
0389     widthHist->SetBins(nBinsX, hist3d.GetXaxis()->GetXbins()->GetArray(),
0390                        nBinsY, hist3d.GetYaxis()->GetXbins()->GetArray());
0391   }
0392
0393   // Loop over all (X,Y) bins
0394   int fitFailures = 0;
0395   for (int i = 1; i <= nBinsX; ++i) {
0396     for (int j = 1; j <= nBinsY; ++j) {
0397       const auto proj = std::unique_ptr<TH1D>(hist3d.ProjectionZ(
0398           std::format("{}_projz_bin_{}_{}", hist3d.GetName(), i, j).c_str(), i,
0399           i, j, j));
0400
0401       if (proj->GetEntries() < minEntriesForFit) {
0402         continue;
0403       }
0404
0405       const TFitResultPtr r = proj->Fit("gaus", fitOption.c_str());
0406       if ((r.Get() == nullptr) || ((r->Status() % 1000) != 0)) {
0407         ++fitFailures;
0408         ACTS_DEBUG("Failed to fit Gaussian for bin "
0409                    << i << ", " << j << ": status "
0410                    << (r.Get() != nullptr ? r->Status() : -1));
0411         continue;
0412       }
0413
0414       // Fill mean
0415       meanHist->SetBinContent(i, j, r->Parameter(1));
0416       meanHist->SetBinError(i, j, r->ParError(1));
0417
0418       // Fill width (sigma)
0419       widthHist->SetBinContent(i, j, r->Parameter(2));
0420       widthHist->SetBinError(i, j, r->ParError(2));
0421     }
0422   }
0423   const double fitFailureFraction =
0424       (nBinsX * nBinsY > 0)
0425           ? static_cast<double>(fitFailures) / (nBinsX * nBinsY)
0426           : 0;
0427
0428   meanHist->GetXaxis()->SetTitle(hist3d.GetXaxis()->GetTitle());
0429   meanHist->GetYaxis()->SetTitle(hist3d.GetYaxis()->GetTitle());
0430   meanHist->GetZaxis()->SetTitle(hist3d.GetZaxis()->GetTitle());
0431
0432   widthHist->GetXaxis()->SetTitle(hist3d.GetXaxis()->GetTitle());
0433   widthHist->GetYaxis()->SetTitle(hist3d.GetYaxis()->GetTitle());
0434   widthHist->GetZaxis()->SetTitle(hist3d.GetZaxis()->GetTitle());
0435
0436   return {std::move(meanHist), std::move(widthHist), fitFailureFraction};
0437 }