rebin.C

TH2D *Rebin(TH2 *h, Int_t n = 1, Int_t m = 1)
{
    Int_t  nxbins = h->GetXaxis()->GetNbins();
    Axis_t xmin   = h->GetXaxis()->GetXmin();
    Axis_t xmax   = h->GetXaxis()->GetXmax();
    Axis_t dx     = (xmax-xmin)/nxbins;
    Int_t  nxbso  = nxbins;
    nxbins /= n;
    dx     *= n;
    xmax = xmin + nxbins*dx;

    Int_t  nybins = h->GetYaxis()->GetNbins();
    Axis_t ymin   = h->GetYaxis()->GetXmin();
    Axis_t ymax   = h->GetYaxis()->GetXmax();
    Axis_t dy     = (ymax-ymin)/nybins;
    Int_t  nybso  = nybins;
    nybins /= m;
    dy     *= m;
    ymax = ymin + nybins*dy;

    Char_t newname[256];
    sprintf(newname,"%s_chbn",h->GetName());
    TH2D *hnew = new TH2D(newname,h->GetTitle(),nxbins,xmin,xmax,
                                                nybins,ymin,ymax);

    Stat_t z;
    Int_t i, j;
    for (i=0; i<=nxbins+1; i++) {
        Int_t ilo, ihi;
        if (i > 0 && i <= nxbins) {
            ilo = n*(i-1) + 1;
            ihi = n*i;
        } else if (i == 0) {
            ilo = ihi = 0;
        } else {
            ilo = n*nxbins + 1;
            ihi = nxbso;
        }
        for (j=0; j<=nybins+1; j++) {
            Int_t jlo, jhi;
            if (j > 0 && j <= nybins) {
                jlo = m*(j-1) + 1;
                jhi = m*j;
            } else if (j == 0) {
                jlo = jhi = 0;
            } else {
                jlo = m*nybins + 1;
                jhi = nybso;
            }
            z  = h->Integral(ilo,ihi,jlo,jhi);
            if (z > 0) {
                hnew->Fill(hnew->GetXaxis()->GetBinCenter(i),
                           hnew->GetYaxis()->GetBinCenter(j),z);
            }
        }       
    }
    return hnew;
}