diff --git a/hist/hist/inc/TH2.h b/hist/hist/inc/TH2.h index 40ef814dde7fb..c783048aaeebf 100644 --- a/hist/hist/inc/TH2.h +++ b/hist/hist/inc/TH2.h @@ -115,7 +115,7 @@ class TH2 : public TH1 { TH2 *RebinX(Int_t ngroup=2, const char *newname="") override; // *MENU* virtual TH2 *RebinY(Int_t ngroup=2, const char *newname=""); // *MENU* TH2 *Rebin(Int_t ngroup=2, const char*newname="", const Double_t *xbins = nullptr) override; // re-implementation of the TH1 function using RebinX - virtual TH2 *Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname=""); // *MENU* + virtual TH2 *Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname="", const Double_t *xbins=nullptr, const Double_t *ybins=nullptr); // *MENU* TProfile *ProfileX(const char *name="_pfx", Int_t firstybin=1, Int_t lastybin=-1, Option_t *option="") const; // *MENU* TProfile *ProfileY(const char *name="_pfy", Int_t firstxbin=1, Int_t lastxbin=-1, Option_t *option="") const; // *MENU* TH1D *ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const; // *MENU* diff --git a/hist/hist/inc/TProfile2D.h b/hist/hist/inc/TProfile2D.h index 564141207c4f3..c7c6e8778e2b6 100644 --- a/hist/hist/inc/TProfile2D.h +++ b/hist/hist/inc/TProfile2D.h @@ -137,7 +137,7 @@ class TProfile2D : public TH2D { TProfile *ProfileY(const char *name="_pfy", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const; // *MENU* void PutStats(Double_t *stats) override; void Reset(Option_t *option="") override; - TProfile2D *Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname="") override; + TProfile2D *Rebin2D(Int_t nxgroup=2, Int_t nygroup=2, const char *newname="", const Double_t *xbins=nullptr, const Double_t *ybins=nullptr) override; TProfile2D *RebinX(Int_t ngroup=2, const char *newname="") override; TProfile2D *RebinY(Int_t ngroup=2, const char *newname="") override; void SavePrimitive(std::ostream &out, Option_t *option = "") override; diff --git a/hist/hist/src/Rebin2DHelpers.h b/hist/hist/src/Rebin2DHelpers.h new file mode 100644 index 0000000000000..83c6dbce018ea --- /dev/null +++ b/hist/hist/src/Rebin2DHelpers.h @@ -0,0 +1,184 @@ +// @(#)root/hist:$Id$ +// Author: Jonas Rembser, CERN 09/2026 + +/************************************************************************* + * Copyright (C) 1995-2026, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef ROOT_Rebin2DHelpers +#define ROOT_Rebin2DHelpers + +// Internal helpers shared by TH2::Rebin2D and TProfile2D::Rebin2D. + +#include "TAxis.h" +#include "TH1.h" +#include "TMath.h" + +#include +#include +#include + +namespace ROOT { +namespace Internal { + +/// Define the axis of the rebinned histogram: either from the user-provided +/// bin edges, or by merging groups of ngroup bins of the old axis. The value +/// of xmax is the (possibly truncated) upper limit for the uniform-bin case. +inline void DefineRebinnedAxis(const TAxis &oldAxis, Int_t ngroup, Int_t nnew, const Double_t *userBins, Double_t xmin, + Double_t xmax, TAxis &newAxis) +{ + if (userBins) { + newAxis.Set(nnew, userBins); + } else if (oldAxis.GetXbins()->GetSize() > 0) { + std::vector edges(nnew + 1); + for (Int_t i = 0; i <= nnew; ++i) + edges[i] = oldAxis.GetBinLowEdge(1 + i * ngroup); + newAxis.Set(nnew, edges.data()); + } else { + newAxis.Set(nnew, xmin, xmax); + } +} + +/// Map each cell of the old axis (including underflow 0 and overflow n+1) to +/// the cell of the new axis that contains its bin center. Old bins outside +/// the new axis range are mapped to the new under-/overflow. +inline std::vector MakeRebinMap(const TAxis &oldAxis, const TAxis &newAxis) +{ + const Int_t nOld = oldAxis.GetNbins(); + std::vector map(nOld + 2); + map[0] = 0; + map[nOld + 1] = newAxis.GetNbins() + 1; + for (Int_t o = 1; o <= nOld; ++o) + map[o] = newAxis.FindFixBin(oldAxis.GetBinCenter(o)); + return map; +} + +/// Warn when a bin edge of the new axis that lies inside the old axis range +/// does not line up with a bin edge of the old axis: the entries of the old +/// bin that is split cannot be distributed correctly. +inline void WarnAboutMisalignedEdges(const TAxis &oldAxis, const TAxis &newAxis, TH1 &hist, const char *where) +{ + for (Int_t b = 0; b <= newAxis.GetNbins(); ++b) { + const Double_t edge = newAxis.GetBinUpEdge(b); // GetBinUpEdge(0) is the axis minimum + if (edge <= oldAxis.GetXmin() || edge >= oldAxis.GetXmax()) + continue; + const Int_t o = oldAxis.FindFixBin(edge); + const Double_t tol = TMath::Max(1.E-8 * oldAxis.GetBinWidth(o), 1.E-16); + if (!TMath::AreEqualAbs(edge, oldAxis.GetBinLowEdge(o), tol) && + !TMath::AreEqualAbs(edge, oldAxis.GetBinUpEdge(o), tol)) { + hist.Warning(where, + "Bin edge %d of rebinned histogram does not match any bin edges of the old histogram. " + "Result can be inconsistent", + b); + } + } +} + +/// The definition of one axis of the rebinned histogram. +struct RebinnedAxisInfo { + Int_t nNewBins = 0; ///< number of bins of the rebinned axis + TAxis newAxis; ///< the rebinned axis + std::vector binMap; ///< map from old cell (0..n+1) to new cell + bool truncated = false; ///< the group count does not divide the old bin count: top bins move to the overflow +}; + +/// Validate the rebinning parameters for one axis and fill the definition of +/// the rebinned axis and the map from old to new bins. For an axis with +/// user-provided bin edges, ngroup is directly the new number of bins, +/// otherwise the old bins are merged in groups of ngroup. Returns false on an +/// invalid group count. +inline bool SetupRebinnedAxis(const TAxis &oldAxis, Int_t ngroup, const Double_t *userBins, char axisName, TH1 &hist, + const char *where, RebinnedAxisInfo &info) +{ + const Int_t nOldBins = oldAxis.GetNbins(); + if (ngroup <= 0 || ngroup > nOldBins) { + hist.Error(where, "Illegal value of n%cgroup=%d", axisName, ngroup); + return false; + } + Double_t newMax = oldAxis.GetXmax(); + if (userBins) { + info.nNewBins = ngroup; + } else { + info.nNewBins = nOldBins / ngroup; + if (info.nNewBins * ngroup != nOldBins) { + hist.Warning(where, "n%cgroup=%d is not an exact divider of n%cbins=%d.", axisName, ngroup, axisName, + nOldBins); + // the top limit is truncated and the top bins move to the overflow + newMax = oldAxis.GetBinUpEdge(info.nNewBins * ngroup); + info.truncated = true; + } + } + DefineRebinnedAxis(oldAxis, ngroup, info.nNewBins, userBins, oldAxis.GetXmin(), newMax, info.newAxis); + if (userBins) + WarnAboutMisalignedEdges(oldAxis, info.newAxis, hist, where); + info.binMap = MakeRebinMap(oldAxis, info.newAxis); + return true; +} + +/// Warn when the range of the new axis extends beyond the old one while the +/// corresponding flow bins hold content: that content stays in the flow bins +/// and is not redistributed into the range of the new axis. The flow bins of +/// the mapped axis are addressed as flowIndex * stride + k * otherStride for +/// the nOther cells of the other axis. +inline void WarnAboutUnusedFlowContent(const TAxis &oldAxis, const RebinnedAxisInfo &info, const Double_t *userBins, + char axisName, const Double_t *bins, Int_t stride, Int_t nOther, + Int_t otherStride, TH1 &hist, const char *where) +{ + if (!userBins) + return; + auto flowContent = [&](Int_t flowIndex) { + Double_t sum = 0.; + for (Int_t k = 0; k < nOther; ++k) + sum += bins[flowIndex * stride + k * otherStride]; + return sum; + }; + if (userBins[0] < oldAxis.GetXmin() && flowContent(0) != 0.) + hist.Warning(where, "underflow entries for %c axis will not be used when rebinning", axisName); + if (userBins[info.nNewBins] > oldAxis.GetXmax() && flowContent(oldAxis.GetNbins() + 1) != 0.) + hist.Warning(where, "overflow entries for %c axis will not be used when rebinning", axisName); +} + +/// Apply the axes of the rebinned histogram, using explicit bin edges if any +/// of the two axes has non-uniform bins. +inline void SetRebinnedBins2D(TH1 &hnew, const TAxis &newXaxis, const TAxis &newYaxis) +{ + const Int_t nx = newXaxis.GetNbins(); + const Int_t ny = newYaxis.GetNbins(); + if (newXaxis.GetXbins()->GetSize() > 0 || newYaxis.GetXbins()->GetSize() > 0) { + std::vector xEdges(nx + 1); + std::vector yEdges(ny + 1); + for (Int_t i = 0; i <= nx; ++i) + xEdges[i] = newXaxis.GetBinUpEdge(i); + for (Int_t i = 0; i <= ny; ++i) + yEdges[i] = newYaxis.GetBinUpEdge(i); + hnew.SetBins(nx, xEdges.data(), ny, yEdges.data()); // changes also errors array (if any) + } else { + hnew.SetBins(nx, newXaxis.GetXmin(), newXaxis.GetXmax(), ny, newYaxis.GetXmin(), newYaxis.GetXmax()); + } +} + +/// Accumulate every old cell (including under- and overflow) into the new +/// cell given by the per-axis bin maps, for each (old array, new array) pair. +/// The new arrays must be zero-initialized by the caller. +inline void MergeRebinnedCells(Int_t nOldX, Int_t nOldY, Int_t nNewX, const std::vector &mapX, + const std::vector &mapY, + std::initializer_list> arrays) +{ + for (Int_t oy = 0; oy < nOldY + 2; ++oy) { + for (Int_t ox = 0; ox < nOldX + 2; ++ox) { + const Int_t oldBin = ox + (nOldX + 2) * oy; + const Int_t newBin = mapX[ox] + (nNewX + 2) * mapY[oy]; + for (auto const &arr : arrays) + arr.second[newBin] += arr.first[oldBin]; + } + } +} + +} // namespace Internal +} // namespace ROOT + +#endif diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx index dfd5df93f51be..04e7ded3b810b 100644 --- a/hist/hist/src/TH2.cxx +++ b/hist/hist/src/TH2.cxx @@ -28,6 +28,9 @@ #include "TVirtualHistPainter.h" #include "snprintf.h" +#include "Rebin2DHelpers.h" + +#include /** \addtogroup Histograms @{ @@ -1627,26 +1630,27 @@ TH2 *TH2::RebinY(Int_t ngroup, const char *newname) } //////////////////////////////////////////////////////////////////////////////// -/// Override TH1::Rebin as TH2::RebinX -/// Rebinning in variable binning as for TH1 is not allowed -/// If a non-null pointer is given an error is flagged +/// Override TH1::Rebin, rebinning only the X axis with the same conventions +/// as the TH1 function (`ngroup` is the number of variable size bins when +/// `xbins` is given). /// see RebinX and Rebin2D -TH2 * TH2::Rebin( Int_t ngroup, const char*newname, const Double_t *xbins) +TH2 *TH2::Rebin(Int_t ngroup, const char *newname, const Double_t *xbins) { - if (xbins != nullptr) { - Error("Rebin","Rebinning a 2-d histogram into variable bins is not supported (it is possible only for 1-d histograms). Return a nullptr"); - return nullptr; - } + if (xbins != nullptr) + return Rebin2D(ngroup, 1, newname, xbins, nullptr); Info("Rebin","Rebinning only the x-axis. Use Rebin2D for rebinning both axes"); return RebinX(ngroup, newname); } + //////////////////////////////////////////////////////////////////////////////// /// Rebin this histogram grouping nxgroup/nygroup bins along the xaxis/yaxis together. /// -/// if newname is not blank a new temporary histogram hnew is created. +/// #### case 1 `xbins`=0 || `ybins`=0 +/// +/// if `newname` is not blank a new temporary histogram hnew is created. /// else the current histogram is modified (default) -/// The parameter nxgroup/nygroup indicate how many bins along the xaxis/yaxis of this +/// The parameters `nxgroup`/`nygroup` indicate how many bins along the xaxis/yaxis of this /// have to me merged into one bin of hnew /// If the original histogram has errors stored (via Sumw2), the resulting /// histograms has new errors correctly calculated. @@ -1660,69 +1664,87 @@ TH2 * TH2::Rebin( Int_t ngroup, const char*newname, const Double_t *xbins) /// // merging 5 bins of h1 along the yaxis in one bin /// ~~~ /// -/// NOTE : If nxgroup/nygroup is not an exact divider of the number of bins, +/// \note If `nxgroup`/`nygroup` is not an exact divider of the number of bins, /// along the xaxis/yaxis the top limit(s) of the rebinned histogram /// is changed to the upper edge of the xbin=newxbins*nxgroup resp. /// ybin=newybins*nygroup and the corresponding bins are added to /// the overflow bin. /// Statistics will be recomputed from the new bin contents. +/// +/// #### case 2 `xbins`!=0 || `ybins`!=0 +/// +/// A new histogram is created and `newname` must be specified. +/// For an axis with a non-null bin-edges array, `nxgroup` (`nygroup`) is the +/// number of bins of the new x-axis (y-axis) and `xbins` (`ybins`) must hold +/// the `nxgroup+1` (`nygroup+1`) edges of the new bins. An axis without an +/// array is rebinned in constant groups as in case 1. +/// The content of each old bin is added to the new bin containing its center; +/// old bins outside the range of the new axes end up in the under-/overflow +/// bins. Errors stored via Sumw2 are correctly recalculated. +/// +/// \note The new bin edges should line up with old bin edges: the entries of +/// an old bin that is split between two new bins are all transferred to the +/// bin containing the old bin center, and a warning is emitted. +/// +/// example: rebinning a TH2F with 100 x 100 bins into 24 x 24 variable bins +/// ~~~ {.cpp} +/// Double_t xbins[25] = {...}; // low-edges plus upper edge of last bin +/// Double_t ybins[25] = {...}; +/// TH2 *hnew = h2->Rebin2D(24, 24, "hnew", xbins, ybins); +/// ~~~ -TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname) +TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname, const Double_t *xbins, const Double_t *ybins) { - Int_t nxbins = fXaxis.GetNbins(); - Int_t nybins = fYaxis.GetNbins(); - Int_t nx = nxbins + 2; // normal bins + underflow and overflow - Int_t ny = nybins + 2; - Double_t xmin = fXaxis.GetXmin(); - Double_t xmax = fXaxis.GetXmax(); - Double_t ymin = fYaxis.GetXmin(); - Double_t ymax = fYaxis.GetXmax(); - if (GetDimension() != 2) { Error("Rebin2D", "Histogram must be TH2. This histogram has %d dimensions.", GetDimension()); return nullptr; } - if ((nxgroup <= 0) || (nxgroup > nxbins)) { - Error("Rebin2D", "Illegal value of nxgroup=%d",nxgroup); - return nullptr; + // something to do? + if (nxgroup == 1 && nygroup == 1 && !xbins && !ybins) { + return (newname && strlen(newname) > 0) ? (TH2 *)Clone(newname) : this; } - if ((nygroup <= 0) || (nygroup > nybins)) { - Error("Rebin2D", "Illegal value of nygroup=%d",nygroup); + if ((!newname || strlen(newname) == 0) && (xbins || ybins)) { + Error("Rebin2D", "if xbins or ybins are specified, newname must be given"); return nullptr; } - Int_t newxbins = nxbins / nxgroup; - Int_t newybins = nybins / nygroup; - Int_t newnx = newxbins + 2; // regular bins + overflow / underflow - Int_t newny = newybins + 2; // regular bins + overflow / underflow + const Int_t nxbins = fXaxis.GetNbins(); + const Int_t nybins = fYaxis.GetNbins(); - // Save old bin contents into a new array - Double_t *oldBins = new Double_t[fNcells]; - for (Int_t i = 0; i < fNcells; ++i) oldBins[i] = RetrieveBinContent(i); - - Double_t* oldErrors = nullptr; - if (fSumw2.fN) { - oldErrors = new Double_t[fNcells]; - for (Int_t i = 0; i < fNcells; ++i) oldErrors[i] = GetBinErrorSqUnchecked(i); + // validate the parameters and define the axes of the rebinned histogram + // and the mapping of old to new bins + ROOT::Internal::RebinnedAxisInfo infoX, infoY; + if (!ROOT::Internal::SetupRebinnedAxis(fXaxis, nxgroup, xbins, 'x', *this, "Rebin2D", infoX) || + !ROOT::Internal::SetupRebinnedAxis(fYaxis, nygroup, ybins, 'y', *this, "Rebin2D", infoY)) { + return nullptr; } + const Int_t newxbins = infoX.nNewBins; + const Int_t newybins = infoY.nNewBins; - // create a clone of the old histogram if newname is specified - TH2* hnew = this; - if (newname && strlen(newname)) { - hnew = (TH2*)Clone(); - hnew->SetName(newname); + // Save old bin contents into a new array + std::vector oldBins(fNcells); + for (Int_t i = 0; i < fNcells; ++i) + oldBins[i] = RetrieveBinContent(i); + + std::vector oldErrors; + if (fSumw2.fN != 0) { + oldErrors.resize(fNcells); + for (Int_t i = 0; i < fNcells; ++i) + oldErrors[i] = GetBinErrorSqUnchecked(i); } - bool resetStat = false; - - // change axis specs and rebuild bin contents array - if(newxbins * nxgroup != nxbins) { - xmax = fXaxis.GetBinUpEdge(newxbins * nxgroup); - resetStat = true; // stats must be reset because top bins will be moved to overflow bin - } - if(newybins * nygroup != nybins) { - ymax = fYaxis.GetBinUpEdge(newybins * nygroup); - resetStat = true; // stats must be reset because top bins will be moved to overflow bin + // rebinning will not redistribute under-/overflow content into the range + // of new axes that extend beyond the old ones + ROOT::Internal::WarnAboutUnusedFlowContent(fXaxis, infoX, xbins, 'X', oldBins.data(), 1, nybins + 2, nxbins + 2, + *this, "Rebin2D"); + ROOT::Internal::WarnAboutUnusedFlowContent(fYaxis, infoY, ybins, 'Y', oldBins.data(), nxbins + 2, nxbins + 2, 1, + *this, "Rebin2D"); + + // create a clone of the old histogram if newname is specified (guaranteed + // when bin edges are passed) + TH2 *hnew = this; + if (newname && strlen(newname) > 0) { + hnew = (TH2 *)Clone(newname); } // save the TAttAxis members (reset by SetBins) for x axis @@ -1750,68 +1772,25 @@ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname) Color_t yTitleColor = fYaxis.GetTitleColor(); Style_t yTitleFont = fYaxis.GetTitleFont(); + ROOT::Internal::SetRebinnedBins2D(*hnew, infoX.newAxis, infoY.newAxis); // changes also errors array (if any) - // copy merged bin contents (ignore under/overflows) - if (nxgroup != 1 || nygroup != 1) { - if(fXaxis.GetXbins()->GetSize() > 0 || fYaxis.GetXbins()->GetSize() > 0){ - // variable bin sizes in x or y, don't treat both cases separately - Double_t *xbins = new Double_t[newxbins + 1]; - for(Int_t i = 0; i <= newxbins; ++i) xbins[i] = fXaxis.GetBinLowEdge(1 + i * nxgroup); - Double_t *ybins = new Double_t[newybins + 1]; - for(Int_t i = 0; i <= newybins; ++i) ybins[i] = fYaxis.GetBinLowEdge(1 + i * nygroup); - hnew->SetBins(newxbins, xbins, newybins, ybins); // changes also errors array (if any) - delete [] xbins; - delete [] ybins; - } else { - hnew->SetBins(newxbins, xmin, xmax, newybins, ymin, ymax); //changes also errors array - } - - // (0, 0): x - underflow; y - underflow - hnew->UpdateBinContent(0, oldBins[0]); - if (oldErrors) hnew->fSumw2[0] = 0; - - // (x, 0): x - regular / overflow; y - underflow - for(Int_t binx = 1, oldbinx = 1; binx < newnx; ++binx, oldbinx += nxgroup){ - Double_t binContent = 0.0, binErrorSq = 0.0; - for (Int_t i = 0; i < nxgroup && (oldbinx + i) < nx; ++i) { - Int_t bin = oldbinx + i; - binContent += oldBins[bin]; - if(oldErrors) binErrorSq += oldErrors[bin]; - } - Int_t newbin = binx; - hnew->UpdateBinContent(newbin, binContent); - if (oldErrors) hnew->fSumw2[newbin] = binErrorSq; - } - - // (0, y): x - underflow; y - regular / overflow - for(Int_t biny = 1, oldbiny = 1; biny < newny; ++biny, oldbiny += nygroup){ - Double_t binContent = 0.0, binErrorSq = 0.0; - for (Int_t j = 0; j < nygroup && (oldbiny + j) < ny; ++j) { - Int_t bin = (oldbiny + j) * nx; - binContent += oldBins[bin]; - if(oldErrors) binErrorSq += oldErrors[bin]; - } - Int_t newbin = biny * newnx; - hnew->UpdateBinContent(newbin, binContent); - if (oldErrors) hnew->fSumw2[newbin] = binErrorSq; - } - - // (x, y): x - regular / overflow; y - regular / overflow - for (Int_t binx = 1, oldbinx = 1; binx < newnx; ++binx, oldbinx += nxgroup) { - for (Int_t biny = 1, oldbiny = 1; biny < newny; ++biny, oldbiny += nygroup) { - Double_t binContent = 0.0, binErrorSq = 0.0; - for (Int_t i = 0; i < nxgroup && (oldbinx + i) < nx; ++i) { - for (Int_t j = 0; j < nygroup && (oldbiny + j) < ny; ++j) { - Int_t bin = oldbinx + i + (oldbiny + j) * nx; - binContent += oldBins[bin]; - if (oldErrors) binErrorSq += oldErrors[bin]; - } - } - Int_t newbin = binx + biny * newnx; - hnew->UpdateBinContent(newbin, binContent); - if (oldErrors) hnew->fSumw2[newbin] = binErrorSq; - } - } + // add the content of each old cell (including under- and overflows) to + // the new cell that contains its bin center + const Int_t newncells = (newxbins + 2) * (newybins + 2); + std::vector newBins(newncells, 0.); + std::vector newErrors; + if (oldErrors.empty()) { + ROOT::Internal::MergeRebinnedCells(nxbins, nybins, newxbins, infoX.binMap, infoY.binMap, + {{oldBins.data(), newBins.data()}}); + } else { + newErrors.resize(newncells, 0.); + ROOT::Internal::MergeRebinnedCells(nxbins, nybins, newxbins, infoX.binMap, infoY.binMap, + {{oldBins.data(), newBins.data()}, {oldErrors.data(), newErrors.data()}}); + } + for (Int_t i = 0; i < newncells; ++i) { + hnew->UpdateBinContent(i, newBins[i]); + if (!oldErrors.empty()) + hnew->fSumw2[i] = newErrors[i]; } // Restore x axis attributes @@ -1839,14 +1818,14 @@ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname) fYaxis.SetTitleColor(yTitleColor); fYaxis.SetTitleFont(yTitleFont); - if (resetStat) hnew->ResetStats(); + // when the group count does not divide the old bin count, the top bins + // moved to the overflow: recompute the statistics from the bin contents + if (infoX.truncated || infoY.truncated) + hnew->ResetStats(); - delete [] oldBins; - if (oldErrors) delete [] oldErrors; return hnew; } - //////////////////////////////////////////////////////////////////////////////// TProfile *TH2::DoProfile(bool onX, const char *name, Int_t firstbin, Int_t lastbin, Option_t *option) const diff --git a/hist/hist/src/TProfile2D.cxx b/hist/hist/src/TProfile2D.cxx index 09e6700e8989e..fa2e5ef736c6d 100644 --- a/hist/hist/src/TProfile2D.cxx +++ b/hist/hist/src/TProfile2D.cxx @@ -16,7 +16,11 @@ #include "TError.h" #include "TClass.h" #include "TProfileHelper.h" +#include "Rebin2DHelpers.h" + +#include #include +#include Bool_t TProfile2D::fgApproximate = kFALSE; @@ -1551,9 +1555,11 @@ void TProfile2D::ExtendAxis(Double_t x, TAxis *axis) //////////////////////////////////////////////////////////////////////////////// /// Rebin this histogram grouping nxgroup/nygroup bins along the xaxis/yaxis together. /// -/// if newname is not blank a new profile hnew is created. +/// ## case 1 `xbins`=0 || `ybins`=0 +/// +/// if `newname` is not blank a new profile hnew is created. /// else the current histogram is modified (default) -/// The parameter nxgroup/nygroup indicate how many bins along the xaxis/yaxis of this +/// The parameters `nxgroup`/`nygroup` indicate how many bins along the xaxis/yaxis of this /// have to be merged into one bin of hnew /// If the original profile has errors stored (via Sumw2), the resulting /// profile has new errors correctly calculated. @@ -1569,292 +1575,116 @@ void TProfile2D::ExtendAxis(Double_t x, TAxis *axis) /// // merging 5 bins of hpxpy along the yaxis in one bin /// ~~~ /// -/// NOTE : If nxgroup/nygroup is not an exact divider of the number of bins, +/// \note If `nxgroup`/`nygroup` is not an exact divider of the number of bins, /// along the xaxis/yaxis the top limit(s) of the rebinned profile /// is changed to the upper edge of the xbin=newxbins*nxgroup resp. /// ybin=newybins*nygroup and the remaining bins are added to /// the overflow bin. /// Statistics will be recomputed from the new bin contents. +/// +/// ## case 2 `xbins`!=0 || `ybins`!=0 +/// +/// A new profile is created and `newname` must be specified. +/// For an axis with a non-null bin-edges array, `nxgroup` (`nygroup`) is the +/// number of bins of the new x-axis (y-axis) and `xbins` (`ybins`) must hold +/// the `nxgroup+1` (`nygroup+1`) edges of the new bins. An axis without an +/// array is rebinned in constant groups as in case 1. +/// The data of each old bin are added to the new bin containing its center; +/// old bins outside the range of the new axes end up in the under-/overflow +/// bins. +/// +/// \note The new bin edges should line up with old bin edges: the entries of +/// an old bin that is split between two new bins are all transferred to the +/// bin containing the old bin center, and a warning is emitted. +/// +/// example: rebinning a TProfile2D with 100 x 100 bins into 24 x 24 variable bins +/// ~~~ {.cpp} +/// Double_t xbins[25] = {...}; // low-edges plus upper edge of last bin +/// Double_t ybins[25] = {...}; +/// TProfile2D *hpnew = hp->Rebin2D(24, 24, "hpnew", xbins, ybins); +/// ~~~ -TProfile2D * TProfile2D::Rebin2D(Int_t nxgroup ,Int_t nygroup,const char * newname ) { +TProfile2D * +TProfile2D::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname, const Double_t *xbins, const Double_t *ybins) +{ //something to do? - if((nxgroup != 1) || (nygroup != 1)){ - Int_t nxbins = fXaxis.GetNbins(); - Int_t nybins = fYaxis.GetNbins(); - Double_t xmin = fXaxis.GetXmin(); - Double_t xmax = fXaxis.GetXmax(); - Double_t ymin = fYaxis.GetXmin(); - Double_t ymax = fYaxis.GetXmax(); - if ((nxgroup <= 0) || (nxgroup > nxbins)) { - Error("Rebin", "Illegal value of nxgroup=%d",nxgroup); - return nullptr; - } - if ((nygroup <= 0) || (nygroup > nybins)) { - Error("Rebin", "Illegal value of nygroup=%d",nygroup); - return nullptr; - } - - Int_t newxbins = nxbins/nxgroup; - Int_t newybins = nybins/nygroup; - - //warning if bins are added to the overflow bin - if(newxbins*nxgroup != nxbins) { - Warning("Rebin", "nxgroup=%d should be an exact divider of nxbins=%d",nxgroup,nxbins); - } - if(newybins*nygroup != nybins) { - Warning("Rebin", "nygroup=%d should be an exact divider of nybins=%d",nygroup,nybins); - } - - //save old bin contents in new arrays - Double_t *oldBins = new Double_t[(nxbins+2)*(nybins+2)]; - Double_t *oldCount = new Double_t[(nxbins+2)*(nybins+2)]; - Double_t *oldErrors = new Double_t[(nxbins+2)*(nybins+2)]; - Double_t *oldBinw2 = (fBinSumw2.fN ? new Double_t[(nxbins+2)*(nybins+2)] : nullptr ); - Double_t *cu1 = GetW(); - Double_t *er1 = GetW2(); - Double_t *en1 = GetB(); - Double_t *ew1 = GetB2(); - for(Int_t ibin=0; ibin < (nxbins+2)*(nybins+2); ibin++){ - oldBins[ibin] = cu1[ibin]; - oldCount[ibin] = en1[ibin]; - oldErrors[ibin] = er1[ibin]; - if (ew1 && fBinSumw2.fN) oldBinw2[ibin] = ew1[ibin]; - } + if ((nxgroup == 1) && (nygroup == 1) && !xbins && !ybins) { + if (newname && (strlen(newname) > 0)) + return (TProfile2D *)Clone(newname); + else + return this; + } - // create a clone of the old profile if newname is specified - TProfile2D *hnew = this; - if(newname && strlen(newname) > 0) { - hnew = (TProfile2D*)Clone(newname); - } + if ((!newname || strlen(newname) == 0) && (xbins || ybins)) { + Error("Rebin2D", "if xbins or ybins are specified, newname must be given"); + return nullptr; + } - // in case of nxgroup/nygroup not an exact divider of nxbins/nybins, - // top limit is changed (see NOTE in method comment) - if(newxbins*nxgroup != nxbins) { - xmax = fXaxis.GetBinUpEdge(newxbins*nxgroup); - hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin - } - if(newybins*nygroup != nybins) { - ymax = fYaxis.GetBinUpEdge(newybins*nygroup); - hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin - } + const Int_t nxbins = fXaxis.GetNbins(); + const Int_t nybins = fYaxis.GetNbins(); - //rebin the axis - if((fXaxis.GetXbins()->GetSize() > 0) || (fYaxis.GetXbins()->GetSize() > 0)){ - Double_t* xbins = new Double_t[newxbins+1]; - Double_t* ybins = new Double_t[newybins+1]; - for(Int_t i=0; i < newxbins+1; i++) - xbins[i] = fXaxis.GetBinLowEdge(1+i*nxgroup); - for(Int_t j=0; j < newybins+1; j++) - ybins[j] = fYaxis.GetBinLowEdge(1+j*nygroup); - hnew->SetBins(newxbins,xbins,newybins,ybins); - delete [] xbins; - delete [] ybins; - } - //fixed bin size - else{ - hnew->SetBins(newxbins,xmin,xmax,newybins,ymin,ymax); - } + // validate the parameters and define the axes of the rebinned profile and + // the mapping of old to new bins + ROOT::Internal::RebinnedAxisInfo infoX, infoY; + if (!ROOT::Internal::SetupRebinnedAxis(fXaxis, nxgroup, xbins, 'x', *this, "Rebin2D", infoX) || + !ROOT::Internal::SetupRebinnedAxis(fYaxis, nygroup, ybins, 'y', *this, "Rebin2D", infoY)) { + return nullptr; + } + const Int_t newxbins = infoX.nNewBins; + const Int_t newybins = infoY.nNewBins; - //merge bins - Double_t *cu2 = hnew->GetW(); - Double_t *er2 = hnew->GetW2(); - Double_t *en2 = hnew->GetB(); + // save old bin contents in new arrays + const Int_t ncells = (nxbins + 2) * (nybins + 2); + std::vector oldBins(GetW(), GetW() + ncells); + std::vector oldErrors(GetW2(), GetW2() + ncells); + std::vector oldCount(GetB(), GetB() + ncells); + std::vector oldBinw2; + if (fBinSumw2.fN) + oldBinw2.assign(GetB2(), GetB2() + ncells); + + // rebinning will not redistribute under-/overflow content into the range + // of new axes that extend beyond the old ones + ROOT::Internal::WarnAboutUnusedFlowContent(fXaxis, infoX, xbins, 'X', oldBins.data(), 1, nybins + 2, nxbins + 2, + *this, "Rebin2D"); + ROOT::Internal::WarnAboutUnusedFlowContent(fYaxis, infoY, ybins, 'Y', oldBins.data(), nxbins + 2, nxbins + 2, 1, + *this, "Rebin2D"); + + // create a clone of the old profile if newname is specified (guaranteed + // when bin edges are passed) + TProfile2D *hnew = this; + if (newname && strlen(newname) > 0) { + hnew = (TProfile2D *)Clone(newname); + } + + // when the group count does not divide the old bin count, the top bins + // move to the overflow and the stats must be recomputed + if (infoX.truncated || infoY.truncated) + hnew->fTsumw = 0; + + // rebin the axes + ROOT::Internal::SetRebinnedBins2D(*hnew, infoX.newAxis, infoY.newAxis); + + // merge bins: add the content of each old cell (including under- and + // overflow) to the new cell that contains its bin center + const Int_t newncells = (newxbins + 2) * (newybins + 2); + Double_t *cu2 = hnew->GetW(); + Double_t *er2 = hnew->GetW2(); + Double_t *en2 = hnew->GetB(); + std::fill(cu2, cu2 + newncells, 0.); + std::fill(er2, er2 + newncells, 0.); + std::fill(en2, en2 + newncells, 0.); + if (fBinSumw2.fN) { Double_t *ew2 = hnew->GetB2(); - Double_t binContent, binCount, binError, binSumw2; - //connection between x and y bin number and linear global bin number: - //global bin = xbin + (nxbins+2) * ybin - Int_t oldxbin = 1; - Int_t oldybin = 1; - //global bin number - Int_t bin; - for(Int_t xbin = 1; xbin <= newxbins; xbin++){ - oldybin = 1; - for(Int_t ybin = 1; ybin <= newybins; ybin++){ - binContent = 0; - binCount = 0; - binError = 0; - binSumw2 = 0; - for(Int_t i=0; i < nxgroup; i++){ - if(oldxbin + i > nxbins) break; - for(Int_t j=0; j < nygroup; j++){ - if(oldybin + j > nybins) break; - bin = oldxbin + i + (nxbins+2)*(oldybin+j); - binContent += oldBins[bin]; - binCount += oldCount[bin]; - binError += oldErrors[bin]; - if(fBinSumw2.fN) binSumw2 += oldBinw2[bin]; - } - } - bin = xbin + (newxbins + 2)*ybin; - cu2[bin] = binContent; - er2[bin] = binError; - en2[bin] = binCount; - if(fBinSumw2.fN) ew2[bin] = binSumw2; - oldybin += nygroup; - } - oldxbin += nxgroup; - } - - //copy the underflow bin in x and y (0,0) - cu2[0] = oldBins[0]; - er2[0] = oldErrors[0]; - en2[0] = oldCount[0]; - if(fBinSumw2.fN) ew2[0] = oldBinw2[0]; - //calculate overflow bin in x and y (newxbins+1,newybins+1) - //therefore the oldxbin and oldybin from above are needed! - binContent = 0; - binCount = 0; - binError = 0; - binSumw2 = 0; - for(Int_t i=oldxbin; i <= nxbins+1; i++){ - for(Int_t j=oldybin; j <= nybins+1; j++){ - //global bin number - bin = i + (nxbins+2)*j; - binContent += oldBins[bin]; - binCount += oldCount[bin]; - binError += oldErrors[bin]; - if(fBinSumw2.fN) binSumw2 += oldBinw2[bin]; - } - } - bin = (newxbins+2)*(newybins+2)-1; - cu2[bin] = binContent; - er2[bin] = binError; - en2[bin] = binCount; - if(fBinSumw2.fN) ew2[bin] = binSumw2; - //calculate overflow bin in x and underflow bin in y (newxbins+1,0) - binContent = 0; - binCount = 0; - binError = 0; - binSumw2 = 0; - for(Int_t i=oldxbin; i <= nxbins+1; i++){ - bin = i; - binContent += oldBins[bin]; - binCount += oldCount[bin]; - binError += oldErrors[bin]; - if(fBinSumw2.fN) binSumw2 += oldBinw2[bin]; - } - bin = newxbins + 1; - cu2[bin] = binContent; - er2[bin] = binError; - en2[bin] = binCount; - if(fBinSumw2.fN) ew2[bin] = binSumw2; - //calculate underflow bin in x and overflow bin in y (0,newybins+1) - binContent = 0; - binCount = 0; - binError = 0; - binSumw2 = 0; - for(Int_t i=oldybin; i <= nybins+1; i++){ - bin = i*(nxbins + 2); - binContent += oldBins[bin]; - binCount += oldCount[bin]; - binError += oldErrors[bin]; - if(fBinSumw2.fN) binSumw2 += oldBinw2[bin]; - } - bin = (newxbins + 2)*(newybins + 1); - cu2[bin] = binContent; - er2[bin] = binError; - en2[bin] = binCount; - if(fBinSumw2.fN) ew2[bin] = binSumw2; - //calculate under/overflow contents in y for the new x bins - Double_t binContentuf, binCountuf, binErroruf, binSumw2uf; - Double_t binContentof, binCountof, binErrorof, binSumw2of; - Int_t ufbin, ofbin; - Int_t oldxbin2 = 1; - for(Int_t xbin = 1; xbin <= newxbins; xbin++){ - binContentuf = 0; - binCountuf = 0; - binErroruf = 0; - binSumw2uf = 0; - binContentof = 0; - binCountof = 0; - binErrorof = 0; - binSumw2of = 0; - for(Int_t i = 0; i < nxgroup; i++){ - //index of under/overflow bin for y in old binning - ufbin = (oldxbin2 + i); - binContentuf += oldBins[ufbin]; - binCountuf += oldCount[ufbin]; - binErroruf += oldErrors[ufbin]; - if(fBinSumw2.fN) binSumw2uf += oldBinw2[ufbin]; - for(Int_t j = oldybin; j <= nybins+1; j++) - { - ofbin = ufbin + j*(nxbins + 2); - binContentof += oldBins[ofbin]; - binCountof += oldCount[ofbin]; - binErrorof += oldErrors[ofbin]; - if(fBinSumw2.fN) binSumw2of += oldBinw2[ofbin]; - } - } - //index of under/overflow bin for y in new binning - ufbin = xbin; - ofbin = ufbin + (newybins + 1)*(newxbins + 2); - cu2[ufbin] = binContentuf; - er2[ufbin] = binErroruf; - en2[ufbin] = binCountuf; - if(fBinSumw2.fN) ew2[ufbin] = binSumw2uf; - cu2[ofbin] = binContentof; - er2[ofbin] = binErrorof; - en2[ofbin] = binCountof; - if(fBinSumw2.fN) ew2[ofbin] = binSumw2of; - - oldxbin2 += nxgroup; - } - //calculate under/overflow contents in x for the new y bins - Int_t oldybin2 = 1; - for(Int_t ybin = 1; ybin <= newybins; ybin++){ - binContentuf = 0; - binCountuf = 0; - binErroruf = 0; - binSumw2uf = 0; - binContentof = 0; - binCountof = 0; - binErrorof = 0; - binSumw2of = 0; - for(Int_t i = 0; i < nygroup; i++){ - //index of under/overflow bin for x in old binning - ufbin = (oldybin2 + i)*(nxbins+2); - binContentuf += oldBins[ufbin]; - binCountuf += oldCount[ufbin]; - binErroruf += oldErrors[ufbin]; - if(fBinSumw2.fN) binSumw2uf += oldBinw2[ufbin]; - for(Int_t j = oldxbin; j <= nxbins+1; j++) - { - ofbin = j + ufbin; - binContentof += oldBins[ofbin]; - binCountof += oldCount[ofbin]; - binErrorof += oldErrors[ofbin]; - if(fBinSumw2.fN) binSumw2of += oldBinw2[ofbin]; - } - } - //index of under/overflow bin for x in new binning - ufbin = ybin * (newxbins + 2); - ofbin = newxbins + 1 + ufbin; - cu2[ufbin] = binContentuf; - er2[ufbin] = binErroruf; - en2[ufbin] = binCountuf; - if(fBinSumw2.fN) ew2[ufbin] = binSumw2uf; - cu2[ofbin] = binContentof; - er2[ofbin] = binErrorof; - en2[ofbin] = binCountof; - if(fBinSumw2.fN) ew2[ofbin] = binSumw2of; - - oldybin2 += nygroup; - } - - delete [] oldBins; - delete [] oldCount; - delete [] oldErrors; - if (oldBinw2) delete [] oldBinw2; - - return hnew; - } - //nxgroup == nygroup == 1 - else{ - if(newname && (strlen(newname) > 0)) - return (TProfile2D*)Clone(newname); - else - return this; + std::fill(ew2, ew2 + newncells, 0.); + ROOT::Internal::MergeRebinnedCells( + nxbins, nybins, newxbins, infoX.binMap, infoY.binMap, + {{oldBins.data(), cu2}, {oldErrors.data(), er2}, {oldCount.data(), en2}, {oldBinw2.data(), ew2}}); + } else { + ROOT::Internal::MergeRebinnedCells(nxbins, nybins, newxbins, infoX.binMap, infoY.binMap, + {{oldBins.data(), cu2}, {oldErrors.data(), er2}, {oldCount.data(), en2}}); } + + return hnew; } //////////////////////////////////////////////////////////////////////////////// diff --git a/hist/hist/test/CMakeLists.txt b/hist/hist/test/CMakeLists.txt index 3a9a1623c6d42..f970aa730ff7f 100644 --- a/hist/hist/test/CMakeLists.txt +++ b/hist/hist/test/CMakeLists.txt @@ -13,6 +13,7 @@ ROOT_ADD_GTEST(testTH2PolyGetNumberOfBins test_TH2Poly_GetNumberOfBins.cxx LIBRA ROOT_ADD_GTEST(testTHn THn.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTH1 test_TH1.cxx LIBRARIES Hist) ROOT_ADD_GTEST(testTH2 test_TH2.cxx LIBRARIES Hist) +ROOT_ADD_GTEST(testRebin2D test_rebin2D.cxx LIBRARIES Hist MathCore) ROOT_ADD_GTEST(testTH3 test_TH3.cxx LIBRARIES Hist) ROOT_ADD_GTEST(testTHStack test_THStack.cxx LIBRARIES Hist) ROOT_ADD_GTEST(testProject3Dname test_Project3D_name.cxx LIBRARIES Hist) diff --git a/hist/hist/test/test_rebin2D.cxx b/hist/hist/test/test_rebin2D.cxx new file mode 100644 index 0000000000000..50da82c025813 --- /dev/null +++ b/hist/hist/test/test_rebin2D.cxx @@ -0,0 +1,218 @@ +// Tests for rebinning TH2 and TProfile2D with variable bin edges (ROOT-5224). +// The rebinned histograms are compared against reference histograms created +// directly with the target binning and filled with the same pseudo-data. + +#include "gtest/gtest.h" + +#include "ROOT/TestSupport.hxx" + +#include "TH2D.h" +#include "TProfile2D.h" +#include "TRandom3.h" + +#include +#include +#include + +namespace { + +void fillOne(TH2D &h, double x, double y, double, double w) +{ + h.Fill(x, y, w); +} + +void fillOne(TProfile2D &p, double x, double y, double z, double w) +{ + p.Fill(x, y, z, w); +} + +// Fill both histograms with the same weighted pseudo-data, including entries +// that end up in the under- and overflow bins. +template +void fillSame(Hist &a, Hist &b) +{ + TRandom3 rng(42); + for (int i = 0; i < 10000; ++i) { + const double x = rng.Uniform(-10., 110.); + const double y = rng.Uniform(-10., 110.); + const double z = rng.Gaus(5., 1.); + const double w = rng.Uniform(0.1, 2.); + fillOne(a, x, y, z, w); + fillOne(b, x, y, z, w); + } +} + +void expectSameBinsAndStats(const TH2 &h, const TH2 &ref, const std::string &ctx) +{ + ASSERT_EQ(h.GetNbinsX(), ref.GetNbinsX()) << ctx; + ASSERT_EQ(h.GetNbinsY(), ref.GetNbinsY()) << ctx; + for (Int_t j = 0; j <= ref.GetNbinsY() + 1; ++j) { + for (Int_t i = 0; i <= ref.GetNbinsX() + 1; ++i) { + const std::string binCtx = ctx + " bin (" + std::to_string(i) + "," + std::to_string(j) + ")"; + const Int_t bin = ref.GetBin(i, j); + const double refContent = ref.GetBinContent(bin); + const double refError = ref.GetBinError(bin); + EXPECT_NEAR(h.GetBinContent(bin), refContent, 1e-6 * std::max(1., std::abs(refContent))) << binCtx; + EXPECT_NEAR(h.GetBinError(bin), refError, 1e-6 * std::max(1., refError)) << binCtx; + } + } + EXPECT_NEAR(h.GetEntries(), ref.GetEntries(), 1e-6) << ctx; + EXPECT_NEAR(h.GetMean(1), ref.GetMean(1), 1e-9) << ctx; + EXPECT_NEAR(h.GetMean(2), ref.GetMean(2), 1e-9) << ctx; + EXPECT_NEAR(h.GetStdDev(1), ref.GetStdDev(1), 1e-9) << ctx; +} + +const double xEdges[5] = {0., 10., 40., 70., 100.}; +const double yEdges[5] = {0., 30., 50., 90., 100.}; + +} // namespace + +// Rebin a uniform TH2 into variable bins on both axes. +TEST(Rebin2DVariable, TH2BothAxes) +{ + TH2D fine("fine", "fine", 100, 0., 100., 100, 0., 100.); + fine.Sumw2(); + TH2D ref("ref", "ref", 4, xEdges, 4, yEdges); + ref.Sumw2(); + fillSame(fine, ref); + + std::unique_ptr hnew{fine.Rebin2D(4, 4, "hnew", xEdges, yEdges)}; + ASSERT_NE(hnew, nullptr); + ASSERT_NE(hnew.get(), &fine); + expectSameBinsAndStats(*hnew, ref, "TH2BothAxes"); +} + +// Variable bins on the x-axis only: the y-axis is grouped by a constant factor. +TEST(Rebin2DVariable, TH2SingleAxis) +{ + TH2D fine("fine", "fine", 100, 0., 100., 100, 0., 100.); + fine.Sumw2(); + TH2D ref("ref", "ref", 4, xEdges, 50, 0., 100.); + ref.Sumw2(); + fillSame(fine, ref); + + std::unique_ptr hnew{fine.Rebin2D(4, 2, "hnew", xEdges, nullptr)}; + ASSERT_NE(hnew, nullptr); + expectSameBinsAndStats(*hnew, ref, "TH2SingleAxis"); + + // TH2::Rebin with bin edges must forward to Rebin2D with the TH1 + // conventions, leaving the y-axis untouched + std::unique_ptr viaRebin{fine.Rebin(4, "viaRebin", xEdges)}; + std::unique_ptr viaRebin2D{fine.Rebin2D(4, 1, "viaRebin2D", xEdges, nullptr)}; + ASSERT_NE(viaRebin, nullptr); + ASSERT_NE(viaRebin2D, nullptr); + expectSameBinsAndStats(*viaRebin, *viaRebin2D, "TH2RebinForwards"); +} + +// Constant-group rebinning of a TH2 that has variable-width axes: the new +// edges are synthesized from the old axis. +TEST(Rebin2DVariable, TH2VariableSourceConstantGroups) +{ + TH2D fine("fine", "fine", 4, xEdges, 4, yEdges); + fine.Sumw2(); + const double xCoarse[3] = {0., 40., 100.}; + const double yCoarse[3] = {0., 50., 100.}; + TH2D ref("ref", "ref", 2, xCoarse, 2, yCoarse); + ref.Sumw2(); + fillSame(fine, ref); + + std::unique_ptr hnew{fine.Rebin2D(2, 2, "hnew")}; + ASSERT_NE(hnew, nullptr); + expectSameBinsAndStats(*hnew, ref, "TH2VariableSourceConstantGroups"); +} + +// Rebin a uniform TProfile2D into variable bins on both axes. +TEST(Rebin2DVariable, Profile2DBothAxes) +{ + TProfile2D fine("fine", "fine", 100, 0., 100., 100, 0., 100.); + fine.Sumw2(); + TProfile2D ref("ref", "ref", 4, xEdges, 4, yEdges); + ref.Sumw2(); + fillSame(fine, ref); + + std::unique_ptr pnew{fine.Rebin2D(4, 4, "pnew", xEdges, yEdges)}; + ASSERT_NE(pnew, nullptr); + expectSameBinsAndStats(*pnew, ref, "Profile2DBothAxes"); + for (Int_t j = 0; j <= ref.GetNbinsY() + 1; ++j) { + for (Int_t i = 0; i <= ref.GetNbinsX() + 1; ++i) { + const Int_t bin = ref.GetBin(i, j); + EXPECT_NEAR(pnew->GetBinEntries(bin), ref.GetBinEntries(bin), 1e-6) + << "Profile2DBothAxes entries bin (" << i << "," << j << ")"; + } + } +} + +// Constant-group rebinning of a TProfile2D that has variable-width axes +// (regression test: this used to pass null bin edges to SetBins). +TEST(Rebin2DVariable, Profile2DVariableSourceConstantGroups) +{ + TProfile2D fine("fine", "fine", 4, xEdges, 4, yEdges); + const double xCoarse[3] = {0., 40., 100.}; + const double yCoarse[3] = {0., 50., 100.}; + TProfile2D ref("ref", "ref", 2, xCoarse, 2, yCoarse); + fillSame(fine, ref); + + std::unique_ptr pnew{fine.Rebin2D(2, 2, "pnew")}; + ASSERT_NE(pnew, nullptr); + expectSameBinsAndStats(*pnew, ref, "Profile2DVariableSourceConstantGroups"); +} + +// When the group count does not divide the number of bins, the top bins move +// to the overflow and the statistics must be recomputed from the bin contents. +TEST(Rebin2DVariable, TH2NonDivisorGroupStats) +{ + TH2D h("h", "h", 10, 0., 10., 10, 0., 10.); + TRandom3 rng(43); + for (int i = 0; i < 10000; ++i) + h.Fill(rng.Uniform(0., 10.), rng.Uniform(0., 10.)); + + { + ROOT::TestSupport::CheckDiagsRAII checkDiag; + checkDiag.requiredDiag(kWarning, "TH2D::Rebin2D", "is not an exact divider", false); + ASSERT_EQ(h.Rebin2D(3, 3), &h); // in-place + } + ASSERT_EQ(h.GetNbinsX(), 3); + ASSERT_DOUBLE_EQ(h.GetXaxis()->GetXmax(), 9.); + + // reference statistics recomputed from the rebinned contents + double sw = 0., swx = 0., swx2 = 0.; + for (Int_t j = 1; j <= h.GetNbinsY(); ++j) { + for (Int_t i = 1; i <= h.GetNbinsX(); ++i) { + const double c = h.GetBinContent(i, j); + const double x = h.GetXaxis()->GetBinCenter(i); + sw += c; + swx += c * x; + swx2 += c * x * x; + } + } + EXPECT_NEAR(h.GetMean(1), swx / sw, 1e-9); + EXPECT_NEAR(h.GetStdDev(1), std::sqrt(swx2 / sw - (swx / sw) * (swx / sw)), 1e-9); +} + +// Variable-bin rebinning requires a new name, and warns when a new bin edge +// does not line up with an old bin edge. +TEST(Rebin2DVariable, Diagnostics) +{ + TH2D h("h", "h", 100, 0., 100., 100, 0., 100.); + { + // an empty name must be rejected like a null one, otherwise Clone("") + // would create a second histogram registered under the original name + ROOT::TestSupport::CheckDiagsRAII checkDiag(kError, "TH2D::Rebin2D", "newname must be given", false); + EXPECT_EQ(h.Rebin2D(4, 4, nullptr, xEdges, yEdges), nullptr); + EXPECT_EQ(h.Rebin2D(4, 4, "", xEdges, yEdges), nullptr); + } + { + ROOT::TestSupport::CheckDiagsRAII checkDiag(kWarning, "TH2D::Rebin2D", "does not match any bin edges", false); + const double misaligned[3] = {0., 10.5, 100.}; + std::unique_ptr hnew{h.Rebin2D(2, 2, "hnew", misaligned, nullptr)}; + EXPECT_NE(hnew, nullptr); + } + { + // a new top edge that splits an old bin between range and overflow must + // warn as well + ROOT::TestSupport::CheckDiagsRAII checkDiag(kWarning, "TH2D::Rebin2D", "does not match any bin edges", false); + const double splitTop[3] = {0., 50., 99.5}; + std::unique_ptr hnew{h.Rebin2D(2, 2, "hnew", splitTop, nullptr)}; + EXPECT_NE(hnew, nullptr); + } +}