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 <cmath>
0012 #include <vector>
0013
0014 #include <TEfficiency.h>
0015 #include <TH1F.h>
0016 #include <TH2F.h>
0017 #include <TProfile.h>
0018 #include <boost/histogram/accumulators/weighted_sum.hpp>
0019
0020 using namespace Acts::Experimental;
0021
0022 namespace {
0023
0024 std::vector<double> extractBinEdges(const AxisVariant& axis) {
0025   assert(axis.size() > 0 && "Axis must have at least one bin");
0026   std::vector<double> edges(axis.size() + 1);
0027   for (int i = 0; i < axis.size(); ++i) {
0028     edges.at(i) = axis.bin(i).lower();
0029   }
0030   edges.back() = axis.bin(axis.size() - 1).upper();
0031
0032   return edges;
0033 }
0034
0035 }  // namespace
0036
0037 namespace ActsPlugins {
0038
0039 TH1F* toRoot(const Histogram1& boostHist) {
0040   const auto& bh = boostHist.histogram();
0041   const auto& axis = bh.axis(0);
0042
0043   // Extract bin edges from boost histogram axis
0044   std::vector<double> edges = extractBinEdges(axis);
0045
0046   // Create ROOT histogram with variable binning
0047   TH1F* rootHist = new TH1F(boostHist.name().c_str(), boostHist.title().c_str(),
0048                             static_cast<int>(axis.size()), edges.data());
0049
0050   // Copy bin contents from boost to ROOT
0051   for (auto&& x : boost::histogram::indexed(bh)) {
0052     // Dereference to get bin content
0053     double content = *x;
0054
0055     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0056     int rootBinIndex = static_cast<int>(x.index(0)) + 1;
0057     rootHist->SetBinContent(rootBinIndex, content);
0058   }
0059
0060   // Set axis titles from axis metadata
0061   rootHist->GetXaxis()->SetTitle(axis.metadata().c_str());
0062
0063   return rootHist;
0064 }
0065
0066 TH2F* toRoot(const Histogram2& boostHist) {
0067   const auto& bh = boostHist.histogram();
0068   const auto& xAxis = bh.axis(0);
0069   const auto& yAxis = bh.axis(1);
0070
0071   // Extract bin edges from X axis
0072   std::vector<double> xEdges = extractBinEdges(xAxis);
0073
0074   // Extract bin edges from Y axis
0075   std::vector<double> yEdges = extractBinEdges(yAxis);
0076
0077   // Create ROOT histogram with 2D variable binning
0078   TH2F* rootHist = new TH2F(boostHist.name().c_str(), boostHist.title().c_str(),
0079                             static_cast<int>(xAxis.size()), xEdges.data(),
0080                             static_cast<int>(yAxis.size()), yEdges.data());
0081
0082   // Copy bin contents from boost to ROOT
0083   for (auto&& x : boost::histogram::indexed(bh)) {
0084     // Dereference to get bin content
0085     double content = *x;
0086
0087     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0088     // indexed() gives us 0-based bin indices for each axis
0089     int rootXBin = static_cast<int>(x.index(0)) + 1;
0090     int rootYBin = static_cast<int>(x.index(1)) + 1;
0091     rootHist->SetBinContent(rootXBin, rootYBin, content);
0092   }
0093
0094   // Set axis titles from axis metadata
0095   rootHist->GetXaxis()->SetTitle(xAxis.metadata().c_str());
0096   rootHist->GetYaxis()->SetTitle(yAxis.metadata().c_str());
0097
0098   return rootHist;
0099 }
0100
0101 TProfile* toRoot(const ProfileHistogram1& boostProfile) {
0102   const auto& bh = boostProfile.histogram();
0103   const auto& axis = bh.axis(0);
0104
0105   // Extract bin edges from boost histogram axis
0106   std::vector<double> edges = extractBinEdges(axis);
0107
0108   // Create ROOT TProfile with variable binning
0109   TProfile* rootProfile =
0110       new TProfile(boostProfile.name().c_str(), boostProfile.title().c_str(),
0111                    static_cast<int>(axis.size()), edges.data());
0112
0113   // Enable sum of weights squared storage for proper error calculation
0114   rootProfile->Sumw2();
0115
0116   // Copy data from boost profile to ROOT profile
0117   // - count(): number of fills
0118   // - value(): mean of y-values
0119   // - variance(): sample variance = sum((y - mean)^2) / (n - 1)
0120
0121   using Accumulator = boost::histogram::accumulators::mean<double>;
0122
0123   for (auto&& x : boost::histogram::indexed(bh)) {
0124     const Accumulator& acc = *x;
0125
0126     // ROOT bin numbering starts at 1 (bin 0 is underflow)
0127     int rootBinIndex = static_cast<int>(x.index(0)) + 1;
0128
0129     double count = acc.count();
0130
0131     if (count == 0) {
0132       continue;
0133     }
0134     double mean = acc.value();
0135     double variance = acc.variance();  // Sample variance from boost
0136
0137     // ROOT TProfile internally stores the following data (all weights=1):
0138     // - fArray[bin]      = sum of (w*x)    = count * mean
0139     // - fBinEntries[bin] = sum of w        = count
0140     // - fSumw2[bin]      = sum of (w*x)^2  = sum of x^2
0141     // - fBinSumw2[bin]   = sum of w^2      = count
0142
0143     double sum = mean * count;
0144
0145     // To compute sum of y^2 from sample variance s^2:
0146     // s^2 = Sum((x - μ)^2) / (n-1) = (Sum(x^2) - n*µ^2) / (n-1)
0147     // Sum(x^2) = (n-1) * s^2 + n*µ^2
0148     double sumOfSquares = (count - 1.0) * variance + count * mean * mean;
0149
0150     // Set bin content (sum) and entries via public interface
0151     rootProfile->SetBinContent(rootBinIndex, sum);
0152     rootProfile->SetBinEntries(rootBinIndex, count);
0153
0154     // Access internal arrays to set sum of squares for error calculation
0155     // This is necessary because ROOT has no public API to set these arrays
0156     TArrayD* sumw2 = rootProfile->GetSumw2();        // fSumw2 array
0157     TArrayD* binSumw2 = rootProfile->GetBinSumw2();  // fBinSumw2 array
0158
0159     assert(sumw2 && "Sumw2 is null");
0160     assert(binSumw2 && "BinSumw2 is null");
0161
0162     // Set sum of (weight * value)^2 = sum of y^2 for unweighted data
0163     sumw2->fArray[rootBinIndex] = sumOfSquares;
0164     // Set sum of weights^2 = count for unweighted data (all weights = 1)
0165     binSumw2->fArray[rootBinIndex] = count;
0166   }
0167
0168   // Set X axis title from axis metadata
0169   rootProfile->GetXaxis()->SetTitle(axis.metadata().c_str());
0170   // Set Y axis title from sampleAxisTitle member
0171   rootProfile->GetYaxis()->SetTitle(boostProfile.sampleAxisTitle().c_str());
0172
0173   return rootProfile;
0174 }
0175
0176 TEfficiency* toRoot(const Efficiency1& boostEff) {
0177   const auto& accepted = boostEff.acceptedHistogram();
0178   const auto& total = boostEff.totalHistogram();
0179   const auto& axis = accepted.axis(0);
0180
0181   // Extract bin edges
0182   std::vector<double> edges = extractBinEdges(axis);
0183
0184   // Create accepted and total TH1F histograms
0185   TH1F* acceptedHist =
0186       new TH1F((boostEff.name() + "_accepted").c_str(), "Passed",
0187                static_cast<int>(axis.size()), edges.data());
0188
0189   TH1F* totalHist = new TH1F((boostEff.name() + "_total").c_str(), "Total",
0190                              static_cast<int>(axis.size()), edges.data());
0191
0192   // Fill histograms with counts
0193   for (int i = 0; i < axis.size(); ++i) {
0194     double acceptedCount = static_cast<double>(accepted.at(i));
0195     double totalCount = static_cast<double>(total.at(i));
0196
0197     acceptedHist->SetBinContent(i + 1, acceptedCount);
0198     totalHist->SetBinContent(i + 1, totalCount);
0199   }
0200
0201   // Create TEfficiency from the two histograms
0202   // TEfficiency takes ownership of the histograms
0203   TEfficiency* rootEff = new TEfficiency(*acceptedHist, *totalHist);
0204   rootEff->SetName(boostEff.name().c_str());
0205   rootEff->SetTitle(boostEff.title().c_str());
0206
0207   // Clean up temporary histograms (TEfficiency made copies)
0208   delete acceptedHist;
0209   delete totalHist;
0210
0211   return rootEff;
0212 }
0213
0214 TEfficiency* toRoot(const Efficiency2& boostEff) {
0215   const auto& accepted = boostEff.acceptedHistogram();
0216   const auto& total = boostEff.totalHistogram();
0217   const auto& xAxis = accepted.axis(0);
0218   const auto& yAxis = accepted.axis(1);
0219
0220   // Extract X bin edges
0221   std::vector<double> xEdges = extractBinEdges(xAxis);
0222
0223   // Extract Y bin edges
0224   std::vector<double> yEdges = extractBinEdges(yAxis);
0225
0226   // Create accepted and total TH2F histograms
0227   TH2F* acceptedHist =
0228       new TH2F((boostEff.name() + "_accepted").c_str(), "Accepted",
0229                static_cast<int>(xAxis.size()), xEdges.data(),
0230                static_cast<int>(yAxis.size()), yEdges.data());
0231
0232   TH2F* totalHist = new TH2F((boostEff.name() + "_total").c_str(), "Total",
0233                              static_cast<int>(xAxis.size()), xEdges.data(),
0234                              static_cast<int>(yAxis.size()), yEdges.data());
0235
0236   // Fill histograms with counts
0237   for (int i = 0; i < xAxis.size(); ++i) {
0238     for (int j = 0; j < yAxis.size(); ++j) {
0239       double acceptedCount = static_cast<double>(accepted.at(i, j));
0240       double totalCount = static_cast<double>(total.at(i, j));
0241
0242       acceptedHist->SetBinContent(i + 1, j + 1, acceptedCount);
0243       totalHist->SetBinContent(i + 1, j + 1, totalCount);
0244     }
0245   }
0246
0247   // Create TEfficiency from the two histograms
0248   // TEfficiency takes ownership of the histograms
0249   TEfficiency* rootEff = new TEfficiency(*acceptedHist, *totalHist);
0250   rootEff->SetName(boostEff.name().c_str());
0251   rootEff->SetTitle(boostEff.title().c_str());
0252
0253   // Clean up temporary histograms (TEfficiency made copies)
0254   delete acceptedHist;
0255   delete totalHist;
0256
0257   return rootEff;
0258 }
0259
0260 }  // namespace ActsPlugins