Implement variable rebinning 2-D histogram classes (ROOT-5224)

hist/hist/inc/TH2.h

@@ -107,7 +107,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*

hist/hist/inc/TProfile2D.h

@@ -131,7 +131,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;

hist/hist/src/TH2.cxx

@@ -1613,7 +1613,7 @@ TH2 *TH2::RebinY(Int_t ngroup, const char *newname)
 /// If a non-null pointer is given an error is flagged
 /// 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");
@@ -1622,12 +1622,15 @@ TH2 * TH2::Rebin( Int_t ngroup, const char*newname, const Double_t *xbins)
    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.
@@ -1641,14 +1644,41 @@ 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 (you should specify `newname`).
+/// The parameter `nxgroup` (`nygroup`) is the number of variable size bins for the x-axis
+/// (y-axis) in the created histogram.
+/// The arrays `xbins` and `ybins` must contain `nxgroup+1` and `nygroups+1` elements
+/// that represent the low-edges of the x and y bins respectively.
+/// If the original histogram has errors stored (via Sumw2), the resulting
+/// histograms has new errors correctly calculated.
+///
+/// \note  The bin edges specified in xbins and ybins should correspond
+/// to bin edges in the original histogram. If a bin edge in the new histogram
+/// is in the middle of a bin in the original histogram, all entries in
+/// the split bin in the original histogram will be transfered to the
+/// lower of the two possible bins in the new histogram. This is
+/// probably not what you want. A warning message is emitted in this
+/// case.
+///
+/// examples: if h2 is an existing TH2F histogram with 100 bins on x-axis
+/// and 100 bins y-axis
+///
+/// ~~~ {.cpp}
+///     Double_t xbins[25] = {...} array of low-edges for x-axis (xbins[25] is the upper edge of last bin)
+///     Double_t ybins[25] = {...} array of low-edges for y-axis (ybins[25] is the upper edge of last bin)
+///     h1->Rebin(24,24,"hnew",xbins,ybins);  //creates a new variable bin size histogram hnew
+/// ~~~
 
-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();
@@ -1671,37 +1701,94 @@ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname)
       Error("Rebin2D", "Illegal value of nygroup=%d",nygroup);
       return nullptr;
    }
+   if (!newname && xbins) {
+       Error("Rebin2D","if xbins is specified, newname must be given");
+       return 0;
+    }
+   if (!newname && ybins) {
+       Error("Rebin2D","if ybins is specified, newname must be given");
+       return 0;
+    }
 
    Int_t newxbins = nxbins / nxgroup;
-   Int_t newybins = nybins / nygroup;
    Int_t newnx = newxbins + 2; // regular bins + overflow / underflow
+   if (!xbins) {
+      Int_t nbgx = nxbins/nxgroup;
+      if (nbgx*nxgroup != nxbins) {
+         Warning("Rebin2D", "nxgroup=%d is not an exact divider of nxbins=%d.",nxgroup,nxbins);
+      }
+   }
+   else {
+   // in the case that xbins is given (rebinning in variable bins), nxgroup is
+   // the new number of bins and number of grouped bins is not constant.
+   // when looping for setting the contents for the new histogram we
+   // need to loop on all bins of original histogram. Then set nxgroup=nxbins
+      newxbins = nxgroup;
+      nxgroup = nxbins;
+   }
+
+   Int_t newybins = nybins / nygroup;
    Int_t newny = newybins + 2; // regular bins + overflow / underflow
+   if (!ybins) {
+      Int_t nbgy = nybins/nygroup;
+      if (nbgy*nygroup != nybins) {
+         Warning("Rebin2D", "nygroup=%d is not an exact divider of nybins=%d.",nygroup,nybins);
+      }
+   }
+   else {
+   // in the case that ybins is given (rebinning in variable bins), nygroup is
+   // the new number of bins and number of grouped bins is not constant.
+   // when looping for setting the contents for the new histogram we
+   // need to loop on all bins of original histogram. Then set nygroup=nybins
+      newybins = nygroup;
+      nygroup = nybins;
+   }
 
    // Save old bin contents into a new array
