diff --git a/bindings/pyroot/pythonizations/python/ROOT/_pythonization/_uhi.py b/bindings/pyroot/pythonizations/python/ROOT/_pythonization/_uhi.py index ec1deb574ef28..73fd3e3a03b04 100644 --- a/bindings/pyroot/pythonizations/python/ROOT/_pythonization/_uhi.py +++ b/bindings/pyroot/pythonizations/python/ROOT/_pythonization/_uhi.py @@ -77,13 +77,14 @@ def __call__(self, hist): return rebin_method(*self.ngroup, newname=hist.GetName()) -def _sum(hist, axis): +def _sum(hist, axis, args=None): """ - Represents a summation operation for histograms, which either computes the integral - (1D histograms) or projects the histogram along specified axes (2D and 3D histograms). + Represents a summation operation for histograms, which either computes the integral (1D histograms) + or projects the histogram along specified axes (projection is only for 2D and 3D histograms). Example: - ans = h[::ROOT.uhi.sum] # Compute the integral for a 1D histogram + ans = h[0:len:ROOT.uhi.sum] # Compute the integral for a 1D histogram excluding flow bins + ans_2 = h[::ROOT.uhi.sum, ::ROOT.uhi.sum] # Compute the integral for a 2D histogram including flow bins h_projected = h[:, ::ROOT.uhi.sum] # Project the Y axis for a 2D histogram h_projected = h[:, :, ::ROOT.uhi.sum] # Project the Z axis for a 3D histogram """ @@ -95,19 +96,28 @@ def _invalid_axis(axis, dim): if isinstance(axis, int): axis = (axis,) if dim == 1: - return hist.Integral() + return hist.Integral(*args) if axis == (0,) else _invalid_axis(axis, dim) if dim == 2: - return hist.ProjectionX() if axis == (0,) else hist.ProjectionY() if axis == (1,) else _invalid_axis(axis, dim) + if axis == (0,): + return hist.ProjectionY() + elif axis == (1,): + return hist.ProjectionX() + elif axis == (0, 1): + return hist.Integral() + else: + return _invalid_axis(axis, dim) if dim == 3: - # It is not possible from the interface to specify the options "yx", "zy", "zx" + # It is not possible from the interface to specify the options "xy", "yz", "xz" project_map = { - (0,): "yz", - (1,): "xz", - (2,): "xy", + (0,): "zy", + (1,): "zx", + (2,): "yx", (0, 1): "z", (0, 2): "y", (1, 2): "x", } + if axis == (0, 1, 2): + return hist.Integral() return hist.Project3D(project_map[axis]) if axis in project_map else _invalid_axis(axis, dim) raise NotImplementedError(f"Summing not implemented for {dim}D histograms") @@ -148,22 +158,26 @@ def _get_axis_len(self, axis, include_flow_bins=False): return _get_axis(self, axis).GetNbins() + (2 if include_flow_bins else 0) -def _process_index_for_axis(self, index, axis): +def _process_index_for_axis(self, index, axis, include_flow_bins=False, is_slice_stop=False): """Process an index for a histogram axis handling callables and index shifting.""" if callable(index): # If the index is a `loc`, `underflow`, `overflow`, or `len` - return _get_axis_len(self, axis) if index is len else index(self, axis) + return _get_axis_len(self, axis) + 1 if index is len else index(self, axis) if isinstance(index, int): # -1 index returns the last valid bin if index == -1: return _overflow(self, axis) - 1 + + if index == _overflow(self, axis): + return index + (1 if include_flow_bins else 0) + # Shift the indices by 1 to align with the UHI convention, # where 0 corresponds to the first bin, unlike ROOT where 0 represents underflow and 1 is the first bin. + nbins = _get_axis_len(self, axis) + (1 if is_slice_stop else 0) index = index + 1 - nbins = _get_axis_len(self, axis) if abs(index) > nbins: - raise IndexError(f"Histogram index {index} out of range for axis {axis}") + raise IndexError(f"Histogram index {index-1} out of range for axis {axis}. Valid range: (0,{nbins})") return index raise index @@ -200,14 +214,12 @@ def _compute_common_index(self, index, include_flow_bins=True): raise IndexError("Only one ellipsis is allowed in the index.") if any(idx is ... for idx in index): - expanded_index = [] - for idx in index: - if idx is ...: - break - expanded_index.append(idx) - # fill remaining dimensions with `slice(None)` - expanded_index.extend([slice(None)] * (dim - len(expanded_index))) - index = tuple(expanded_index) + ellipsis_pos = index.index(...) + index = ( + index[:ellipsis_pos] + + (slice(None),) * (dim - len(index) + 1) + + index[ellipsis_pos + 1:] + ) if len(index) != dim: raise IndexError(f"Expected {dim} indices, got {len(index)}") @@ -224,34 +236,40 @@ def _resolve_slice_indices(self, index, axis, include_flow_bins=True): """Resolve slice start and stop indices for a given axis""" start, stop = index.start, index.stop start = ( - _process_index_for_axis(self, start, axis) + _process_index_for_axis(self, start, axis, include_flow_bins) if start is not None else _underflow(self, axis) + (0 if include_flow_bins else 1) ) stop = ( - _process_index_for_axis(self, stop, axis) + _process_index_for_axis(self, stop, axis, include_flow_bins, is_slice_stop=True) if stop is not None else _overflow(self, axis) + (1 if include_flow_bins else 0) ) if start < _underflow(self, axis) or stop > (_overflow(self, axis) + 1) or start > stop: - raise IndexError(f"Slice indices {start, stop} out of range for axis {axis}") + raise IndexError(f"Slice indices {start, stop} out of range for axis {axis}. Valid range: {_underflow(self, axis), _overflow(self, axis) + 1}") return start, stop -def _apply_actions(hist, actions): +def _apply_actions(hist, actions, index, unprocessed_index, original_hist): """Apply rebinning or summing actions to the histogram, returns a new histogram""" if not actions or all(a is None for a in actions): return hist - - if any(a is _sum for a in actions): - sum_axes = tuple(i for i, a in enumerate(actions) if a is _sum) - hist = _sum(hist, sum_axes) + + if any(a is _sum or a is sum for a in actions): + sum_axes = tuple(i for i, a in enumerate(actions) if a is _sum or a is sum) + if original_hist.GetDimension() == 1: + start, stop = index[0].start, index[0].stop + include_oflow = True if unprocessed_index.stop is None else False + args = [start, stop - (1 if not include_oflow else 0)] + hist = _sum(original_hist, sum_axes, args) + else: + hist = _sum(hist, sum_axes) if any(isinstance(a, _rebin) for a in actions): rebins = [a.ngroup if isinstance(a, _rebin) else 1 for a in actions if a is not _sum] hist = _rebin(rebins)(hist) - if any(a is not None and not (isinstance(a, _rebin) or a is _sum) for a in actions): + if any(a is not None and not (isinstance(a, _rebin) or a is _sum or a is sum) for a in actions): raise ValueError(f"Unsupported action detected in actions {actions}") return hist @@ -261,102 +279,41 @@ def _get_processed_slices(self, index): """Process slices and extract actions for each axis""" if len(index) != self.GetDimension(): raise IndexError(f"Expected {self.GetDimension()} indices, got {len(index)}") - processed_slices, out_of_range_indices, actions = [], [], [None] * self.GetDimension() + processed_slices, actions = [], [None] * self.GetDimension() for axis, idx in enumerate(index): - axis_bins = range(_overflow(self, axis) + 1) if isinstance(idx, slice): slice_range = range(idx.start, idx.stop) processed_slices.append(slice_range) - uflow = [b for b in axis_bins if b < idx.start] - oflow = [b for b in axis_bins if b >= idx.stop] - out_of_range_indices.append((uflow, oflow)) actions[axis] = idx.step + elif isinstance(idx, int): + # A single value v is like v:v+1:sum, example: h2 = h[v, a:b] + processed_slices.append(range(idx, idx + 1)) + actions[axis] = _sum else: processed_slices.append([idx]) - return processed_slices, out_of_range_indices, actions - - -def _get_slice_indices(slices): - """ - This function uses numpy's meshgrid to create a grid of indices from the input slices, - and reshapes the grid into a list of all possible index combinations. - - Example: - slices = [range(2), range(3)] - # This represents two dimensions: - # - The first dimension has indices [0, 1] - # - The second dimension has indices [0, 1, 2] - - result = _get_slice_indices(slices) - # result: - # [[0, 0], - # [0, 1], - # [0, 2], - # [1, 0], - # [1, 1], - # [1, 2]] - """ - import numpy as np - - grids = np.meshgrid(*slices, indexing="ij") - return np.array(grids).reshape(len(slices), -1).T - - -def _set_flow_bins(self, target_hist, out_of_range_indices): - """ - Accumulate content from bins outside the slice range into flow bins. - """ - dim = self.GetDimension() - uflow_bin = tuple(_underflow(self, axis) for axis in range(dim)) - oflow_bin = tuple(_overflow(self, axis) for axis in range(dim)) - flow_sum = 0 - - for axis, (underflow_indices, overflow_indices) in enumerate(out_of_range_indices): - all_axes = [range(_overflow(self, j)) for j in range(dim)] - - def sum_bin_content(indices_list, target_bin): - current_val = target_hist.GetBinContent(*target_bin) - temp_axes = list(all_axes) - temp_axes[axis] = indices_list - for idx in _get_slice_indices(temp_axes): - current_val += self.GetBinContent(*tuple(map(int, idx))) - target_hist.SetBinContent(*target_bin, current_val) - return current_val - - flow_sum += sum_bin_content(underflow_indices, uflow_bin) - flow_sum += sum_bin_content(overflow_indices, oflow_bin) - - return flow_sum + return processed_slices, actions -def _slice_get(self, index): +def _slice_get(self, index, unprocessed_index): """ This method creates a new histogram containing only the data from the specified slice. Steps: - Process the slices and extract the actions for each axis. - - Clone the original histogram and reset its contents. - - Set the bin content for each index in the slice. - - Update the number of entries in the cloned histogram (also updates the statistics). + - Get a new sliced histogram. - Apply any rebinning or summing actions to the resulting histogram. """ - processed_slices, out_of_range_indices, actions = _get_processed_slices(self, index) - slice_indices = _get_slice_indices(processed_slices) - with _temporarily_disable_add_directory(): - target_hist = self.Clone() - target_hist.Reset() - - for indices in slice_indices: - indices = tuple(map(int, indices)) - target_hist.SetBinContent(*indices, self.GetBinContent(self.GetBin(*indices))) + import ROOT - flow_sum = _set_flow_bins(self, target_hist, out_of_range_indices) + processed_slices, actions = _get_processed_slices(self, index) + start_stop = [(r.start, r.stop) for r in processed_slices] + args_vec = ROOT.std.vector('Int_t')([item for pair in start_stop for item in pair]) - target_hist.SetEntries(target_hist.GetEffectiveEntries() + flow_sum) + target_hist = ROOT.Internal.Slice(self, args_vec) - return _apply_actions(target_hist, actions) + return _apply_actions(target_hist, actions, index, unprocessed_index, self) def _slice_set(self, index, unprocessed_index, value): @@ -367,42 +324,49 @@ def _slice_set(self, index, unprocessed_index, value): """ import numpy as np + import ROOT + + if isinstance(value, (list, range)): + value = np.array(value) + # Depending on the shape of the array provided, we can set or not the flow bins # Setting with a scalar does not set the flow bins - include_flow_bins = not ( - (isinstance(value, np.ndarray) and value.shape == _shape(self, include_flow_bins=False)) or np.isscalar(value) + include_flow_bins = ( + (isinstance(value, np.ndarray) and value.shape != _shape(_slice_get(self, index, unprocessed_index), include_flow_bins=False)) or np.isscalar(value) ) if not include_flow_bins: index = _compute_common_index(self, unprocessed_index, include_flow_bins=False) - processed_slices, _, actions = _get_processed_slices(self, index) - slice_indices = _get_slice_indices(processed_slices) - if isinstance(value, np.ndarray): - if value.size != len(slice_indices): - raise ValueError(f"Expected {len(slice_indices)} bin values, got {value.size}") - - expected_shape = tuple(len(slice_range) for slice_range in processed_slices) - if value.shape != expected_shape: - raise ValueError(f"Shape mismatch: expected {expected_shape}, got {value.shape}") - - for indices, val in zip(slice_indices, value.ravel()): - _setbin(self, self.GetBin(*map(int, indices)), val) - elif np.isscalar(value): - for indices in slice_indices: - _setbin(self, self.GetBin(*map(int, indices)), value) + processed_slices, actions = _get_processed_slices(self, index) + start_stop = [(r.start, r.stop) for r in processed_slices] + slice_shape = tuple(stop - start for start, stop in start_stop) + args_vec = ROOT.std.vector('Int_t')([item for pair in start_stop for item in pair]) + + if np.isscalar(value): + value = ROOT.std.vector('Double_t')([value] * np.prod(slice_shape)) else: - raise TypeError(f"Unsupported value type: {type(value).__name__}") - - _apply_actions(self, actions) + try: + value = np.asanyarray(value) + if value.size != np.prod(slice_shape): + try: + value = np.broadcast_to(value, slice_shape) + except ValueError: + raise ValueError(f"Expected {np.prod(slice_shape)} bin values, got {value.size}") + value = ROOT.std.vector('Double_t')(value.flatten().astype(np.float64)) + except AttributeError: + raise TypeError(f"Unsupported value type: {type(value).__name__}") + + ROOT.Internal.SetSliceContent(self, value, args_vec) + + _apply_actions(self, actions, index, unprocessed_index, self) def _getitem(self, index): uhi_index = _compute_common_index(self, index) if all(isinstance(i, int) for i in uhi_index): return self.GetBinContent(*uhi_index) - if any(isinstance(i, slice) for i in uhi_index): - return _slice_get(self, uhi_index) + return _slice_get(self, uhi_index, index) def _setitem(self, index, value): diff --git a/bindings/pyroot/pythonizations/test/uhi_indexing.py b/bindings/pyroot/pythonizations/test/uhi_indexing.py index ccadde29703e1..fa59ad33f86c9 100644 --- a/bindings/pyroot/pythonizations/test/uhi_indexing.py +++ b/bindings/pyroot/pythonizations/test/uhi_indexing.py @@ -4,12 +4,7 @@ import pytest import ROOT -from ROOT._pythonization._uhi import ( - _get_axis, - _get_processed_slices, - _get_slice_indices, - _shape, -) +from ROOT._pythonization._uhi import _get_axis, _get_processed_slices, _overflow, _shape, _underflow from ROOT.uhi import loc, overflow, rebin, sum, underflow @@ -32,6 +27,12 @@ def _iterate_bins(hist): yield tuple(filter(None, (i, j, k))) +def _get_slice_indices(slices): + import numpy as np + + grids = np.meshgrid(*slices, indexing="ij") + return np.array(grids).reshape(len(slices), -1).T + class TestTH1Indexing: def test_access_with_bin_number(self, hist_setup): for index in [0, 8]: @@ -52,7 +53,7 @@ def test_access_flow_bins(self, hist_setup): def test_access_with_len(self, hist_setup): len_indices = (len,) * hist_setup.GetDimension() - bin_counts = (_get_axis(hist_setup, i).GetNbins() for i in range(hist_setup.GetDimension())) + bin_counts = (_get_axis(hist_setup, i).GetNbins() + 1 for i in range(hist_setup.GetDimension())) assert hist_setup[len_indices] == hist_setup.GetBinContent(*bin_counts) def test_access_with_ellipsis(self, hist_setup): @@ -110,17 +111,34 @@ def test_setting_with_scalar(self, hist_setup): def _test_slices_match(self, hist_setup, slice_ranges, processed_slices): dim = hist_setup.GetDimension() - slices, _, _ = _get_processed_slices(hist_setup, processed_slices[dim]) + slices, _ = _get_processed_slices(hist_setup, processed_slices[dim]) expected_indices = _get_slice_indices(slices) sliced_hist = hist_setup[tuple(slice_ranges[dim])] for bin_indices in expected_indices: bin_indices = tuple(map(int, bin_indices)) - assert sliced_hist.GetBinContent(*bin_indices) == hist_setup.GetBinContent(*bin_indices) - - for bin_indices in _iterate_bins(hist_setup): - if list(bin_indices) not in expected_indices.tolist(): - assert sliced_hist.GetBinContent(*bin_indices) == 0 + shifted_indices = [] + is_flow_bin = False + for i, idx in enumerate(bin_indices): + shift = slice_ranges[dim][i].start + if callable(shift): + shift = shift(hist_setup, i) + elif shift is None: + shift = 1 + else: + shift += 1 + + shifted_idx = idx - shift + 1 + if shifted_idx <= 0 or shifted_idx == _overflow(hist_setup, i): + is_flow_bin = True + break + + shifted_indices.append(shifted_idx) + + if is_flow_bin: + continue + + assert sliced_hist.GetBinContent(*tuple(shifted_indices)) == hist_setup.GetBinContent(*bin_indices) def test_slicing_with_endpoints(self, hist_setup): if _special_setting(hist_setup): @@ -144,13 +162,13 @@ def test_slicing_without_endpoints(self, hist_setup): processed_slices = { 1: [slice(0, 8)], - 2: [slice(0, 8), slice(4, 11)], - 3: [slice(0, 8), slice(4, 11), slice(3, 6)], + 2: [slice(0, 8), slice(0, 8)], + 3: [slice(0, 8), slice(0, 8), slice(3, 6)], } slice_ranges = { 1: [slice(None, 7)], - 2: [slice(None, 7), slice(3, None)], - 3: [slice(None, 7), slice(3, None), slice(2, 5)], + 2: [slice(None, 7), slice(None, 7)], + 3: [slice(None, 7), slice(None, 7), slice(2, 5)], } self._test_slices_match(hist_setup, slice_ranges, processed_slices) @@ -160,17 +178,17 @@ def test_slicing_with_data_coordinates(self, hist_setup): processed_slices = { 1: [slice(hist_setup.FindBin(2), 11)], - 2: [slice(hist_setup.FindBin(2), 11), slice(hist_setup.FindBin(3), 11)], + 2: [slice(hist_setup.FindBin(2)-1, 11), slice(2, 11)], 3: [ slice(hist_setup.FindBin(2), 11), - slice(hist_setup.FindBin(3), 11), - slice(hist_setup.FindBin(1.5), 11), + slice(2, 11), + slice(2, 11), ], } slice_ranges = { 1: [slice(loc(2), None)], - 2: [slice(loc(2), None), slice(loc(3), None)], - 3: [slice(loc(2), None), slice(loc(3), None), slice(loc(1.5), None)], + 2: [slice(loc(2), None), slice(3, None)], + 3: [slice(loc(2), None), slice(3, None), slice(3, None)], } self._test_slices_match(hist_setup, slice_ranges, processed_slices) @@ -191,7 +209,10 @@ def test_slicing_over_everything_with_action_sum(self, hist_setup): dim = hist_setup.GetDimension() if dim == 1: - integral = hist_setup[::sum] + full_integral = hist_setup[::sum] + assert full_integral == hist_setup.Integral(_underflow(hist_setup, 0), _overflow(hist_setup, 0)) + + integral = hist_setup[0:len:sum] assert integral == hist_setup.Integral() if dim == 2: @@ -225,7 +246,7 @@ def test_slicing_with_action_rebin_and_sum(self, hist_setup): if dim == 1: sliced_hist_rebin = hist_setup[5 : 9 : rebin(2)] assert isinstance(sliced_hist_rebin, ROOT.TH1) - assert sliced_hist_rebin.GetNbinsX() == hist_setup.GetNbinsX() // 2 + assert sliced_hist_rebin.GetNbinsX() == 2 sliced_hist_sum = hist_setup[5:9:sum] assert isinstance(sliced_hist_sum, float) @@ -237,10 +258,10 @@ def test_slicing_with_action_rebin_and_sum(self, hist_setup): assert sliced_hist.GetNbinsX() == hist_setup.GetNbinsX() // 2 if dim == 3: - sliced_hist = hist_setup[:: rebin(2), ::sum, 5 : 9 : rebin(3)] + sliced_hist = hist_setup[:: rebin(2), ::sum, 3 : 9 : rebin(3)] assert isinstance(sliced_hist, ROOT.TH2) assert sliced_hist.GetNbinsX() == hist_setup.GetNbinsX() // 2 - assert sliced_hist.GetNbinsY() == hist_setup.GetNbinsZ() // 3 + assert sliced_hist.GetNbinsY() == 2 def test_slicing_with_dict_syntax(self, hist_setup): if _special_setting(hist_setup): @@ -262,16 +283,22 @@ def test_integral_full_slice(self, hist_setup): assert hist_setup.Integral() == pytest.approx(sliced_hist.Integral(), rel=10e-6) def test_statistics_slice(self, hist_setup): - if _special_setting(hist_setup): + if _special_setting(hist_setup) or isinstance(hist_setup, (ROOT.TH1C, ROOT.TH2C, ROOT.TH3C)): pytest.skip("Setting cannot be tested here") + # Check if slicing over everything preserves the statistics + sliced_hist_full = hist_setup[...] + + assert hist_setup.GetEffectiveEntries() == sliced_hist_full.GetEffectiveEntries() + assert hist_setup.Integral() == sliced_hist_full.Integral() + # Check if slicing over a range updates the statistics dim = hist_setup.GetDimension() - [_get_axis(hist_setup, i).SetRange(3, 5) for i in range(dim)] + [_get_axis(hist_setup, i).SetRange(3, 7) for i in range(dim)] slice_indices = tuple(slice(2, 7) for _ in range(dim)) sliced_hist = hist_setup[slice_indices] - assert hist_setup.Integral() == pytest.approx(sliced_hist.Integral(), rel=1e-6) + assert hist_setup.Integral() == sliced_hist.Integral() assert hist_setup.GetMean() == pytest.approx(sliced_hist.GetMean(), abs=1e-3) assert hist_setup.GetStdDev() == pytest.approx(sliced_hist.GetStdDev(), abs=1e-3) assert hist_setup.GetEffectiveEntries() == sliced_hist.GetEffectiveEntries() diff --git a/hist/hist/inc/TH1.h b/hist/hist/inc/TH1.h index afdddd087e055..727e5f546e30e 100644 --- a/hist/hist/inc/TH1.h +++ b/hist/hist/inc/TH1.h @@ -44,6 +44,7 @@ #include #include +#include class TF1; class TH1D; @@ -56,6 +57,88 @@ class TVirtualHistPainter; class TRandom; +namespace{ + /* + merci https://stackoverflow.com/a/73551511 + */ + template + struct unpack_caller { + private: + template + void call(FuncType &f, const std::vector &args, std::index_sequence) { + f(args[I]...); + } + + public: + template + void operator()(FuncType &f, const std::vector &args) { + call(f, args, std::make_index_sequence{}); + } + }; + + std::vector> getSliceIndices(const std::vector>& sliceEdges) { + const auto ndim = sliceEdges.size(); + if (ndim == 0) return {}; + + std::vector> slices(ndim); + for (size_t i = 0; i < ndim; ++i) { + for (int val = sliceEdges[i].first; val < sliceEdges[i].second; ++val) { + slices[i].push_back(val); + } + } + + size_t totalCombinations = 1; + for (const auto& slice : slices) { + totalCombinations *= slice.size(); + } + + std::vector> result(totalCombinations, std::vector(slices.size())); + for (size_t dim = 0; dim < slices.size(); ++dim) { + size_t repeat = 1; + for (size_t i = dim + 1; i < slices.size(); ++i) { + repeat *= slices[i].size(); + } + + size_t index = 0; + for (size_t i = 0; i < totalCombinations; ++i) { + result[i][dim] = slices[dim][(index / repeat) % slices[dim].size()]; + ++index; + } + } + + return result; + } +} + +namespace ROOT::Internal { + template + auto Slice(const T &histo, std::vector& args) + { + T slicedHisto = histo; + slicedHisto.SliceHistoInPlace(args); + return slicedHisto; + } + + template + auto SetSliceContent(T &histo, const std::vector &input, const std::vector& args) + { + switch (args.size() / 2) { + case 1: + histo.template SetSliceContent<1>(input, args); + break; + case 2: + histo.template SetSliceContent<2>(input, args); + break; + case 3: + histo.template SetSliceContent<3>(input, args); + break; + default: + throw std::invalid_argument("Invalid number of arguments for SetSliceContent"); + } + } +} // namespace ROOT::Internal + + class TH1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker { public: @@ -137,6 +220,47 @@ class TH1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker { TH1(const TH1&) = delete; TH1& operator=(const TH1&) = delete; + void SliceHistoInPlace(std::vector& args); + + template + friend auto ROOT::Internal::Slice(const T &histo, std::vector& args); + + template + void SetSliceContent(const std::vector &input, const std::vector& args) { + if (Dim != fDimension) { + throw std::invalid_argument("Number of edges in the specified slice does not match histogram dimension."); + } + + // Extract low/up edges in pairs + std::vector> sliceEdges; + for (size_t i = 0; i < Dim; ++i) { + sliceEdges.emplace_back(args[i * 2], args[i * 2 + 1]); + } + + // Get the indices to set + auto sliceIndices = getSliceIndices(sliceEdges); + + if (input.size() != sliceIndices.size()) { + throw std::invalid_argument("Number of provided values does not match number of bins to set."); + } + + unpack_caller caller; + for (size_t i = 0; i < sliceIndices.size(); ++i) { + auto globalBin = 0; + // Compute the bin to set + auto getGlobalBin = [&](auto... idx) { + globalBin = this->GetBin(idx...); + }; + caller(getGlobalBin, sliceIndices[i]); + + // Set the bin content + SetBinContent(globalBin, input[i]); + } + } + + template + friend auto ROOT::Internal::SetSliceContent(T &histo, const std::vector &input, const std::vector& args); + protected: TH1(); TH1(const char *name,const char *title,Int_t nbinsx,Double_t xlow,Double_t xup); diff --git a/hist/hist/src/TH1.cxx b/hist/hist/src/TH1.cxx index d252261cc3bd2..c6e85cc78f265 100644 --- a/hist/hist/src/TH1.cxx +++ b/hist/hist/src/TH1.cxx @@ -22,6 +22,7 @@ #include #include #include +#include #include "TROOT.h" #include "TBuffer.h" @@ -9430,6 +9431,90 @@ TH1* TH1::TransformHisto(TVirtualFFT *fft, TH1* h_output, Option_t *option) return hout; } +//////////////////////////////////////////////////////////////////////////////// +/// TODO docstring + +void TH1::SliceHistoInPlace(std::vector& args){ + auto originalHisto = static_cast(this->Clone()); + this->Reset(); + const auto ndim = args.size() / 2; + if (ndim != static_cast(fDimension)) { + throw std::invalid_argument("Number of dimensions in slice does not match histogram dimension."); + } + + constexpr Int_t kMaxDim = 3; + assert(ndim <= kMaxDim && "Number of dimensions exceeds the maximum allowed dimensions"); + + Int_t nBins[kMaxDim] = {0}; + Double_t binLowEdge[kMaxDim] = {0.0}, binUpEdge[kMaxDim] = {0.0}; + std::vector shiftedBins(ndim*2); + + // Configure all axes + for (size_t i = 0; i < ndim; ++i) { + const auto &axis = (i == 0) ? originalHisto->fXaxis : (i == 1) ? originalHisto->fYaxis : originalHisto->fZaxis; + args[i * 2] = std::max(1, args[i * 2]); + args[i * 2 + 1] = std::min(axis.GetNbins() + 1, args[i * 2 + 1]); + nBins[i] = args[i * 2 + 1] - args[i * 2]; + binLowEdge[i] = axis.GetBinLowEdge(args[i * 2]); + binUpEdge[i] = axis.GetBinLowEdge(args[i * 2 + 1]); + shiftedBins[i * 2] = 1; + shiftedBins[i * 2 + 1] = nBins[i] + 1; + } + + // Configure name, title and bins + this->SetNameTitle(Form("%s_slice", originalHisto->GetName()), originalHisto->GetTitle()); + + if (ndim == 1) { + this->SetBins(nBins[0], binLowEdge[0], binUpEdge[0]); + } else if (ndim == 2) { + this->SetBins(nBins[0], binLowEdge[0], binUpEdge[0], + nBins[1], binLowEdge[1], binUpEdge[1]); + } else if (ndim == 3) { + this->SetBins(nBins[0], binLowEdge[0], binUpEdge[0], + nBins[1], binLowEdge[1], binUpEdge[1], + nBins[2], binLowEdge[2], binUpEdge[2]); + } + + std::vector sliceContents; + std::unordered_set processedBins; + Double_t underflowContent = 0.0, overflowContent = 0.0; + + auto processBinContent = [&](Int_t x, Int_t y, Int_t z) { + const Int_t globalBin = originalHisto->GetBin(x, y, z); + if (!processedBins.insert(globalBin).second) return; + + const auto content = originalHisto->RetrieveBinContent(globalBin); + const bool isUnderflow = (x < args[0] || (ndim > 1 && y < args[2]) || (ndim > 2 && z < args[4])); + const bool isOverflow = (x >= args[1] || (ndim > 1 && y >= args[3]) || (ndim > 2 && z >= args[5])); + + if (isUnderflow) { + underflowContent += content; + } else if (isOverflow) { + overflowContent += content; + } else { + sliceContents.push_back(content); + } + }; + + for (Int_t x = 0; x <= originalHisto->fXaxis.GetNbins() + 1; ++x) { + for (Int_t y = 0; y <= (ndim > 1 ? originalHisto->fYaxis.GetNbins() + 1 : 1); ++y) { + for (Int_t z = 0; z <= (ndim > 2 ? originalHisto->fZaxis.GetNbins() + 1 : 1); ++z) { + processBinContent(x, y, z); + } + } + } + + // Set the contents of the slice histogram + ROOT::Internal::SetSliceContent(*this, sliceContents, shiftedBins); + + // Set the flow bins + this->SetBinContent(0, underflowContent); + auto overflowBin = (ndim == 1) ? this->GetBin(fXaxis.GetNbins() + 1) + : (ndim == 2) ? this->GetBin(fXaxis.GetNbins() + 1, fYaxis.GetNbins() + 1) + : this->GetBin(fXaxis.GetNbins() + 1, fYaxis.GetNbins() + 1, fZaxis.GetNbins() + 1); + this->SetBinContent(overflowBin, overflowContent); +} + //////////////////////////////////////////////////////////////////////////////// /// Print value overload