+   Double_t entries = fEntries;
    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) {
+   if (fSumw2.fN != 0) {
       oldErrors = new Double_t[fNcells];
       for (Int_t i = 0; i < fNcells; ++i) oldErrors[i] = GetBinErrorSqUnchecked(i);
    }
 
+   // rebin will not include underflow/overflow if new axis range is larger than old axis range
+   if (xbins) {
+      if (xbins[0] < fXaxis.GetXmin() && oldBins[0] != 0 )
+         Warning("Rebin2D","underflow entries for X axis will not be used when rebinning");
+      if (xbins[newxbins] > fXaxis.GetXmax() && oldBins[nxbins+1] != 0 )
+         Warning("Rebin2D","overflow entries for X axis will not be used when rebinning");
+   }
+   if (ybins) {
+       if (ybins[0] < fYaxis.GetXmin() && oldBins[0] != 0 )
+          Warning("Rebin2D","underflow entries for Y axis will not be used when rebinning");
+       if (ybins[newybins] > fYaxis.GetXmax() && oldBins[nybins+1] != 0 )
+          Warning("Rebin2D","overflow entries for Y axis will not be used when rebinning");
+    }
+
    // create a clone of the old histogram if newname is specified
    TH2* hnew = this;
-   if (newname && strlen(newname)) {
-      hnew = (TH2*)Clone();
-      hnew->SetName(newname);
+   if ((newname && strlen(newname)) || xbins || ybins) {
+      hnew = (TH2*)Clone(newname);
    }
 
+   //reset can extend bit to avoid an axis extension in SetBinContent
+   UInt_t oldExtendBitMask = hnew->SetCanExtend(kNoAxis);
+
+   // save original statistics
+   Double_t stat[kNstat];
+   GetStats(stat);
    bool resetStat = false;
 
    // change axis specs and rebuild bin contents array
-   if(newxbins * nxgroup != nxbins) {
+   if(!xbins && (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) {
+   if(!ybins && (newybins * nygroup != nybins)) {
       ymax = fYaxis.GetBinUpEdge(newybins * nygroup);
       resetStat = true; // stats must be reset because top bins will be moved to overflow bin
    }
@@ -1734,15 +1821,17 @@ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname)
 
    // copy merged bin contents (ignore under/overflows)
    if (nxgroup != 1 || nygroup != 1) {
-      if(fXaxis.GetXbins()->GetSize() > 0 || fYaxis.GetXbins()->GetSize() > 0){
+      if((!xbins && fXaxis.GetXbins()->GetSize() > 0) || (!ybins && 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;
+         Double_t *xbinsTmp = new Double_t[newxbins + 1];
+         for(Int_t i = 0; i <= newxbins; ++i) xbinsTmp[i] = fXaxis.GetBinLowEdge(1 + i * nxgroup);
+         Double_t *ybinsTmp = new Double_t[newybins + 1];
+         for(Int_t i = 0; i <= newybins; ++i) ybinsTmp[i] = fYaxis.GetBinLowEdge(1 + i * nygroup);
+         hnew->SetBins(newxbins, xbinsTmp, newybins, ybinsTmp); // changes also errors array (if any)
+         delete [] xbinsTmp;
+         delete [] ybinsTmp;
+      } else if(xbins && ybins) {
+          hnew->SetBins(newxbins,xbins, newybins, ybins);
       } else {
          hnew->SetBins(newxbins, xmin, xmax, newybins, ymin, ymax); //changes also errors array
       }
@@ -1820,7 +1909,11 @@ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *newname)
    fYaxis.SetTitleColor(yTitleColor);
    fYaxis.SetTitleFont(yTitleFont);
 
-   if (resetStat) hnew->ResetStats();
+   hnew->SetCanExtend(oldExtendBitMask); // restore previous state
+
+   // restore statistics and entries modified by SetBinContent
+   hnew->SetEntries(entries);
+   if (!resetStat) hnew->PutStats(stat);
 
    delete [] oldBins;
    if (oldErrors) delete [] oldErrors;

hist/hist/src/TProfile2D.cxx

@@ -1522,9 +1522,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.
@@ -1540,14 +1542,33 @@ 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 (you should specify `newname`).
+///  The parameter `nxgroup` (`nygroup`) is the number of variable size bins for the x-axis
+///  (y-axis) in the created profile.
+///  The arrays `xbins` and `ybins` must contain `nxgroup+1` and `nygroup+1` elements that
+///  represent the low-edge of the x and y bins respectively.
+///  The data of the old bins are added to the new bin which contains the bin center
+///  of the old bins. It is possible that information from the old binning are attached
+///  to the under-/overflow bins of the new binning.
+///
+///  examples: if hp is an existing TProfile2D with 100 bins on x-axis
+///  and 100 bins y-axis
+///
+/// ~~~ {.cpp}
+///      Double_t xbins[25] = {...} array of low-edges for x-axis (xbins[25] is the upper edge of last bin)
+///      Double_t ybins[25] = {...} array of low-edges for y-axis (ybins[25] is the upper edge of last bin)
+///      hp->Rebin(24,24,"hpnew",xbins,ybins);  //creates a new variable bin size profile hpnew
+/// ~~~
 
-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();
@@ -1557,23 +1578,46 @@ TProfile2D * TProfile2D::Rebin2D(Int_t nxgroup ,Int_t nygroup,const char * newna
       Double_t ymin  = fYaxis.GetXmin();
       Double_t ymax  = fYaxis.GetXmax();
       if ((nxgroup <= 0) || (nxgroup > nxbins)) {
-         Error("Rebin", "Illegal value of nxgroup=%d",nxgroup);
+         Error("Rebin2D", "Illegal value of nxgroup=%d",nxgroup);
          return nullptr;
       }
       if ((nygroup <= 0) || (nygroup > nybins)) {
-         Error("Rebin", "Illegal value of nygroup=%d",nygroup);
+         Error("Rebin2D", "Illegal value of nygroup=%d",nygroup);
          return nullptr;
       }
 
       Int_t newxbins = nxbins/nxgroup;
-      Int_t newybins = nybins/nygroup;
+      if (!xbins) {
+         Int_t nbg = nxbins/nxgroup;
+         //warning if bins are added to the overflow bin
+         if (nbg*nxgroup != nxbins) {
+            Warning("Rebin2D", "nxgroup=%d should be an exact divider of nxbins=%d",nxgroup,nxbins);
+         }
+      }
+      else {
+         // in the case of xbins given (rebinning in variable bins) ngroup is the new number of bins.
+         // and number of grouped bins is not constant.
+         // when looping for setting the contents for the new histogram we
+         // need to loop on all bins of original histogram. Set then nxgroup=nxbins
+         newxbins = nxgroup;
+         nxgroup = nxbins;
+      }
 
-      //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);
+      Int_t newybins = nybins/nygroup;
+      if (!ybins) {
+         Int_t nbg = nybins/nygroup;
+         //warning if bins are added to the overflow bin
+         if (nbg*nygroup != nybins) {
+            Warning("Rebin2D", "nygroup=%d should be an exact divider of nybins=%d",nygroup,nybins);
+         }
       }
-      if(newybins*nygroup != nybins) {
-         Warning("Rebin", "nygroup=%d should be an exact divider of nybins=%d",nygroup,nybins);
+      else {
+         // in the case of ybins given (rebinning in variable bins) ngroup is the new number of bins.
+         // and number of grouped bins is not constant.
+         // when looping for setting the contents for the new histogram we
+         // need to loop on all bins of original histogram. Set then nygroup=nybins
+         newybins = nygroup;
+         nygroup = nybins;
       }
 
       //save old bin contents in new arrays
@@ -1594,35 +1638,38 @@ TProfile2D * TProfile2D::Rebin2D(Int_t nxgroup ,Int_t nygroup,const char * newna
 
       // create a clone of the old profile if newname is specified
       TProfile2D *hnew = this;
-      if(newname && strlen(newname) > 0) {
+      if((newname && strlen(newname) > 0) || xbins || ybins) {
          hnew = (TProfile2D*)Clone(newname);
       }
 
       // 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) {
+      if(!xbins && (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) {
+      if(!ybins && (newybins*nygroup != nybins)) {
          ymax = fYaxis.GetBinUpEdge(newybins*nygroup);
          hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin
       }
 
       //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];
+      if((!xbins && (fXaxis.GetXbins()->GetSize() > 0)) || (!ybins && (fYaxis.GetXbins()->GetSize() > 0))){
+         // for rebinning of variable bins in a constant group
+         Double_t* xbinsTmp = new Double_t[newxbins+1];
+         Double_t* ybinsTmp = new Double_t[newybins+1];
          for(Int_t i=0; i < newxbins+1; i++)
-            xbins[i] = fXaxis.GetBinLowEdge(1+i*nxgroup);
+            xbinsTmp[i] = fXaxis.GetBinLowEdge(1+i*nxgroup);
          for(Int_t j=0; j < newybins+1; j++)
-            ybins[j] = fYaxis.GetBinLowEdge(1+j*nygroup);
+            ybinsTmp[j] = fYaxis.GetBinLowEdge(1+j*nygroup);
          hnew->SetBins(newxbins,xbins,newybins,ybins);
-         delete [] xbins;
-         delete [] ybins;
-      }
+         delete [] xbinsTmp;
+         delete [] ybinsTmp;
+      // when rebinning in variable bins
+      } else if (xbins && ybins) {
+          hnew->SetBins(newxbins,xbins,newybins,ybins);
       //fixed bin size
-      else{
+      } else{
          hnew->SetBins(newxbins,xmin,xmax,newybins,ymin,ymax);
       }