    // @(#)root/hist:$Name:  $:$Id: TH2.cxx,v 1.76 2005/06/24 20:44:07 brun Exp $
    // Author: Rene Brun   26/12/94

    /*************************************************************************
     * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
     * All rights reserved.                                                  *
     *                                                                       *
     * For the licensing terms see $ROOTSYS/LICENSE.                         *
     * For the list of contributors see $ROOTSYS/README/CREDITS.             *
     *************************************************************************/

    #include "TROOT.h"
    #include "TH2.h"
    #include "TVirtualPad.h"
    #include "TF2.h"
    #include "TProfile.h"
    #include "TRandom.h"
    #include "TMatrixFBase.h"
    #include "TMatrixDBase.h"
    #include "THLimitsFinder.h"

    ClassImp(TH2)

    //______________________________________________________________________________
    //
    // Service class for 2-Dim histogram classes
    //
    //  TH2C a 2-D histogram with one byte per cell (char)
    //  TH2S a 2-D histogram with two bytes per cell (short integer)
    //  TH2I a 2-D histogram with four bytes per cell (32 bits integer)
    //  TH2F a 2-D histogram with four bytes per cell (float)
    //  TH2D a 2-D histogram with eight bytes per cell (double)
    //

    //______________________________________________________________________________
    TH2::TH2()
    {
       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       SetBinsLength(9);
    }

    //______________________________________________________________________________
    TH2::TH2(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup
                                         ,Int_t nbinsy,Axis_t ylow,Axis_t yup)
         :TH1(name,title,nbinsx,xlow,xup)
    {
       // see comments in the TH1 base class constructors

       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       if (nbinsy <= 0) nbinsy = 1;
       fYaxis.Set(nbinsy,ylow,yup);
       fNcells      = fNcells*(nbinsy+2); // fNCells is set in the TH1 constructor
    }

    //______________________________________________________________________________
    TH2::TH2(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
                                         ,Int_t nbinsy,Axis_t ylow,Axis_t yup)
         :TH1(name,title,nbinsx,xbins)
    {
       // see comments in the TH1 base class constructors
       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       if (nbinsy <= 0) nbinsy = 1;
       fYaxis.Set(nbinsy,ylow,yup);
       fNcells      = fNcells*(nbinsy+2); // fNCells is set in the TH1 constructor
    }

    //______________________________________________________________________________
    TH2::TH2(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup
                                         ,Int_t nbinsy,const Double_t *ybins)
         :TH1(name,title,nbinsx,xlow,xup)
    {
       // see comments in the TH1 base class constructors
       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       if (nbinsy <= 0) nbinsy = 1;
       if (ybins) fYaxis.Set(nbinsy,ybins);
       else       fYaxis.Set(nbinsy,0,1);
       fNcells      = fNcells*(nbinsy+2); // fNCells is set in the TH1 constructor
    }

    //______________________________________________________________________________
    TH2::TH2(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
                                               ,Int_t nbinsy,const Double_t *ybins)
         :TH1(name,title,nbinsx,xbins)
    {
       // see comments in the TH1 base class constructors
       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       if (nbinsy <= 0) nbinsy = 1;
       if (ybins) fYaxis.Set(nbinsy,ybins);
       else       fYaxis.Set(nbinsy,0,1);
       fNcells      = fNcells*(nbinsy+2); // fNCells is set in the TH1 constructor
    }

    //______________________________________________________________________________
    TH2::TH2(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins
                                               ,Int_t nbinsy,const Float_t *ybins)
         :TH1(name,title,nbinsx,xbins)
    {
       // see comments in the TH1 base class constructors
       fDimension   = 2;
       fScalefactor = 1;
       fTsumwy      = fTsumwy2 = fTsumwxy = 0;
       if (nbinsy <= 0) nbinsy = 1;
       if (ybins) fYaxis.Set(nbinsy,ybins);
       else       fYaxis.Set(nbinsy,0,1);
       fNcells      = fNcells*(nbinsy+2); // fNCells is set in the TH1 constructor.
    }

    //______________________________________________________________________________
    TH2::TH2(const TH2 &h) : TH1()
    {
       // Copy constructor.
       // The list of functions is not copied. (Use Clone if needed)
       Copy((TObject&)h);
    }

    //______________________________________________________________________________
    TH2::~TH2()
    {
    }

    //______________________________________________________________________________
    Int_t TH2::BufferEmpty(Int_t action)
    {
    // Fill histogram with all entries in the buffer.
    // action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
    // action =  0 histogram is filled from the buffer
    // action =  1 histogram is filled and buffer is deleted
    //             The buffer is automatically deleted when the number of entries
    //             in the buffer is greater than the number of entries in the histogram

       // do we need to compute the bin size?
       if (!fBuffer) return 0;
       Int_t nbentries = (Int_t)fBuffer[0];
       if (!nbentries) return 0;
       Double_t *buffer = fBuffer;
       if (nbentries < 0) {
          if (action == 0) return 0;
          nbentries  = -nbentries;
          fBuffer=0;
          Reset();
          fBuffer = buffer;
       }
       if (TestBit(kCanRebin) || fXaxis.GetXmax() <= fXaxis.GetXmin() || fYaxis.GetXmax() <= fYaxis.GetXmin()) {
          //find min, max of entries in buffer
          Double_t xmin = fBuffer[2];
          Double_t xmax = xmin;
          Double_t ymin = fBuffer[3];
          Double_t ymax = ymin;
          for (Int_t i=1;i<nbentries;i++) {
             Double_t x = fBuffer[3*i+2];
             if (x < xmin) xmin = x;
             if (x > xmax) xmax = x;
             Double_t y = fBuffer[3*i+3];
             if (y < ymin) ymin = y;
             if (y > ymax) ymax = y;
          }
          if (fXaxis.GetXmax() <= fXaxis.GetXmin() || fYaxis.GetXmax() <= fYaxis.GetXmin()) {
             THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax,ymin,ymax);
          } else {
             fBuffer = 0;
             Int_t keep = fBufferSize; fBufferSize = 0;
             if (xmin <  fXaxis.GetXmin()) RebinAxis(xmin,"X");
             if (xmax >= fXaxis.GetXmax()) RebinAxis(xmax,"X");
             if (ymin <  fYaxis.GetXmin()) RebinAxis(ymin,"Y");
             if (ymax >= fYaxis.GetXmax()) RebinAxis(ymax,"Y");
             fBuffer = buffer;
             fBufferSize = keep;
          }
       }

       fBuffer = 0;
       for (Int_t i=0;i<nbentries;i++) {
          Fill(buffer[3*i+2],buffer[3*i+3],buffer[3*i+1]);
       }
       fBuffer = buffer;

       if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
       else {
          if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
          else                              fBuffer[0] = 0;
       }
       return nbentries;
    }

    //______________________________________________________________________________
    Int_t TH2::BufferFill(Axis_t x, Axis_t y, Stat_t w)
    {
    // accumulate arguments in buffer. When buffer is full, empty the buffer
    // fBuffer[0] = number of entries in buffer
    // fBuffer[1] = w of first entry
    // fBuffer[2] = x of first entry
    // fBuffer[3] = y of first entry

       if (!fBuffer) return -3;
       Int_t nbentries = (Int_t)fBuffer[0];
       if (nbentries < 0) {
          nbentries  = -nbentries;
          fBuffer[0] =  nbentries;
          if (fEntries > 0) {
             Double_t *buffer = fBuffer; fBuffer=0;
             Reset();
             fBuffer = buffer;
          }
       }
       if (3*nbentries+3 >= fBufferSize) {
          BufferEmpty(1);
          return Fill(x,y,w);
       }
       fBuffer[3*nbentries+1] = w;
       fBuffer[3*nbentries+2] = x;
       fBuffer[3*nbentries+3] = y;
       fBuffer[0] += 1;
       return -3;
    }
    Double_t TH2::Chi2Test(TH1 *h, Option_t *option, Int_t constraint)
    {
      //The Chi2 (Pearson's) test for differences between h and this histogram.
      //a small value of prob indicates a significant difference between the distributions
      //
      //if the data was collected in such a way that the number of entries
      //in the first histogram is necessarily equal to the number of entries
      //in the second, the parameter _constraint_ must be made 1. Default is 0.
      //any additional constraints on the data lower the number of degrees of freedom
      //(i.e. increase constraint to more positive values) in accordance with
      //their number
      //
      ///options:
      //"O" -overflows included
      //"U" - underflows included
      //
      //"P" - print information about number of degrees of freedom and
      //the value of chi2
      //by default underflows and overflows are not included

      //algorithm taken from "Numerical Recipes in C++"
      // implementation by Anna Kreshuk

      Int_t df;
      Double_t chsq = 0;
      Double_t prob;
      Double_t temp;
      Double_t koef1,koef2;
      Double_t nen1, nen2;
      Double_t bin1, bin2;
      Int_t i_start, i_end, j_start, j_end;

      TString opt = option;
      opt.ToUpper();

      TAxis *xaxis1 = this->GetXaxis();
      TAxis *xaxis2 = h->GetXaxis();
      TAxis *yaxis1 = this->GetYaxis();
      TAxis *yaxis2 = h->GetYaxis();

      Int_t nbinx1   = xaxis1->GetNbins();
      Int_t nbinx2   = xaxis2->GetNbins();
      Int_t nbiny1   = yaxis1->GetNbins();
      Int_t nbiny2   = yaxis2->GetNbins();

      //check dimensions
      if (this->GetDimension() != 2 || h->GetDimension() != 2) {
         Error("Chi2Test","Histograms must be 2-D");
         return 0;
      }
      //check that the histograms are not empty
      nen1 = this->GetEntries();
      nen2 = h->GetEntries();
      if((nen1==0)||(nen2==0)){
         Error("Chi2Test","one of the histograms is empty");
         return 0;
      }
      //check channels
      if (nbinx1 != nbinx2) {
         Error("Chi2Test", "different number of x channels 1 - %d and 2 - %d", nbinx1, nbinx2);
         return 0;
      }
      if (nbiny1!=nbiny2) {
         Error("Chi2Test", "different number of y channels 1 - %d and 2 - %d", nbiny1, nbiny2);
         return 0;
      }

      //check binning

      Double_t difprec = 1e-5;
      Double_t diff1 = TMath::Abs(xaxis1->GetXmin() - xaxis2->GetXmin());
      Double_t diff2 = TMath::Abs(xaxis1->GetXmax() - xaxis2->GetXmax());
      if (diff1 > difprec || diff2 > difprec) {
         Error("Chi2Test","histograms with different binning along X");
         return 0;
      }
      diff1 = TMath::Abs(yaxis1->GetXmin() - yaxis2->GetXmin());
      diff2 = TMath::Abs(yaxis1->GetXmax() - yaxis2->GetXmax());
      if (diff1 > difprec || diff2 > difprec) {
         Error("Chi2Test","histograms with different binning along Y");
         return 0;
      }

      //check options
      i_start=1; j_start=1;
      i_end = nbinx1; j_end = nbiny1;
      df = nbinx1 * nbiny1-constraint; //total number of bins
      if (opt.Contains("U")) {
         i_start = 0;
         j_start = 0;
         df      +=nbinx1 + nbiny1;
      }
      if (opt.Contains("O")) {
         i_end = nbinx1+1;
         j_end = nbiny1+1;
         df   += nbinx1 + nbiny1;
      }

      //the test

      if (TMath::Abs(nen1-nen2) > difprec){
         koef1 = TMath::Sqrt(nen2/nen1);
         koef2 = TMath::Sqrt(nen1/nen2);
      }
      else{
         koef1 = 1;
         koef2 = 1;
      }

      for (Int_t i=i_start; i<i_end; i++){
         for (Int_t j=j_start; j<j_end; j++){
            bin1 = this->GetCellContent(i, j);
            bin2 = h->GetCellContent(i, j);
            if (bin1==0 && bin2==0){
               --df; //no data means one less degree of freedom
            } else{
               temp  = koef1*bin1 - koef2 *bin2;
               chsq += temp*temp/(bin1+bin2);
            }
         }
      }

       prob = TMath::Prob(0.5*chsq, Int_t(0.5*df));

      if (opt.Contains("P")){
         Printf("the value of chi2 = %f\n", chsq);
         Printf("the number of degrees of freedom = %d\n", df);

      }

      return prob;
    }

    //______________________________________________________________________________
    void TH2::Copy(TObject &obj) const
    {
       TH1::Copy(obj);
       ((TH2&)obj).fScalefactor = fScalefactor;
       ((TH2&)obj).fTsumwy      = fTsumwy;
       ((TH2&)obj).fTsumwy2     = fTsumwy2;
       ((TH2&)obj).fTsumwxy     = fTsumwxy;
    }

    //______________________________________________________________________________
    Int_t TH2::Fill(Axis_t x,Axis_t y)
    {
    //*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y by 1*-*-*-*-*-*-*-*-*-*
    //*-*                  ==================================
    //*-*
    //*-* if x or/and y is less than the low-edge of the corresponding axis first bin,
    //*-*   the Underflow cell is incremented.
    //*-* if x or/and y is greater than the upper edge of corresponding axis last bin,
    //*-*   the Overflow cell is incremented.
    //*-*
    //*-* If the storage of the sum of squares of weights has been triggered,
    //*-* via the function Sumw2, then the sum of the squares of weights is incremented
    //*-* by 1in the cell corresponding to x,y.
    //*-*
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

       if (fBuffer) return BufferFill(x,y,1);

       Int_t binx, biny, bin;
       fEntries++;
       binx = fXaxis.FindBin(x);
       biny = fYaxis.FindBin(y);
       bin  = biny*(fXaxis.GetNbins()+2) + binx;
       AddBinContent(bin);
       if (fSumw2.fN) ++fSumw2.fArray[bin];
       if (binx == 0 || binx > fXaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       if (biny == 0 || biny > fYaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       ++fTsumw;
       ++fTsumw2;
       fTsumwx  += x;
       fTsumwx2 += x*x;
       fTsumwy  += y;
       fTsumwy2 += y*y;
       fTsumwxy += x*y;
       return bin;
    }

    //______________________________________________________________________________
    Int_t TH2::Fill(Axis_t x, Axis_t y, Stat_t w)
    {
    //*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y by a weight w*-*-*-*-*-*
    //*-*                  ===========================================
    //*-*
    //*-* if x or/and y is less than the low-edge of the corresponding axis first bin,
    //*-*   the Underflow cell is incremented.
    //*-* if x or/and y is greater than the upper edge of corresponding axis last bin,
    //*-*   the Overflow cell is incremented.
    //*-*
    //*-* If the storage of the sum of squares of weights has been triggered,
    //*-* via the function Sumw2, then the sum of the squares of weights is incremented
    //*-* by w^2 in the cell corresponding to x,y.
    //*-*
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

       if (fBuffer) return BufferFill(x,y,w);

       Int_t binx, biny, bin;
       fEntries++;
       binx = fXaxis.FindBin(x);
       biny = fYaxis.FindBin(y);
       bin  = biny*(fXaxis.GetNbins()+2) + binx;
       AddBinContent(bin,w);
       if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
       if (binx == 0 || binx > fXaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       if (biny == 0 || biny > fYaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       Stat_t z= (w > 0 ? w : -w);
       fTsumw   += z;
       fTsumw2  += z*z;
       fTsumwx  += z*x;
       fTsumwx2 += z*x*x;
       fTsumwy  += z*y;
       fTsumwy2 += z*y*y;
       fTsumwxy += z*x*y;
       return bin;
    }

    //______________________________________________________________________________
    Int_t TH2::Fill(const char *namex, const char *namey, Stat_t w)
    {
    // Increment cell defined by namex,namey by a weight w
    //
    // if x or/and y is less than the low-edge of the corresponding axis first bin,
    //   the Underflow cell is incremented.
    // if x or/and y is greater than the upper edge of corresponding axis last bin,
    //   the Overflow cell is incremented.
    //
    // If the storage of the sum of squares of weights has been triggered,
    // via the function Sumw2, then the sum of the squares of weights is incremented
    // by w^2 in the cell corresponding to x,y.
    //

       Int_t binx, biny, bin;
       fEntries++;
       binx = fXaxis.FindBin(namex);
       biny = fYaxis.FindBin(namey);
       bin  = biny*(fXaxis.GetNbins()+2) + binx;
       AddBinContent(bin,w);
       if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
       if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
       if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
       Axis_t x = fXaxis.GetBinCenter(binx);
       Axis_t y = fYaxis.GetBinCenter(biny);
       Stat_t z= (w > 0 ? w : -w);
       fTsumw   += z;
       fTsumw2  += z*z;
       fTsumwx  += z*x;
       fTsumwx2 += z*x*x;
       fTsumwy  += z*y;
       fTsumwy2 += z*y*y;
       fTsumwxy += z*x*y;
       return bin;
    }

    //______________________________________________________________________________
    Int_t TH2::Fill(const char *namex, Axis_t y, Stat_t w)
    {
    // Increment cell defined by namex,y by a weight w
    //
    // if x or/and y is less than the low-edge of the corresponding axis first bin,
    //   the Underflow cell is incremented.
    // if x or/and y is greater than the upper edge of corresponding axis last bin,
    //   the Overflow cell is incremented.
    //
    // If the storage of the sum of squares of weights has been triggered,
    // via the function Sumw2, then the sum of the squares of weights is incremented
    // by w^2 in the cell corresponding to x,y.
    //

       Int_t binx, biny, bin;
       fEntries++;
       binx = fXaxis.FindBin(namex);
       biny = fYaxis.FindBin(y);
       bin  = biny*(fXaxis.GetNbins()+2) + binx;
       AddBinContent(bin,w);
       if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
       if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
       if (biny == 0 || biny > fYaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       Axis_t x = fXaxis.GetBinCenter(binx);
       Stat_t z= (w > 0 ? w : -w);
       fTsumw   += z;
       fTsumw2  += z*z;
       fTsumwx  += z*x;
       fTsumwx2 += z*x*x;
       fTsumwy  += z*y;
       fTsumwy2 += z*y*y;
       fTsumwxy += z*x*y;
       return bin;
    }

    //______________________________________________________________________________
    Int_t TH2::Fill(Axis_t x, const char *namey, Stat_t w)
    {
    // Increment cell defined by x,namey by a weight w
    //
    // if x or/and y is less than the low-edge of the corresponding axis first bin,
    //   the Underflow cell is incremented.
    // if x or/and y is greater than the upper edge of corresponding axis last bin,
    //   the Overflow cell is incremented.
    //
    // If the storage of the sum of squares of weights has been triggered,
    // via the function Sumw2, then the sum of the squares of weights is incremented
    // by w^2 in the cell corresponding to x,y.
    //

       Int_t binx, biny, bin;
       fEntries++;
       binx = fXaxis.FindBin(x);
       biny = fYaxis.FindBin(namey);
       bin  = biny*(fXaxis.GetNbins()+2) + binx;
       AddBinContent(bin,w);
       if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
       if (binx == 0 || binx > fXaxis.GetNbins()) {
          if (!fgStatOverflows) return -1;
       }
       if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
       Axis_t y = fYaxis.GetBinCenter(biny);
       Stat_t z= (w > 0 ? w : -w);
       fTsumw   += z;
       fTsumw2  += z*z;
       fTsumwx  += z*x;
       fTsumwx2 += z*x*x;
       fTsumwy  += z*y;
       fTsumwy2 += z*y*y;
       fTsumwxy += z*x*y;
       return bin;
    }

    //______________________________________________________________________________
    void TH2::FillN(Int_t ntimes, const Axis_t *x, const Axis_t *y, const Double_t *w, Int_t stride)
    {
    //*-*-*-*-*-*-*Fill a 2-D histogram with an array of values and weights*-*-*-*
    //*-*          ========================================================
    //*-*
    //*-* ntimes:  number of entries in arrays x and w (array size must be ntimes*stride)
    //*-* x:       array of x values to be histogrammed
    //*-* y:       array of y values to be histogrammed
    //*-* w:       array of weights
    //*-* stride:  step size through arrays x, y and w
    //*-*
    //*-* If the storage of the sum of squares of weights has been triggered,
    //*-* via the function Sumw2, then the sum of the squares of weights is incremented
    //*-* by w[i]^2 in the cell corresponding to x[i],y[i].
    //*-* if w is NULL each entry is assumed a weight=1
    //*-*
    //*-* NB: function only valid for a TH2x object
    //*-*
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
       Int_t binx, biny, bin, i;
       fEntries += ntimes;
       Double_t ww = 1;
       ntimes *= stride;
       for (i=0;i<ntimes;i+=stride) {
          binx = fXaxis.FindBin(x[i]);
          biny = fYaxis.FindBin(y[i]);
          bin  = biny*(fXaxis.GetNbins()+2) + binx;
          if (w) ww = w[i];
          AddBinContent(bin,ww);
          if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
          if (binx == 0 || binx > fXaxis.GetNbins()) {
             if (!fgStatOverflows) continue;
          }
          if (biny == 0 || biny > fYaxis.GetNbins()) {
             if (!fgStatOverflows) continue;
          }
          Stat_t z= (ww > 0 ? ww : -ww);
          fTsumw   += z;
          fTsumw2  += z*z;
          fTsumwx  += z*x[i];
          fTsumwx2 += z*x[i]*x[i];
          fTsumwy  += z*y[i];
          fTsumwy2 += z*y[i]*y[i];
          fTsumwxy += z*x[i]*y[i];
       }
    }

    //______________________________________________________________________________
    void TH2::FillRandom(const char *fname, Int_t ntimes)
    {
    //*-*-*-*-*-*-*Fill histogram following distribution in function fname*-*-*-*
    //*-*          =======================================================
    //*-*
    //*-*   The distribution contained in the function fname (TF2) is integrated
    //*-*   over the channel contents.
    //*-*   It is normalized to 1.
    //*-*   Getting one random number implies:
    //*-*     - Generating a random number between 0 and 1 (say r1)
    //*-*     - Look in which bin in the normalized integral r1 corresponds to
    //*-*     - Fill histogram channel
    //*-*   ntimes random numbers are generated
    //*-*
    //*-*  One can also call TF2::GetRandom2 to get a random variate from a function.
    //*-*
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*

       Int_t bin, binx, biny, ibin, loop;
       Double_t r1, x, y, xv[2];
    //*-*- Search for fname in the list of ROOT defined functions
       TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
       if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }

    //*-*- Allocate temporary space to store the integral and compute integral
       Int_t nbinsx = GetNbinsX();
       Int_t nbinsy = GetNbinsY();
       Int_t nbins  = nbinsx*nbinsy;

       Double_t *integral = new Double_t[nbins+1];
       ibin = 0;
       integral[ibin] = 0;
       for (biny=1;biny<=nbinsy;biny++) {
          xv[1] = fYaxis.GetBinCenter(biny);
          for (binx=1;binx<=nbinsx;binx++) {
             xv[0] = fXaxis.GetBinCenter(binx);
             ibin++;
             integral[ibin] = integral[ibin-1] + f1->Eval(xv[0],xv[1]);
          }
       }

    //*-*- Normalize integral to 1
       if (integral[nbins] == 0 ) {
          Error("FillRandom", "Integral = zero"); return;
       }
       for (bin=1;bin<=nbins;bin++)  integral[bin] /= integral[nbins];

    //*-*--------------Start main loop ntimes
       for (loop=0;loop<ntimes;loop++) {
          r1 = gRandom->Rndm(loop);
          ibin = TMath::BinarySearch(nbins,&integral[0],r1);
          biny = ibin/nbinsx;
          binx = 1 + ibin - nbinsx*biny;
          biny++;
          x    = fXaxis.GetBinCenter(binx);
          y    = fYaxis.GetBinCenter(biny);
          Fill(x,y, 1.);
      }
      delete [] integral;
    }

    //______________________________________________________________________________
    void TH2::FillRandom(TH1 *h, Int_t ntimes)
    {
    //*-*-*-*-*-*-*Fill histogram following distribution in histogram h*-*-*-*
    //*-*          ====================================================
    //*-*
    //*-*   The distribution contained in the histogram h (TH2) is integrated
    //*-*   over the channel contents.
    //*-*   It is normalized to 1.
    //*-*   Getting one random number implies:
    //*-*     - Generating a random number between 0 and 1 (say r1)
    //*-*     - Look in which bin in the normalized integral r1 corresponds to
    //*-*     - Fill histogram channel
    //*-*   ntimes random numbers are generated
    //*-*
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*

       if (!h) { Error("FillRandom", "Null histogram"); return; }
       if (fDimension != h->GetDimension()) {
          Error("FillRandom", "Histograms with different dimensions"); return;
       }

       if (h->ComputeIntegral() == 0) return;

       Int_t loop;
       Axis_t x,y;
       TH2 *h2 = (TH2*)h;
       for (loop=0;loop<ntimes;loop++) {
          h2->GetRandom2(x,y);
          Fill(x,y,1.);
       }
    }

    //______________________________________________________________________________
    void TH2::FitSlicesX(TF1 *f1, Int_t binmin, Int_t binmax, Int_t cut, Option_t *option)
    {
    // Project slices along X in case of a 2-D histogram, then fit each slice
    // with function f1 and make a histogram for each fit parameter
    // Only bins along Y between binmin and binmax are considered.
    // if f1=0, a gaussian is assumed
    // Before invoking this function, one can set a subrange to be fitted along X
    // via f1->SetRange(xmin,xmax)
    // The argument option (default="QNR") can be used to change the fit options.
    //     "Q"  means Quiet mode
    //     "N"  means do not show the result of the fit
    //     "R"  means fit the function in the specified function range
    //     "G2" merge 2 consecutive bins along X
    //     "G3" merge 3 consecutive bins along X
    //     "G4" merge 4 consecutive bins along X
    //     "G5" merge 5 consecutive bins along X
    //
    // Note that the generated histograms are added to the list of objects
    // in the current directory. It is the user's responsability to delete
    // these histograms.
    //
    //  Example: Assume a 2-d histogram h2
    //   Root > h2->FitSlicesX(); produces 4 TH1D histograms
    //          with h2_0 containing parameter 0(Constant) for a Gaus fit
    //                    of each bin in Y projected along X
    //          with h2_1 containing parameter 1(Mean) for a gaus fit
    //          with h2_2 containing parameter 2(RMS)  for a gaus fit
    //          with h2_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
    //
    //   Root > h2->FitSlicesX(0,15,22,10);
    //          same as above, but only for bins 15 to 22 along Y
    //          and only for bins in Y for which the corresponding projection
    //          along X has more than cut bins filled.
    //
    //  NOTE: To access the generated histograms in the current directory, do eg:
    //     TH1D *h2_1 = (TH1D*)gDirectory->Get("h2_1");

       Int_t nbins  = fYaxis.GetNbins();
       if (binmin < 1) binmin = 1;
       if (binmax > nbins) binmax = nbins;
       if (binmax < binmin) {binmin = 1; binmax = nbins;}
       TString opt = option;
       opt.ToLower();
       Int_t ngroup = 1;
       if (opt.Contains("g2")) {ngroup = 2; opt.ReplaceAll("g2","");}
       if (opt.Contains("g3")) {ngroup = 3; opt.ReplaceAll("g3","");}
       if (opt.Contains("g4")) {ngroup = 4; opt.ReplaceAll("g4","");}
       if (opt.Contains("g5")) {ngroup = 5; opt.ReplaceAll("g5","");}

       //default is to fit with a gaussian
       if (f1 == 0) {
          f1 = (TF1*)gROOT->GetFunction("gaus");
          if (f1 == 0) f1 = new TF1("gaus","gaus",fXaxis.GetXmin(),fXaxis.GetXmax());
          else         f1->SetRange(fXaxis.GetXmin(),fXaxis.GetXmax());
       }
       Int_t npar = f1->GetNpar();
       if (npar <= 0) return;
       Double_t *parsave = new Double_t[npar];
       f1->GetParameters(parsave);

       //Create one histogram for each function parameter
       Int_t ipar;
       TH1D **hlist = new TH1D*[npar];
       char *name   = new char[2000];
       char *title  = new char[2000];
       const TArrayD *bins = fYaxis.GetXbins();
       for (ipar=0;ipar<npar;ipar++) {
          sprintf(name,"%s_%d",GetName(),ipar);
          sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
          delete gDirectory->FindObject(name);
          if (bins->fN == 0) {
             hlist[ipar] = new TH1D(name,title, nbins, fYaxis.GetXmin(), fYaxis.GetXmax());
          } else {
             hlist[ipar] = new TH1D(name,title, nbins,bins->fArray);
          }
          hlist[ipar]->GetXaxis()->SetTitle(fYaxis.GetTitle());
       }
       sprintf(name,"%s_chi2",GetName());
       delete gDirectory->FindObject(name);
       TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fYaxis.GetXmin(), fYaxis.GetXmax());
       hchi2->GetXaxis()->SetTitle(fYaxis.GetTitle());

       //Loop on all bins in Y, generate a projection along X
       Int_t bin;
       Int_t nentries;
       for (bin=binmin;bin<=binmax;bin += ngroup) {
          TH1D *hpx = ProjectionX("_temp",bin,bin+ngroup-1,"e");
          if (hpx == 0) continue;
          nentries = Int_t(hpx->GetEntries());
          if (nentries == 0 || nentries < cut) {delete hpx; continue;}
          f1->SetParameters(parsave);
          hpx->Fit(f1,opt.Data());
          Int_t npfits = f1->GetNumberFitPoints();
          if (npfits > npar && npfits >= cut) {
             Int_t binx = bin + ngroup/2;
             for (ipar=0;ipar<npar;ipar++) {
                hlist[ipar]->Fill(fYaxis.GetBinCenter(binx),f1->GetParameter(ipar));
                hlist[ipar]->SetBinError(binx,f1->GetParError(ipar));
             }
             hchi2->Fill(fYaxis.GetBinCenter(binx),f1->GetChisquare()/(npfits-npar));
          }
          delete hpx;
       }
       delete [] parsave;
       delete [] name;
       delete [] title;
       delete [] hlist;
    }

    //______________________________________________________________________________
    void TH2::FitSlicesY(TF1 *f1, Int_t binmin, Int_t binmax, Int_t cut, Option_t *option)
    {
    // Project slices along Y in case of a 2-D histogram, then fit each slice
    // with function f1 and make a histogram for each fit parameter
    // Only bins along X between binmin and binmax are considered.
    // if f1=0, a gaussian is assumed
    // Before invoking this function, one can set a subrange to be fitted along Y
    // via f1->SetRange(ymin,ymax)
    // The argument option (default="QNR") can be used to change the fit options.
    //     "Q"  means Quiet mode
    //     "N"  means do not show the result of the fit
    //     "R"  means fit the function in the specified function range
    //     "G2" merge 2 consecutive bins along Y
    //     "G3" merge 3 consecutive bins along Y
    //     "G4" merge 4 consecutive bins along Y
    //     "G5" merge 5 consecutive bins along Y
    //
    // Note that the generated histograms are added to the list of objects
    // in the current directory. It is the user's responsability to delete
    // these histograms.
    //
    //  Example: Assume a 2-d histogram h2
    //   Root > h2->FitSlicesY(); produces 4 TH1D histograms
    //          with h2_0 containing parameter 0(Constant) for a Gaus fit
    //                    of each bin in X projected along Y
    //          with h2_1 containing parameter 1(Mean) for a gaus fit
    //          with h2_2 containing parameter 2(RMS)  for a gaus fit
    //          with h2_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
    //
    //   Root > h2->FitSlicesY(0,15,22,10);
    //          same as above, but only for bins 15 to 22 along X
    //          and only for bins in X for which the corresponding projection
    //          along Y has more than cut bins filled.
    //
    //  NOTE: To access the generated histograms in the current directory, do eg:
    //     TH1D *h2_1 = (TH1D*)gDirectory->Get("h2_1");
    //
    // A complete example of this function is given in

tutorial:fitslicesy.C

    // with the following output:
    //

/* */ //

Int_t nbins = fXaxis.GetNbins(); if (binmin < 1) binmin = 1; if (binmax > nbins) binmax = nbins; if (binmax < binmin) {binmin = 1; binmax = nbins;} TString opt = option; opt.ToLower(); Int_t ngroup = 1; if (opt.Contains("g2")) {ngroup = 2; opt.ReplaceAll("g2","");} if (opt.Contains("g3")) {ngroup = 3; opt.ReplaceAll("g3","");} if (opt.Contains("g4")) {ngroup = 4; opt.ReplaceAll("g4","");} if (opt.Contains("g5")) {ngroup = 5; opt.ReplaceAll("g5","");} //default is to fit with a gaussian if (f1 == 0) { f1 = (TF1*)gROOT->GetFunction("gaus"); if (f1 == 0) f1 = new TF1("gaus","gaus",fYaxis.GetXmin(),fYaxis.GetXmax()); else f1->SetRange(fYaxis.GetXmin(),fYaxis.GetXmax()); } Int_t npar = f1->GetNpar(); if (npar <= 0) return; Double_t *parsave = new Double_t[npar]; f1->GetParameters(parsave); //Create one histogram for each function parameter Int_t ipar; TH1D **hlist = new TH1D*[npar]; char *name = new char[2000]; char *title = new char[2000]; const TArrayD *bins = fXaxis.GetXbins(); for (ipar=0;ipar<npar;ipar++) { sprintf(name,"%s_%d",GetName(),ipar); sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar)); delete gDirectory->FindObject(name); if (bins->fN == 0) { hlist[ipar] = new TH1D(name,title, nbins, fXaxis.GetXmin(), fXaxis.GetXmax()); } else { hlist[ipar] = new TH1D(name,title, nbins,bins->fArray); } hlist[ipar]->GetXaxis()->SetTitle(fXaxis.GetTitle()); } sprintf(name,"%s_chi2",GetName()); delete gDirectory->FindObject(name); TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fXaxis.GetXmin(), fXaxis.GetXmax()); hchi2->GetXaxis()->SetTitle(fXaxis.GetTitle()); //Loop on all bins in X, generate a projection along Y Int_t bin; Int_t nentries; for (bin=binmin;bin<=binmax;bin += ngroup) { TH1D *hpy = ProjectionY("_temp",bin,bin+ngroup-1,"e"); if (hpy == 0) continue; nentries = Int_t(hpy->GetEntries()); if (nentries == 0 || nentries < cut) {delete hpy; continue;} f1->SetParameters(parsave); hpy->Fit(f1,option); Int_t npfits = f1->GetNumberFitPoints(); if (npfits > npar && npfits >= cut) { Int_t biny = bin + ngroup/2; for (ipar=0;ipar<npar;ipar++) { hlist[ipar]->Fill(fXaxis.GetBinCenter(biny),f1->GetParameter(ipar)); hlist[ipar]->SetBinError(biny,f1->GetParError(ipar)); } hchi2->Fill(fXaxis.GetBinCenter(biny),f1->GetChisquare()/(npfits-npar)); } delete hpy; } delete [] parsave; delete [] name; delete [] title; delete [] hlist; } //______________________________________________________________________________ Stat_t TH2::GetCorrelationFactor(Int_t axis1, Int_t axis2) const { //*-*-*-*-*-*-*-*Return correlation factor between axis1 and axis2*-*-*-*-* //*-* ==================================================== if (axis1 < 1 || axis2 < 1 || axis1 > 2 || axis2 > 2) { Error("GetCorrelationFactor","Wrong parameters"); return 0; } if (axis1 == axis2) return 1; Stat_t rms1 = GetRMS(axis1); if (rms1 == 0) return 0; Stat_t rms2 = GetRMS(axis2); if (rms2 == 0) return 0; return GetCovariance(axis1,axis2)/rms1/rms2; } //______________________________________________________________________________ Stat_t TH2::GetCovariance(Int_t axis1, Int_t axis2) const { //*-*-*-*-*-*-*-*Return covariance between axis1 and axis2*-*-*-*-* //*-* ==================================================== if (axis1 < 1 || axis2 < 1 || axis1 > 2 || axis2 > 2) { Error("GetCovariance","Wrong parameters"); return 0; } Stat_t stats[7]; GetStats(stats); Stat_t sumw = stats[0]; //Stat_t sumw2 = stats[1]; Stat_t sumwx = stats[2]; Stat_t sumwx2 = stats[3]; Stat_t sumwy = stats[4]; Stat_t sumwy2 = stats[5]; Stat_t sumwxy = stats[6]; if (sumw == 0) return 0; if (axis1 == 1 && axis2 == 1) { return TMath::Abs(sumwx2/sumw - sumwx/sumw*sumwx/sumw); } if (axis1 == 2 && axis2 == 2) { return TMath::Abs(sumwy2/sumw - sumwy/sumw*sumwy/sumw); } return sumwxy/sumw - sumwx/sumw*sumwy/sumw; } //______________________________________________________________________________ void TH2::GetRandom2(Axis_t &x, Axis_t &y) { // return 2 random numbers along axis x and y distributed according // the cellcontents of a 2-dim histogram Int_t nbinsx = GetNbinsX(); Int_t nbinsy = GetNbinsY(); Int_t nbins = nbinsx*nbinsy; Double_t integral; if (fIntegral) { if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral(); } else { integral = ComputeIntegral(); if (integral == 0 || fIntegral == 0) return; } Float_t r1 = gRandom->Rndm(); Int_t ibin = TMath::BinarySearch(nbins,fIntegral,r1); Int_t biny = ibin/nbinsx; Int_t binx = ibin - nbinsx*biny; x = fXaxis.GetBinLowEdge(binx+1) +fXaxis.GetBinWidth(binx+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]); y = fYaxis.GetBinLowEdge(biny+1) + fYaxis.GetBinWidth(biny+1)*gRandom->Rndm(); } //______________________________________________________________________________ void TH2::GetStats(Stat_t *stats) const { // fill the array stats from the contents of this histogram // The array stats must be correctly dimensionned in the calling program. // stats[0] = sumw // stats[1] = sumw2 // stats[2] = sumwx // stats[3] = sumwx2 // stats[4] = sumwy // stats[5] = sumwy2 // stats[6] = sumwxy if (fBuffer) ((TH2*)this)->BufferEmpty(); Int_t bin, binx, biny; Stat_t w,err; Double_t x,y; if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange) || fYaxis.TestBit(TAxis::kAxisRange)) { for (bin=0;bin<7;bin++) stats[bin] = 0; for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) { y = fYaxis.GetBinCenter(biny); for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) { bin = GetBin(binx,biny); x = fXaxis.GetBinCenter(binx); w = TMath::Abs(GetBinContent(bin)); err = TMath::Abs(GetBinError(bin)); stats[0] += w; stats[1] += err*err; stats[2] += w*x; stats[3] += w*x*x; stats[4] += w*y; stats[5] += w*y*y; stats[6] += w*x*y; } } } else { stats[0] = fTsumw; stats[1] = fTsumw2; stats[2] = fTsumwx; stats[3] = fTsumwx2; stats[4] = fTsumwy; stats[5] = fTsumwy2; stats[6] = fTsumwxy; } } //______________________________________________________________________________ Stat_t TH2::Integral(Option_t *option) const { //Return integral of bin contents. Only bins in the bins range are considered. // By default the integral is computed as the sum of bin contents in the range. // if option "width" is specified, the integral is the sum of // the bin contents multiplied by the bin width in x and in y. return Integral(fXaxis.GetFirst(),fXaxis.GetLast(), fYaxis.GetFirst(),fYaxis.GetLast(),option); } //______________________________________________________________________________ Stat_t TH2::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Option_t *option) const { //Return integral of bin contents in range [binx1,binx2],[biny1,biny2] // for a 2-D histogram // By default the integral is computed as the sum of bin contents in the range. // if option "width" is specified, the integral is the sum of // the bin contents multiplied by the bin width in x and in y. if (fBuffer) ((TH2*)this)->BufferEmpty(); Int_t nbinsx = GetNbinsX(); Int_t nbinsy = GetNbinsY(); if (binx1 < 0) binx1 = 0; if (binx2 > nbinsx+1) binx2 = nbinsx+1; if (binx2 < binx1) binx2 = nbinsx; if (biny1 < 0) biny1 = 0; if (biny2 > nbinsy+1) biny2 = nbinsy+1; if (biny2 < biny1) biny2 = nbinsy; Stat_t integral = 0; //*-*- Loop on bins in specified range TString opt = option; opt.ToLower(); Bool_t width = kFALSE; if (opt.Contains("width")) width = kTRUE; Int_t bin, binx, biny; for (biny=biny1;biny<=biny2;biny++) { for (binx=binx1;binx<=binx2;binx++) { bin = binx +(nbinsx+2)*biny; if (width) integral += GetBinContent(bin)*fXaxis.GetBinWidth(binx)*fYaxis.GetBinWidth(biny); else integral += GetBinContent(bin); } } return integral; } //______________________________________________________________________________ Double_t TH2::KolmogorovTest(TH1 *h2, Option_t *option) const { // Statistical test of compatibility in shape between // THIS histogram and h2, using Kolmogorov test. // Default: Ignore under- and overflow bins in comparison // // option is a character string to specify options // "U" include Underflows in test // "O" include Overflows // "N" include comparison of normalizations // "D" Put out a line of "Debug" printout // // The returned function value is the probability of test // (much less than one means NOT compatible) // // Code adapted by Rene Brun from original HBOOK routine HDIFF TString opt = option; opt.ToUpper(); Double_t prb = 0; TH1 *h1 = (TH1*)this; if (h2 == 0) return 0; TAxis *xaxis1 = h1->GetXaxis(); TAxis *xaxis2 = h2->GetXaxis(); TAxis *yaxis1 = h1->GetYaxis(); TAxis *yaxis2 = h2->GetYaxis(); Int_t ncx1 = xaxis1->GetNbins(); Int_t ncx2 = xaxis2->GetNbins(); Int_t ncy1 = yaxis1->GetNbins(); Int_t ncy2 = yaxis2->GetNbins(); // Check consistency of dimensions if (h1->GetDimension() != 2 || h2->GetDimension() != 2) { Error("KolmogorovTest","Histograms must be 2-D\n"); return 0; } // Check consistency in number of channels if (ncx1 != ncx2) { Error("KolmogorovTest","Number of channels in X is different, %d and %d\n",ncx1,ncx2); return 0; } if (ncy1 != ncy2) { Error("KolmogorovTest","Number of channels in Y is different, %d and %d\n",ncy1,ncy2); return 0; } // Check consistency in channel edges Bool_t afunc1 = kFALSE; Bool_t afunc2 = kFALSE; Double_t difprec = 1e-5; Double_t diff1 = TMath::Abs(xaxis1->GetXmin() - xaxis2->GetXmin()); Double_t diff2 = TMath::Abs(xaxis1->GetXmax() - xaxis2->GetXmax()); if (diff1 > difprec || diff2 > difprec) { Error("KolmogorovTest","histograms with different binning along X"); return 0; } diff1 = TMath::Abs(yaxis1->GetXmin() - yaxis2->GetXmin()); diff2 = TMath::Abs(yaxis1->GetXmax() - yaxis2->GetXmax()); if (diff1 > difprec || diff2 > difprec) { Error("KolmogorovTest","histograms with different binning along Y"); return 0; } // Should we include Uflows, Oflows? Int_t ibeg = 1, jbeg = 1; Int_t iend = ncx1, jend = ncy1; if (opt.Contains("U")) {ibeg = 0; jbeg = 0;} if (opt.Contains("O")) {iend = ncx1+1; jend = ncy1+1;} Int_t i,j; Double_t hsav; Double_t sum1 = 0; Double_t tsum1 = 0; for (i=0;i<=ncx1+1;i++) { for (j=0;j<=ncy1+1;j++) { hsav = h1->GetCellContent(i,j); tsum1 += hsav; if (i >= ibeg && i <= iend && j >= jbeg && j <= jend) sum1 += hsav; } } Double_t sum2 = 0; Double_t tsum2 = 0; for (i=0;i<=ncx1+1;i++) { for (j=0;j<=ncy1+1;j++) { hsav = h2->GetCellContent(i,j); tsum2 += hsav; if (i >= ibeg && i <= iend && j >= jbeg && j <= jend) sum2 += hsav; } } // Check that both scatterplots contain events if (sum1 == 0) { Error("KolmogorovTest","Integral is zero for h1=%s\n",h1->GetName()); return 0; } if (sum2 == 0) { Error("KolmogorovTest","Integral is zero for h2=%s\n",h2->GetName()); return 0; } // Check that scatterplots are not weighted or saturated Double_t num1 = h1->GetEntries(); Double_t num2 = h2->GetEntries(); if (num1 != tsum1) { Warning("KolmogorovTest","Saturation or weighted events for h1=%s, num1=%g, tsum1=%g\n",h1->GetName(),num1,tsum1); } if (num2 != tsum2) { Warning("KolmogorovTest","Saturation or weighted events for h2=%s, num2=%g, tsum2=%g\n",h2->GetName(),num2,tsum2); } // Find first Kolmogorov distance Double_t s1 = 1/sum1; Double_t s2 = 1/sum2; Double_t dfmax = 0; Double_t rsum1=0, rsum2=0; for (i=ibeg;i<=iend;i++) { for (j=jbeg;j<=jend;j++) { rsum1 += s1*h1->GetCellContent(i,j); rsum2 += s2*h2->GetCellContent(i,j); dfmax = TMath::Max(dfmax, TMath::Abs(rsum1-rsum2)); } } // Find second Kolmogorov distance Double_t dfmax2 = 0; rsum1=0, rsum2=0; for (j=jbeg;j<=jend;j++) { for (i=ibeg;i<=iend;i++) { rsum1 += s1*h1->GetCellContent(i,j); rsum2 += s2*h2->GetCellContent(i,j); dfmax2 = TMath::Max(dfmax2, TMath::Abs(rsum1-rsum2)); } } // Get Kolmogorov probability Double_t factnm; if (afunc1) factnm = TMath::Sqrt(sum2); else if (afunc2) factnm = TMath::Sqrt(sum1); else factnm = TMath::Sqrt(sum1*sum2/(sum1+sum2)); Double_t z = dfmax*factnm; Double_t z2 = dfmax2*factnm; prb = TMath::KolmogorovProb(0.5*(z+z2)); Double_t prb1=0, prb2=0; Double_t resum1, resum2, chi2, d12; if (opt.Contains("N")) { //Combine probabilities for shape and normalization, prb1 = prb; resum1 = sum1; if (afunc1) resum1 = 0; resum2 = sum2; if (afunc2) resum2 = 0; d12 = sum1-sum2; chi2 = d12*d12/(resum1+resum2); prb2 = TMath::Prob(chi2,1); // see Eadie et al., section 11.6.2 if (prb > 0 && prb2 > 0) prb = prb*prb2*(1-TMath::Log(prb*prb2)); else prb = 0; } // debug printout if (opt.Contains("D")) { printf(" Kolmo Prob h1 = %s, sum1=%g\n",h1->GetName(),sum1); printf(" Kolmo Prob h2 = %s, sum2=%g\n",h2->GetName(),sum2); printf(" Kolmo Probabil = %f, Max Dist = %g\n",prb,dfmax); if (opt.Contains("N")) printf(" Kolmo Probabil = %f for shape alone, =%f for normalisation alone\n",prb1,prb2); } // This numerical error condition should never occur: if (TMath::Abs(rsum1-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h1=%s\n",h1->GetName()); if (TMath::Abs(rsum2-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h2=%s\n",h2->GetName()); return prb; } //______________________________________________________________________________ Long64_t TH2::Merge(TCollection *list) { //Add all histograms in the collection to this histogram. //This function computes the min/max for the axes, //compute a new number of bins, if necessary, //add bin contents, errors and statistics. //If overflows are present and limits are different the function will fail. //The function returns the total number of entries in the result histogram //if the merge is successfull, -1 otherwise. // //IMPORTANT remark. The 2 axis x and y may have different number //of bins and different limits, BUT the largest bin width must be //a multiple of the smallest bin width and the upper limit must also //be a multiple of the bin width. if (!list) return 0; if (list->IsEmpty()) return (Int_t) GetEntries(); TList inlist; TH1* hclone = (TH1*)Clone("FirstClone"); BufferEmpty(1); // To remove buffer. Reset(); // BufferEmpty sets limits so we can't use it later. SetEntries(0); inlist.Add(hclone); inlist.AddAll(list); TAxis newXAxis; TAxis newYAxis; Bool_t initialLimitsFound = kFALSE; Bool_t same = kTRUE; Bool_t allHaveLimits = kTRUE; TIter next(&inlist); while (TObject *o = next()) { TH2* h = dynamic_cast<TH2*> (o); if (!h) { Error("Add","Attempt to add object of class: %s to a %s", o->ClassName(),this->ClassName()); return -1; } Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax(); allHaveLimits = allHaveLimits && hasLimits; if (hasLimits) { h->BufferEmpty(); if (!initialLimitsFound) { initialLimitsFound = kTRUE; newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax()); newYAxis.Set(h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(), h->GetYaxis()->GetXmax()); } else { if (!RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) { Error("Merge", "Cannot merge histograms - limits are inconsistent:\n " "first: (%d, %f, %f), second: (%d, %f, %f)", newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(), h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax()); } if (!RecomputeAxisLimits(newYAxis, *(h->GetYaxis()))) { Error("Merge", "Cannot merge histograms - limits are inconsistent:\n " "first: (%d, %f, %f), second: (%d, %f, %f)", newYAxis.GetNbins(), newYAxis.GetXmin(), newYAxis.GetXmax(), h->GetYaxis()->GetNbins(), h->GetYaxis()->GetXmin(), h->GetYaxis()->GetXmax()); } } } } next.Reset(); same = same && SameLimitsAndNBins(newXAxis, *GetXaxis()) && SameLimitsAndNBins(newYAxis, *GetYaxis()); if (!same && initialLimitsFound) SetBins(newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(), newYAxis.GetNbins(), newYAxis.GetXmin(), newYAxis.GetXmax()); if (!allHaveLimits) { // fill this histogram with all the data from buffers of histograms without limits while (TH2* h = dynamic_cast<TH2*> (next())) { if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) { // no limits Int_t nbentries = (Int_t)h->fBuffer[0]; for (Int_t i = 0; i < nbentries; i++) Fill(h->fBuffer[3*i + 2], h->fBuffer[3*i + 3], h->fBuffer[3*i + 1]); // Entries from buffers have to be filled one by one // because FillN doesn't resize histograms. } } if (!initialLimitsFound) return (Int_t) GetEntries(); // all histograms have been processed next.Reset(); } //merge bin contents and errors const Int_t kNstat = 7; Stat_t stats[kNstat], totstats[kNstat]; for (Int_t i=0;i<kNstat;i++) {totstats[i] = stats[i] = 0;} GetStats(totstats); Stat_t nentries = GetEntries(); Int_t binx, biny, ix, iy, nx, ny, bin, ibin; Double_t cu; Int_t nbix = fXaxis.GetNbins(); Bool_t canRebin=TestBit(kCanRebin); ResetBit(kCanRebin); // reset, otherwise setting the under/overflow will rebin while (TH1* h=(TH1*)next()) { // process only if the histogram has limits; otherwise it was processed before if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) { // import statistics h->GetStats(stats); for (Int_t i = 0; i < kNstat; i++) totstats[i] += stats[i]; nentries += h->GetEntries(); nx = h->GetXaxis()->GetNbins(); ny = h->GetYaxis()->GetNbins(); for (biny = 0; biny <= ny + 1; biny++) { iy = fYaxis.FindBin(h->GetYaxis()->GetBinCenter(biny)); for (binx = 0; binx <= nx + 1; binx++) { ix = fXaxis.FindBin(h->GetXaxis()->GetBinCenter(binx)); bin = binx +(nx+2)*biny; ibin = ix +(nbix+2)*iy; cu = h->GetBinContent(bin); if ((!same) && (binx == 0 || binx == nx + 1 || biny == 0 || biny == ny + 1)) { if (cu != 0) { Error("Merge", "Cannot merge histograms - the histograms have" " different limits and undeflows/overflows are present." " The initial histogram is now broken!"); return -1; } } AddBinContent(ibin,cu); if (fSumw2.fN) { Double_t error1 = h->GetBinError(bin); fSumw2.fArray[ibin] += error1*error1; } } } } } if (canRebin) SetBit(kCanRebin); //copy merged stats PutStats(totstats); SetEntries(nentries); if (hclone) delete hclone; return (Long64_t)nentries; } //______________________________________________________________________________ TH2 *TH2::RebinX(Int_t ngroup, const char *newname) { // Rebin only the X axis // see Rebin2D return Rebin2D(ngroup, 1, newname); } //______________________________________________________________________________ TH2 *TH2::RebinY(Int_t ngroup, const char *newname) { // Rebin only the Y axis // see Rebin2D return Rebin2D(1, ngroup, newname); } //______________________________________________________________________________ TH2 *TH2::Rebin2D(Int_t nxgroup, Int_t nygroup, const char *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. // else the current histogram is modified (default) // The parameter 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. // // examples: if hpxpy is an existing TH2 histogram with 40 x 40 bins // hpxpy->Rebin(); // merges two bins along the xaxis and yaxis in one in hpxpy // // Carefull: previous contents of hpxpy are lost // hpxpy->RebinX(5); //merges five bins along the xaxis in one in hpxpy // TH2 *hnew = hpxpy->RebinY(5,"hnew"); // creates a new histogram hnew // // 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, // 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. Int_t i,j,xbin,ybin; Int_t nxbins = fXaxis.GetNbins(); Int_t nybins = fYaxis.GetNbins(); Axis_t xmin = fXaxis.GetXmin(); Axis_t xmax = fXaxis.GetXmax(); Axis_t ymin = fYaxis.GetXmin(); Axis_t ymax = fYaxis.GetXmax(); if ((nxgroup <= 0) || (nxgroup > nxbins)) { Error("Rebin", "Illegal value of nxgroup=%d",nxgroup); return 0; } if ((nygroup <= 0) || (nygroup > nybins)) { Error("Rebin", "Illegal value of nygroup=%d",nygroup); return 0; } Int_t newxbins = nxbins/nxgroup; Int_t newybins = nybins/nygroup; // Save old bin contents into a new array Double_t entries = fEntries; Double_t *oldBins = new Double_t[nxbins*nybins]; for (xbin = 0; xbin < nxbins; xbin++) { for (ybin = 0; ybin < nybins; ybin++) { oldBins[xbin*nybins+ybin] = GetBinContent(xbin+1, ybin+1); } } Double_t *oldErrors = 0; if (fSumw2.fN != 0) { oldErrors = new Double_t[nxbins*nybins]; for (xbin = 0; xbin < nxbins; xbin++) { for (ybin = 0; ybin < nybins; ybin++) { oldErrors[xbin*nybins+ybin] = GetBinError(xbin+1, ybin+1); } } } // create a clone of the old histogram if newname is specified TH2 *hnew = this; if (newname && strlen(newname)) { hnew = (TH2*)Clone(); hnew->SetName(newname); } // change axis specs and rebuild bin contents array 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 } // save the TAttAxis members (reset by SetBins) for x axis Int_t NXdivisions = fXaxis.GetNdivisions(); Color_t XAxisColor = fXaxis.GetAxisColor(); Color_t XLabelColor = fXaxis.GetLabelColor(); Style_t XLabelFont = fXaxis.GetLabelFont(); Float_t XLabelOffset = fXaxis.GetLabelOffset(); Float_t XLabelSize = fXaxis.GetLabelSize(); Float_t XTickLength = fXaxis.GetTickLength(); Float_t XTitleOffset = fXaxis.GetTitleOffset(); Float_t XTitleSize = fXaxis.GetTitleSize(); Color_t XTitleColor = fXaxis.GetTitleColor(); Style_t XTitleFont = fXaxis.GetTitleFont(); // save the TAttAxis members (reset by SetBins) for y axis Int_t NYdivisions = fYaxis.GetNdivisions(); Color_t YAxisColor = fYaxis.GetAxisColor(); Color_t YLabelColor = fYaxis.GetLabelColor(); Style_t YLabelFont = fYaxis.GetLabelFont(); Float_t YLabelOffset = fYaxis.GetLabelOffset(); Float_t YLabelSize = fYaxis.GetLabelSize(); Float_t YTickLength = fYaxis.GetTickLength(); Float_t YTitleOffset = fYaxis.GetTitleOffset(); Float_t YTitleSize = fYaxis.GetTitleSize(); Color_t YTitleColor = fYaxis.GetTitleColor(); Style_t YTitleFont = fYaxis.GetTitleFont(); // 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 Axis_t *xbins = new Axis_t[newxbins+1]; for(i = 0; i <= newxbins; ++i) xbins[i] = fXaxis.GetBinLowEdge(1+i*nxgroup); Axis_t *ybins = new Axis_t[newybins+1]; for(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 } Double_t binContent, binError; Int_t oldxbin = 0; for (xbin = 0; xbin <= newxbins; xbin++) { Int_t oldybin = 0; for (ybin = 0; ybin <= newybins; ybin++) { binContent = 0; binError = 0; for (i = 0; i < nxgroup; i++) { if (oldxbin+i >= nxbins) break; for (j =0; j < nygroup; j++) { if (oldybin+j >= nybins) break; binContent += oldBins[oldybin+j + (oldxbin+i)*nybins]; if (oldErrors) binError += oldErrors[oldybin+ j + (oldxbin+i)*nybins]*oldErrors[oldybin + j + (oldxbin+i)*nybins]; } } hnew->SetBinContent(xbin+1,ybin+1, binContent); if (oldErrors) hnew->SetBinError(xbin+1,ybin+1,TMath::Sqrt(binError)); oldybin += nygroup; } oldxbin += nxgroup; } } // Restore x axis attributes fXaxis.SetNdivisions(NXdivisions); fXaxis.SetAxisColor(XAxisColor); fXaxis.SetLabelColor(XLabelColor); fXaxis.SetLabelFont(XLabelFont); fXaxis.SetLabelOffset(XLabelOffset); fXaxis.SetLabelSize(XLabelSize); fXaxis.SetTickLength(XTickLength); fXaxis.SetTitleOffset(XTitleOffset); fXaxis.SetTitleSize(XTitleSize); fXaxis.SetTitleColor(XTitleColor); fXaxis.SetTitleFont(XTitleFont); // Restore y axis attributes fYaxis.SetNdivisions(NYdivisions); fYaxis.SetAxisColor(YAxisColor); fYaxis.SetLabelColor(YLabelColor); fYaxis.SetLabelFont(YLabelFont); fYaxis.SetLabelOffset(YLabelOffset); fYaxis.SetLabelSize(YLabelSize); fYaxis.SetTickLength(YTickLength); fYaxis.SetTitleOffset(YTitleOffset); fYaxis.SetTitleSize(YTitleSize); fYaxis.SetTitleColor(YTitleColor); fYaxis.SetTitleFont(YTitleFont); hnew->SetEntries(entries); //was modified by SetBinContent delete [] oldBins; if (oldErrors) delete [] oldErrors; return hnew; } //______________________________________________________________________________ TProfile *TH2::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option) const { //*-*-*-*-*Project a 2-D histogram into a profile histogram along X*-*-*-*-*-* //*-* ======================================================== // // The projection is made from the channels along the Y axis // ranging from firstybin to lastybin included. // By default, bins 1 to ny are included // When all bins are included, the number of entries in the projection // is set to the number of entries of the 2-D histogram, otherwise // the number of entries is incremented by 1 for all non empty cells. // // if option "d" is specified, the profile is drawn in the current pad. // // Using a TCutG object, it is possible to select a sub-range of a 2-D histogram. // One must create a graphical cut (mouse or C++) and specify the name // of the cut between [] in the option. // For example, with a TCutG named "cutg", one can call: // myhist->ProfileX(" ",firstybin,lastybin,"[cutg]"); // To invert the cut, it is enough to put a "-" in front of its name: // myhist->ProfileX(" ",firstybin,lastybin,"[-cutg]"); // It is possible to apply several cuts ("," means logical AND): // myhist->ProfileX(" ",firstybin,lastybin,[cutg1,cutg2]"); // // NOTE that if a TProfile named name exists in the current directory or pad, // the histogram is reset and filled again with the current contents of the TH2. // The X axis attributes of the TH2 are copied to the X axis of the profile. TString opt = option; Int_t nx = fXaxis.GetNbins(); Int_t ny = fYaxis.GetNbins(); if (firstybin < 0) firstybin = 1; if (lastybin < 0) lastybin = ny; if (lastybin > ny+1) lastybin = ny; // Create the profile histogram char *pname = (char*)name; if (name && strcmp(name,"_pfx") == 0) { Int_t nch = strlen(GetName()) + 5; pname = new char[nch]; sprintf(pname,"%s%s",GetName(),name); } TProfile *h1=0; //check if a profile with identical name exist TObject *h1obj = gROOT->FindObject(pname); if (h1obj && h1obj->InheritsFrom("TProfile")) { h1 = (TProfile*)h1obj; h1->Reset(); } Int_t ncuts = 0; if (opt.Contains("[")) { ((TH2 *)this)->GetPainter(); if (fPainter) ncuts = fPainter->MakeCuts((char*)opt.Data()); } opt.ToLower(); //must be called after MakeCuts if (!h1) { const TArrayD *bins = fXaxis.GetXbins(); if (bins->fN == 0) { h1 = new TProfile(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),option); } else { h1 = new TProfile(pname,GetTitle(),nx,bins->fArray,option); } } if (pname != name) delete [] pname; // Copy attributes h1->GetXaxis()->ImportAttributes(this->GetXaxis()); h1->SetLineColor(this->GetLineColor()); h1->SetFillColor(this->GetFillColor()); h1->SetMarkerColor(this->GetMarkerColor()); h1->SetMarkerStyle(this->GetMarkerStyle()); // Fill the profile histogram Double_t cont; for (Int_t binx =0;binx<=nx+1;binx++) { for (Int_t biny=firstybin;biny<=lastybin;biny++) { if (ncuts) { if (!fPainter->IsInside(binx,biny)) continue; } cont = GetCellContent(binx,biny); if (cont) { h1->Fill(fXaxis.GetBinCenter(binx),fYaxis.GetBinCenter(biny), cont); } } } if ((firstybin <=1 && lastybin >= ny) && !ncuts) h1->SetEntries(fEntries); if (opt.Contains("d")) { TVirtualPad *padsav = gPad; TVirtualPad *pad = gROOT->GetSelectedPad(); if (pad) pad->cd(); char optin[100]; strcpy(optin,opt.Data()); char *d = (char*)strstr(optin,"d"); if (d) {*d = ' '; if (*(d+1) == 0) *d=0;} char *e = (char*)strstr(optin,"e"); if (e) {*e = ' '; if (*(e+1) == 0) *e=0;} if (!gPad->FindObject(h1)) { h1->Draw(optin); } else { h1->Paint(optin); } if (padsav) padsav->cd(); } return h1; } //______________________________________________________________________________ TProfile *TH2::ProfileY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option) const { //*-*-*-*-*Project a 2-D histogram into a profile histogram along Y*-*-*-*-*-* //*-* ======================================================== // // The projection is made from the channels along the X axis // ranging from firstxbin to lastxbin included. // By default, bins 1 to nx are included // When all bins are included, the number of entries in the projection // is set to the number of entries of the 2-D histogram, otherwise // the number of entries is incremented by 1 for all non empty cells. // // if option "d" is specified, the profile is drawn in the current pad. // // Using a TCutG object, it is possible to select a sub-range of a 2-D histogram. // One must create a graphical cut (mouse or C++) and specify the name // of the cut between [] in the option. // For example, with a TCutG named "cutg", one can call: // myhist->ProfileY(" ",firstybin,lastybin,"[cutg]"); // To invert the cut, it is enough to put a "-" in front of its name: // myhist->ProfileY(" ",firstybin,lastybin,"[-cutg]"); // It is possible to apply several cuts: // myhist->ProfileY(" ",firstybin,lastybin,[cutg1,cutg2]"); // // NOTE that if a TProfile named name exists in the current directory or pad, // the histogram is reset and filled again with the current contents of the TH2. // The Y axis attributes of the TH2 are copied to the X axis of the profile. TString opt = option; Int_t nx = fXaxis.GetNbins(); Int_t ny = fYaxis.GetNbins(); if (firstxbin < 0) firstxbin = 1; if (lastxbin < 0) lastxbin = nx; if (lastxbin > nx+1) lastxbin = nx; // Create the projection histogram char *pname = (char*)name; if (name && strcmp(name,"_pfy") == 0) { Int_t nch = strlen(GetName()) + 5; pname = new char[nch]; sprintf(pname,"%s%s",GetName(),name); } TProfile *h1=0; //check if a profile with identical name exist TObject *h1obj = gROOT->FindObject(pname); if (h1obj && h1obj->InheritsFrom("TProfile")) { h1 = (TProfile*)h1obj; h1->Reset(); } Int_t ncuts = 0; if (opt.Contains("[")) { ((TH2 *)this)->GetPainter(); if (fPainter) ncuts = fPainter->MakeCuts((char*)opt.Data()); } opt.ToLower(); //must be called after MakeCuts if (!h1) { const TArrayD *bins = fYaxis.GetXbins(); if (bins->fN == 0) { h1 = new TProfile(pname,GetTitle(),ny,fYaxis.GetXmin(),fYaxis.GetXmax(),option); } else { h1 = new TProfile(pname,GetTitle(),ny,bins->fArray,option); } } if (pname != name) delete [] pname; // Copy attributes h1->GetXaxis()->ImportAttributes(this->GetYaxis()); h1->SetLineColor(this->GetLineColor()); h1->SetFillColor(this->GetFillColor()); h1->SetMarkerColor(this->GetMarkerColor()); h1->SetMarkerStyle(this->GetMarkerStyle()); // Fill the profile histogram Double_t cont; for (Int_t biny =0;biny<=ny+1;biny++) { for (Int_t binx=firstxbin;binx<=lastxbin;binx++) { if (ncuts) { if (!fPainter->IsInside(binx,biny)) continue; } cont = GetCellContent(binx,biny); if (cont) { h1->Fill(fYaxis.GetBinCenter(biny),fXaxis.GetBinCenter(binx), cont); } } } if ((firstxbin <=1 && lastxbin >= nx) && !ncuts) h1->SetEntries(fEntries); if (opt.Contains("d")) { TVirtualPad *padsav = gPad; TVirtualPad *pad = gROOT->GetSelectedPad(); if (pad) pad->cd(); char optin[100]; strcpy(optin,opt.Data()); char *d = (char*)strstr(optin,"d"); if (d) {*d = ' '; if (*(d+1) == 0) *d=0;} char *e = (char*)strstr(optin,"e"); if (e) {*e = ' '; if (*(e+1) == 0) *e=0;} if (!gPad->FindObject(h1)) { h1->Draw(optin); } else { h1->Paint(optin); } if (padsav) padsav->cd(); } return h1; } //______________________________________________________________________________ TH1D *TH2::ProjectionX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option) const { //*-*-*-*-*Project a 2-D histogram into a 1-D histogram along X*-*-*-*-*-*-* //*-* ==================================================== // // The projection is always of the type TH1D. // The projection is made from the channels along the Y axis // ranging from firstybin to lastybin included. // By default, bins 1 to ny are included // When all bins are included, the number of entries in the projection // is set to the number of entries of the 2-D histogram, otherwise // the number of entries is incremented by 1 for all non empty cells. // // To make the projection in X of the underflow bin in Y, use firstybin=lastybin=0; // To make the projection in X of the overflow bin in Y, use firstybin=lastybin=ny+1; // // if option "e" is specified, the errors are computed. // if option "d" is specified, the projection is drawn in the current pad. // // Using a TCutG object, it is possible to select a sub-range of a 2-D histogram. // One must create a graphical cut (mouse or C++) and specify the name // of the cut between [] in the option. // For example, with a TCutG named "cutg", one can call: // myhist->ProjectionX(" ",firstybin,lastybin,"[cutg]"); // To invert the cut, it is enough to put a "-" in front of its name: // myhist->ProjectionX(" ",firstybin,lastybin,"[-cutg]"); // It is possible to apply several cuts: // myhist->ProjectionX(" ",firstybin,lastybin,[cutg1,cutg2]"); // // NOTE that if a TH1D named name exists in the current directory or pad, // the histogram is reset and filled again with the current contents of the TH2. // The X axis attributes of the TH2 are copied to the X axis of the projection. TString opt = option; Int_t nx = fXaxis.GetNbins(); Int_t ny = fYaxis.GetNbins(); if (firstybin < 0) firstybin = 1; if (lastybin < 0) lastybin = ny; if (lastybin > ny+1) lastybin = ny; // Create the projection histogram char *pname = (char*)name; if (name && strcmp(name,"_px") == 0) { Int_t nch = strlen(GetName()) + 4; pname = new char[nch]; sprintf(pname,"%s%s",GetName(),name); } TH1D *h1=0; //check if histogram with identical name exist TObject *h1obj = gROOT->FindObject(pname); if (h1obj && h1obj->InheritsFrom("TH1D")) { h1 = (TH1D*)h1obj; h1->Reset(); } Int_t ncuts = 0; if (opt.Contains("[")) { ((TH2 *)this)->GetPainter(); if (fPainter) ncuts = fPainter->MakeCuts((char*)opt.Data()); } opt.ToLower(); //must be called after MakeCuts if (!h1) { const TArrayD *bins = fXaxis.GetXbins(); if (bins->fN == 0) { h1 = new TH1D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax()); } else { h1 = new TH1D(pname,GetTitle(),nx,bins->fArray); } if (opt.Contains("e")) h1->Sumw2(); } if (pname != name) delete [] pname; // Copy attributes h1->GetXaxis()->ImportAttributes(this->GetXaxis()); h1->SetLineColor(this->GetLineColor()); h1->SetFillColor(this->GetFillColor()); h1->SetMarkerColor(this->GetMarkerColor()); h1->SetMarkerStyle(this->GetMarkerStyle()); // Fill the projected histogram Double_t cont,err,err2; for (Int_t binx =0;binx<=nx+1;binx++) { err2 = 0; for (Int_t biny=firstybin;biny<=lastybin;biny++) { if (ncuts) { if (!fPainter->IsInside(binx,biny)) continue; } cont = GetCellContent(binx,biny); err = GetCellError(binx,biny); err2 += err*err; if (cont) { h1->Fill(fXaxis.GetBinCenter(binx), cont); } } if (h1->GetSumw2N()) h1->SetBinError(binx,TMath::Sqrt(err2)); } if ((firstybin <=1 && lastybin >= ny) && !ncuts) h1->SetEntries(fEntries); if (opt.Contains("d")) { TVirtualPad *padsav = gPad; TVirtualPad *pad = gROOT->GetSelectedPad(); if (pad) pad->cd(); char optin[100]; strcpy(optin,opt.Data()); char *d = (char*)strstr(optin,"d"); if (d) {*d = ' '; if (*(d+1) == 0) *d=0;} char *e = (char*)strstr(optin,"e"); if (e) {*e = ' '; if (*(e+1) == 0) *e=0;} if (!gPad->FindObject(h1)) { h1->Draw(optin); } else { h1->Paint(optin); } if (padsav) padsav->cd(); } return h1; } //______________________________________________________________________________ TH1D *TH2::ProjectionY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option) const { //*-*-*-*-*Project a 2-D histogram into a 1-D histogram along Y*-*-*-*-*-*-* //*-* ==================================================== // // The projection is always of the type TH1D. // The projection is made from the channels along the X axis // ranging from firstxbin to lastxbin included. // By default, bins 1 to nx are included // When all bins are included, the number of entries in the projection // is set to the number of entries of the 2-D histogram, otherwise // the number of entries is incremented by 1 for all non empty cells. // // To make the projection in Y of the underflow bin in X, use firstxbin=lastxbin=0; // To make the projection in Y of the overflow bin in X, use firstxbin=lastxbin=nx+1; // // if option "e" is specified, the errors are computed. // if option "d" is specified, the projection is drawn in the current pad. // // Using a TCutG object, it is possible to select a sub-range of a 2-D histogram. // One must create a graphical cut (mouse or C++) and specify the name // of the cut between [] in the option. // For example, with a TCutG named "cutg", one can call: // myhist->ProjectionY(" ",firstxbin,lastxbin,"[cutg]"); // To invert the cut, it is enough to put a "-" in front of its name: // myhist->ProjectionY(" ",firstxbin,lastxbin,"[-cutg]"); // It is possible to apply several cuts: // myhist->ProjectionY(" ",firstxbin,lastxbin,[cutg1,cutg2]"); // // NOTE that if a TH1D named name exists in the current directory or pad, // the histogram is reset and filled again with the current contents of the TH2. // The Y axis attributes of the TH2 are copied to the X axis of the projection. TString opt = option; Int_t nx = fXaxis.GetNbins(); Int_t ny = fYaxis.GetNbins(); if (firstxbin < 0) firstxbin = 1; if (lastxbin < 0) lastxbin = nx; if (lastxbin > nx+1) lastxbin = nx; // Create the projection histogram char *pname = (char*)name; if (name && strcmp(name,"_py") == 0) { Int_t nch = strlen(GetName()) + 4; pname = new char[nch]; sprintf(pname,"%s%s",GetName(),name); } TH1D *h1=0; //check if histogram with identical name exist TObject *h1obj = gROOT->FindObject(pname); if (h1obj && h1obj->InheritsFrom("TH1D")) { h1 = (TH1D*)h1obj; h1->Reset(); } Int_t ncuts = 0; if (opt.Contains("[")) { ((TH2 *)this)->GetPainter(); if (fPainter) ncuts = fPainter->MakeCuts((char*)opt.Data()); } opt.ToLower(); //must be called after MakeCuts if (!h1) { const TArrayD *bins = fYaxis.GetXbins(); if (bins->fN == 0) { h1 = new TH1D(pname,GetTitle(),ny,fYaxis.GetXmin(),fYaxis.GetXmax()); } else { h1 = new TH1D(pname,GetTitle(),ny,bins->fArray); } if (opt.Contains("e")) h1->Sumw2(); } if (pname != name) delete [] pname; // Copy attributes h1->GetXaxis()->ImportAttributes(this->GetYaxis()); h1->SetLineColor(this->GetLineColor()); h1->SetFillColor(this->GetFillColor()); h1->SetMarkerColor(this->GetMarkerColor()); h1->SetMarkerStyle(this->GetMarkerStyle()); // Fill the projected histogram Double_t cont,err,err2; for (Int_t biny =0;biny<=ny+1;biny++) { err2 = 0; for (Int_t binx=firstxbin;binx<=lastxbin;binx++) { if (ncuts) { if (!fPainter->IsInside(binx,biny)) continue; } cont = GetCellContent(binx,biny); err = GetCellError(binx,biny); err2 += err*err; if (cont) { h1->Fill(fYaxis.GetBinCenter(biny), cont); } } if (h1->GetSumw2N()) h1->SetBinError(biny,TMath::Sqrt(err2)); } if ((firstxbin <=1 && lastxbin >= nx) && !ncuts) h1->SetEntries(fEntries); if (opt.Contains("d")) { TVirtualPad *padsav = gPad; TVirtualPad *pad = gROOT->GetSelectedPad(); if (pad) pad->cd(); char optin[100]; strcpy(optin,opt.Data()); char *d = (char*)strstr(optin,"d"); if (d) {*d = ' '; if (*(d+1) == 0) *d=0;} char *e = (char*)strstr(optin,"e"); if (e) {*e = ' '; if (*(e+1) == 0) *e=0;} if (!gPad->FindObject(h1)) { h1->Draw(optin); } else { h1->Paint(optin); } if (padsav) padsav->cd(); } return h1; } //______________________________________________________________________________ void TH2::PutStats(Stat_t *stats) { // Replace current statistics with the values in array stats TH1::PutStats(stats); fTsumwy = stats[4]; fTsumwy2 = stats[5]; fTsumwxy = stats[6]; } //______________________________________________________________________________ void TH2::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH1::Reset(option); TString opt = option; opt.ToUpper(); if (opt.Contains("ICE")) return; fTsumwy = 0; fTsumwy2 = 0; fTsumwxy = 0; } //______________________________________________________________________________ void TH2::Streamer(TBuffer &R__b) { // Stream an object of class TH2. if (R__b.IsReading()) { UInt_t R__s, R__c; Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v > 2) { TH2::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c); return; } //====process old versions before automatic schema evolution TH1::Streamer(R__b); R__b >> fScalefactor; R__b >> fTsumwy; R__b >> fTsumwy2; R__b >> fTsumwxy; //====end of old versions } else { TH2::Class()->WriteBuffer(R__b,this); } } ClassImp(TH2C) //______________________________________________________________________________ // TH2C methods //______________________________________________________________________________ TH2C::TH2C(): TH2() { } //______________________________________________________________________________ TH2C::~TH2C() { } //______________________________________________________________________________ TH2C::TH2C(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup) { TArrayC::Set(fNcells); //if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize); } //______________________________________________________________________________ TH2C::TH2C(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xbins,nbinsy,ylow,yup) { TArrayC::Set(fNcells); } //______________________________________________________________________________ TH2C::TH2C(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ybins) { TArrayC::Set(fNcells); } //______________________________________________________________________________ TH2C::TH2C(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayC::Set(fNcells); } //______________________________________________________________________________ TH2C::TH2C(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins ,Int_t nbinsy,const Float_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayC::Set(fNcells); } //______________________________________________________________________________ TH2C::TH2C(const TH2C &h2c) : TH2(), TArrayC() { ((TH2C&)h2c).Copy(*this); } //______________________________________________________________________________ void TH2C::AddBinContent(Int_t bin) { //*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== if (fArray[bin] < 127) fArray[bin]++; } //______________________________________________________________________________ void TH2C::AddBinContent(Int_t bin, Stat_t w) { //*-*-*-*-*-*-*-*-*-*Increment bin content by w*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== Int_t newval = fArray[bin] + Int_t(w); if (newval > -128 && newval < 128) {fArray[bin] = Char_t(newval); return;} if (newval < -127) fArray[bin] = -127; if (newval > 127) fArray[bin] = 127; } //______________________________________________________________________________ void TH2C::Copy(TObject &newth2) const { TH2::Copy((TH2C&)newth2); TArrayC::Copy((TH2C&)newth2); } //______________________________________________________________________________ TH1 *TH2C::DrawCopy(Option_t *option) const { TString opt = option; opt.ToLower(); if (gPad && !opt.Contains("same")) gPad->Clear(); TH2C *newth2 = (TH2C*)Clone(); newth2->SetDirectory(0); newth2->SetBit(kCanDelete); newth2->AppendPad(option); return newth2; } //______________________________________________________________________________ Stat_t TH2C::GetBinContent(Int_t bin) const { if (fBuffer) ((TH2C*)this)->BufferEmpty(); if (bin < 0) bin = 0; if (bin >= fNcells) bin = fNcells-1; if (!fArray) return 0; return Stat_t (fArray[bin]); } //______________________________________________________________________________ void TH2C::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH2::Reset(option); TArrayC::Reset(); } //______________________________________________________________________________ void TH2C::SetBinContent(Int_t bin, Stat_t content) { // Set bin content if (bin < 0) return; if (bin >= fNcells) return; fArray[bin] = Char_t (content); fEntries++; } //______________________________________________________________________________ void TH2C::SetBinsLength(Int_t n) { // Set total number of bins including under/overflow // Reallocate bin contents array if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2); fNcells = n; TArrayC::Set(n); } //______________________________________________________________________________ void TH2C::Streamer(TBuffer &R__b) { // Stream an object of class TH2C. if (R__b.IsReading()) { UInt_t R__s, R__c; Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v > 2) { TH2C::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c); return; } //====process old versions before automatic schema evolution if (R__v < 2) { R__b.ReadVersion(); TH1::Streamer(R__b); TArrayC::Streamer(R__b); R__b.ReadVersion(); R__b >> fScalefactor; R__b >> fTsumwy; R__b >> fTsumwy2; R__b >> fTsumwxy; } else { TH2::Streamer(R__b); TArrayC::Streamer(R__b); R__b.CheckByteCount(R__s, R__c, TH2C::IsA()); } //====end of old versions } else { TH2C::Class()->WriteBuffer(R__b,this); } } //______________________________________________________________________________ TH2C& TH2C::operator=(const TH2C &h1) { if (this != &h1) ((TH2C&)h1).Copy(*this); return *this; } //______________________________________________________________________________ TH2C operator*(Float_t c1, TH2C &h1) { TH2C hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2C operator+(TH2C &h1, TH2C &h2) { TH2C hnew = h1; hnew.Add(&h2,1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2C operator-(TH2C &h1, TH2C &h2) { TH2C hnew = h1; hnew.Add(&h2,-1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2C operator*(TH2C &h1, TH2C &h2) { TH2C hnew = h1; hnew.Multiply(&h2); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2C operator/(TH2C &h1, TH2C &h2) { TH2C hnew = h1; hnew.Divide(&h2); hnew.SetDirectory(0); return hnew; } ClassImp(TH2S) //______________________________________________________________________________ // TH2S methods //______________________________________________________________________________ TH2S::TH2S(): TH2() { } //______________________________________________________________________________ TH2S::~TH2S() { } //______________________________________________________________________________ TH2S::TH2S(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup) { TArrayS::Set(fNcells); //if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize); } //______________________________________________________________________________ TH2S::TH2S(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xbins,nbinsy,ylow,yup) { TArrayS::Set(fNcells); } //______________________________________________________________________________ TH2S::TH2S(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ybins) { TArrayS::Set(fNcells); } //______________________________________________________________________________ TH2S::TH2S(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayS::Set(fNcells); } //______________________________________________________________________________ TH2S::TH2S(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins ,Int_t nbinsy,const Float_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayS::Set(fNcells); } //______________________________________________________________________________ TH2S::TH2S(const TH2S &h2s) : TH2(), TArrayS() { ((TH2S&)h2s).Copy(*this); } //______________________________________________________________________________ void TH2S::AddBinContent(Int_t bin) { //*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== if (fArray[bin] < 32767) fArray[bin]++; } //______________________________________________________________________________ void TH2S::AddBinContent(Int_t bin, Stat_t w) { //*-*-*-*-*-*-*-*-*-*Increment bin content by w*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== Int_t newval = fArray[bin] + Int_t(w); if (newval > -32768 && newval < 32768) {fArray[bin] = Short_t(newval); return;} if (newval < -32767) fArray[bin] = -32767; if (newval > 32767) fArray[bin] = 32767; } //______________________________________________________________________________ void TH2S::Copy(TObject &newth2) const { TH2::Copy((TH2S&)newth2); TArrayS::Copy((TH2S&)newth2); } //______________________________________________________________________________ TH1 *TH2S::DrawCopy(Option_t *option) const { TString opt = option; opt.ToLower(); if (gPad && !opt.Contains("same")) gPad->Clear(); TH2S *newth2 = (TH2S*)Clone(); newth2->SetDirectory(0); newth2->SetBit(kCanDelete); newth2->AppendPad(option); return newth2; } //______________________________________________________________________________ Stat_t TH2S::GetBinContent(Int_t bin) const { if (fBuffer) ((TH2C*)this)->BufferEmpty(); if (bin < 0) bin = 0; if (bin >= fNcells) bin = fNcells-1; if (!fArray) return 0; return Stat_t (fArray[bin]); } //______________________________________________________________________________ void TH2S::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH2::Reset(option); TArrayS::Reset(); } //______________________________________________________________________________ void TH2S::SetBinContent(Int_t bin, Stat_t content) { // Set bin content if (bin < 0) return; if (bin >= fNcells) return; fArray[bin] = Short_t (content); fEntries++; } //______________________________________________________________________________ void TH2S::SetBinsLength(Int_t n) { // Set total number of bins including under/overflow // Reallocate bin contents array if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2); fNcells = n; TArrayS::Set(n); } //______________________________________________________________________________ void TH2S::Streamer(TBuffer &R__b) { // Stream an object of class TH2S. if (R__b.IsReading()) { UInt_t R__s, R__c; Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v > 2) { TH2S::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c); return; } //====process old versions before automatic schema evolution if (R__v < 2) { R__b.ReadVersion(); TH1::Streamer(R__b); TArrayS::Streamer(R__b); R__b.ReadVersion(); R__b >> fScalefactor; R__b >> fTsumwy; R__b >> fTsumwy2; R__b >> fTsumwxy; } else { TH2::Streamer(R__b); TArrayS::Streamer(R__b); R__b.CheckByteCount(R__s, R__c, TH2S::IsA()); } //====end of old versions } else { TH2S::Class()->WriteBuffer(R__b,this); } } //______________________________________________________________________________ TH2S& TH2S::operator=(const TH2S &h1) { if (this != &h1) ((TH2S&)h1).Copy(*this); return *this; } //______________________________________________________________________________ TH2S operator*(Float_t c1, TH2S &h1) { TH2S hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2S operator+(TH2S &h1, TH2S &h2) { TH2S hnew = h1; hnew.Add(&h2,1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2S operator-(TH2S &h1, TH2S &h2) { TH2S hnew = h1; hnew.Add(&h2,-1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2S operator*(TH2S &h1, TH2S &h2) { TH2S hnew = h1; hnew.Multiply(&h2); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2S operator/(TH2S &h1, TH2S &h2) { TH2S hnew = h1; hnew.Divide(&h2); hnew.SetDirectory(0); return hnew; } ClassImp(TH2I) //______________________________________________________________________________ // TH2I methods //______________________________________________________________________________ TH2I::TH2I(): TH2() { } //______________________________________________________________________________ TH2I::~TH2I() { } //______________________________________________________________________________ TH2I::TH2I(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup) { TArrayI::Set(fNcells); //if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize); } //______________________________________________________________________________ TH2I::TH2I(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xbins,nbinsy,ylow,yup) { TArrayI::Set(fNcells); } //______________________________________________________________________________ TH2I::TH2I(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ybins) { TArrayI::Set(fNcells); } //______________________________________________________________________________ TH2I::TH2I(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayI::Set(fNcells); } //______________________________________________________________________________ TH2I::TH2I(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins ,Int_t nbinsy,const Float_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayI::Set(fNcells); } //______________________________________________________________________________ TH2I::TH2I(const TH2I &h2i) : TH2(), TArrayI() { ((TH2I&)h2i).Copy(*this); } //______________________________________________________________________________ void TH2I::AddBinContent(Int_t bin) { //*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== if (fArray[bin] < 2147483647) fArray[bin]++; } //______________________________________________________________________________ void TH2I::AddBinContent(Int_t bin, Stat_t w) { //*-*-*-*-*-*-*-*-*-*Increment bin content by w*-*-*-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== Int_t newval = fArray[bin] + Int_t(w); if (newval > -2147483647 && newval < 2147483647) {fArray[bin] = Int_t(newval); return;} if (newval < -2147483647) fArray[bin] = -2147483647; if (newval > 2147483647) fArray[bin] = 2147483647; } //______________________________________________________________________________ void TH2I::Copy(TObject &newth2) const { TH2::Copy((TH2I&)newth2); TArrayI::Copy((TH2I&)newth2); } //______________________________________________________________________________ TH1 *TH2I::DrawCopy(Option_t *option) const { TString opt = option; opt.ToLower(); if (gPad && !opt.Contains("same")) gPad->Clear(); TH2I *newth2 = (TH2I*)Clone(); newth2->SetDirectory(0); newth2->SetBit(kCanDelete); newth2->AppendPad(option); return newth2; } //______________________________________________________________________________ Stat_t TH2I::GetBinContent(Int_t bin) const { if (fBuffer) ((TH2C*)this)->BufferEmpty(); if (bin < 0) bin = 0; if (bin >= fNcells) bin = fNcells-1; if (!fArray) return 0; return Stat_t (fArray[bin]); } //______________________________________________________________________________ void TH2I::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH2::Reset(option); TArrayI::Reset(); } //______________________________________________________________________________ void TH2I::SetBinContent(Int_t bin, Stat_t content) { // Set bin content if (bin < 0) return; if (bin >= fNcells) return; fArray[bin] = Int_t (content); fEntries++; } //______________________________________________________________________________ void TH2I::SetBinsLength(Int_t n) { // Set total number of bins including under/overflow // Reallocate bin contents array if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2); fNcells = n; TArrayI::Set(n); } //______________________________________________________________________________ TH2I& TH2I::operator=(const TH2I &h1) { if (this != &h1) ((TH2I&)h1).Copy(*this); return *this; } //______________________________________________________________________________ TH2I operator*(Float_t c1, TH2I &h1) { TH2I hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2I operator+(TH2I &h1, TH2I &h2) { TH2I hnew = h1; hnew.Add(&h2,1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2I operator-(TH2I &h1, TH2I &h2) { TH2I hnew = h1; hnew.Add(&h2,-1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2I operator*(TH2I &h1, TH2I &h2) { TH2I hnew = h1; hnew.Multiply(&h2); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2I operator/(TH2I &h1, TH2I &h2) { TH2I hnew = h1; hnew.Divide(&h2); hnew.SetDirectory(0); return hnew; } ClassImp(TH2F) //______________________________________________________________________________ // TH2F methods //______________________________________________________________________________ TH2F::TH2F(): TH2() { } //______________________________________________________________________________ TH2F::~TH2F() { } //______________________________________________________________________________ TH2F::TH2F(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup) { TArrayF::Set(fNcells); //if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize); } //______________________________________________________________________________ TH2F::TH2F(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xbins,nbinsy,ylow,yup) { TArrayF::Set(fNcells); } //______________________________________________________________________________ TH2F::TH2F(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ybins) { TArrayF::Set(fNcells); } //______________________________________________________________________________ TH2F::TH2F(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayF::Set(fNcells); } //______________________________________________________________________________ TH2F::TH2F(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins ,Int_t nbinsy,const Float_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayF::Set(fNcells); } //______________________________________________________________________________ TH2F::TH2F(const TMatrixFBase &m) :TH2("TMatrixFBase","",m.GetNcols(),m.GetColLwb(),1+m.GetColUpb(),m.GetNrows(),m.GetRowLwb(),1+m.GetRowUpb()) { TArrayF::Set(fNcells); Int_t ilow = m.GetRowLwb(); Int_t iup = m.GetRowUpb(); Int_t jlow = m.GetColLwb(); Int_t jup = m.GetColUpb(); for (Int_t i=ilow;i<=iup;i++) { for (Int_t j=jlow;j<=jup;j++) { SetCellContent(j-jlow+1,i-ilow+1,m(i,j)); } } } //______________________________________________________________________________ TH2F::TH2F(const TH2F &h2f) : TH2(), TArrayF() { ((TH2F&)h2f).Copy(*this); } //______________________________________________________________________________ void TH2F::Copy(TObject &newth2) const { TH2::Copy((TH2F&)newth2); TArrayF::Copy((TH2F&)newth2); } //______________________________________________________________________________ TH1 *TH2F::DrawCopy(Option_t *option) const { TString opt = option; opt.ToLower(); if (gPad && !opt.Contains("same")) gPad->Clear(); TH2F *newth2 = (TH2F*)Clone(); newth2->SetDirectory(0); newth2->SetBit(kCanDelete); newth2->AppendPad(option); return newth2; } //______________________________________________________________________________ Stat_t TH2F::GetBinContent(Int_t bin) const { if (fBuffer) ((TH2C*)this)->BufferEmpty(); if (bin < 0) bin = 0; if (bin >= fNcells) bin = fNcells-1; if (!fArray) return 0; return Stat_t (fArray[bin]); } //______________________________________________________________________________ void TH2F::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH2::Reset(option); TArrayF::Reset(); } //______________________________________________________________________________ void TH2F::SetBinContent(Int_t bin, Stat_t content) { // Set bin content if (bin < 0) return; if (bin >= fNcells) return; fArray[bin] = Float_t (content); fEntries++; } //______________________________________________________________________________ void TH2F::SetBinsLength(Int_t n) { // Set total number of bins including under/overflow // Reallocate bin contents array if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2); fNcells = n; TArrayF::Set(n); } //______________________________________________________________________________ void TH2F::Streamer(TBuffer &R__b) { // Stream an object of class TH2F. if (R__b.IsReading()) { UInt_t R__s, R__c; Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v > 2) { TH2F::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c); return; } //====process old versions before automatic schema evolution if (R__v < 2) { R__b.ReadVersion(); TH1::Streamer(R__b); TArrayF::Streamer(R__b); R__b.ReadVersion(); R__b >> fScalefactor; R__b >> fTsumwy; R__b >> fTsumwy2; R__b >> fTsumwxy; } else { TH2::Streamer(R__b); TArrayF::Streamer(R__b); R__b.CheckByteCount(R__s, R__c, TH2F::IsA()); } //====end of old versions } else { TH2F::Class()->WriteBuffer(R__b,this); } } //______________________________________________________________________________ TH2F& TH2F::operator=(const TH2F &h1) { if (this != &h1) ((TH2F&)h1).Copy(*this); return *this; } //______________________________________________________________________________ TH2F operator*(Float_t c1, TH2F &h1) { TH2F hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2F operator*(TH2F &h1, Float_t c1) { TH2F hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2F operator+(TH2F &h1, TH2F &h2) { TH2F hnew = h1; hnew.Add(&h2,1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2F operator-(TH2F &h1, TH2F &h2) { TH2F hnew = h1; hnew.Add(&h2,-1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2F operator*(TH2F &h1, TH2F &h2) { TH2F hnew = h1; hnew.Multiply(&h2); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2F operator/(TH2F &h1, TH2F &h2) { TH2F hnew = h1; hnew.Divide(&h2); hnew.SetDirectory(0); return hnew; } ClassImp(TH2D) //______________________________________________________________________________ // TH2D methods //______________________________________________________________________________ TH2D::TH2D(): TH2() { } //______________________________________________________________________________ TH2D::~TH2D() { } //______________________________________________________________________________ TH2D::TH2D(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup) { TArrayD::Set(fNcells); //if (xlow >= xup || ylow >= yup) SetBuffer(fgBufferSize); } //______________________________________________________________________________ TH2D::TH2D(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,Axis_t ylow,Axis_t yup) :TH2(name,title,nbinsx,xbins,nbinsy,ylow,yup) { TArrayD::Set(fNcells); } //______________________________________________________________________________ TH2D::TH2D(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xlow,xup,nbinsy,ybins) { TArrayD::Set(fNcells); } //______________________________________________________________________________ TH2D::TH2D(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins ,Int_t nbinsy,const Double_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayD::Set(fNcells); } //______________________________________________________________________________ TH2D::TH2D(const char *name,const char *title,Int_t nbinsx,const Float_t *xbins ,Int_t nbinsy,const Float_t *ybins) :TH2(name,title,nbinsx,xbins,nbinsy,ybins) { TArrayD::Set(fNcells); } //______________________________________________________________________________ TH2D::TH2D(const TMatrixDBase &m) :TH2("TMatrixDBase","",m.GetNcols(),m.GetColLwb(),1+m.GetColUpb(),m.GetNrows(),m.GetRowLwb(),1+m.GetRowUpb()) { TArrayD::Set(fNcells); Int_t ilow = m.GetRowLwb(); Int_t iup = m.GetRowUpb(); Int_t jlow = m.GetColLwb(); Int_t jup = m.GetColUpb(); for (Int_t i=ilow;i<=iup;i++) { for (Int_t j=jlow;j<=jup;j++) { SetCellContent(j-jlow+1,i-ilow+1,m(i,j)); } } } //______________________________________________________________________________ TH2D::TH2D(const TH2D &h2d) : TH2(), TArrayD() { ((TH2D&)h2d).Copy(*this); } //______________________________________________________________________________ void TH2D::Copy(TObject &newth2) const { TH2::Copy((TH2D&)newth2); TArrayD::Copy((TH2D&)newth2); } //______________________________________________________________________________ TH1 *TH2D::DrawCopy(Option_t *option) const { TString opt = option; opt.ToLower(); if (gPad && !opt.Contains("same")) gPad->Clear(); TH2D *newth2 = (TH2D*)Clone(); newth2->SetDirectory(0); newth2->SetBit(kCanDelete); newth2->AppendPad(option); return newth2; } //______________________________________________________________________________ Stat_t TH2D::GetBinContent(Int_t bin) const { if (fBuffer) ((TH2C*)this)->BufferEmpty(); if (bin < 0) bin = 0; if (bin >= fNcells) bin = fNcells-1; if (!fArray) return 0; return Stat_t (fArray[bin]); } //______________________________________________________________________________ void TH2D::Reset(Option_t *option) { //*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-* //*-* =========================================== TH2::Reset(option); TArrayD::Reset(); } //______________________________________________________________________________ void TH2D::SetBinContent(Int_t bin, Stat_t content) { // Set bin content if (bin < 0) return; if (bin >= fNcells) return; fArray[bin] = Double_t (content); fEntries++; } //______________________________________________________________________________ void TH2D::SetBinsLength(Int_t n) { // Set total number of bins including under/overflow // Reallocate bin contents array if (n < 0) n = (fXaxis.GetNbins()+2)*(fYaxis.GetNbins()+2); fNcells = n; TArrayD::Set(n); } //______________________________________________________________________________ void TH2D::Streamer(TBuffer &R__b) { // Stream an object of class TH2D. if (R__b.IsReading()) { UInt_t R__s, R__c; Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v > 2) { TH2D::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c); return; } //====process old versions before automatic schema evolution if (R__v < 2) { R__b.ReadVersion(); TH1::Streamer(R__b); TArrayD::Streamer(R__b); R__b.ReadVersion(); R__b >> fScalefactor; R__b >> fTsumwy; R__b >> fTsumwy2; R__b >> fTsumwxy; } else { TH2::Streamer(R__b); TArrayD::Streamer(R__b); R__b.CheckByteCount(R__s, R__c, TH2D::IsA()); } //====end of old versions } else { TH2D::Class()->WriteBuffer(R__b,this); } } //______________________________________________________________________________ TH2D& TH2D::operator=(const TH2D &h1) { if (this != &h1) ((TH2D&)h1).Copy(*this); return *this; } //______________________________________________________________________________ TH2D operator*(Float_t c1, TH2D &h1) { TH2D hnew = h1; hnew.Scale(c1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2D operator+(TH2D &h1, TH2D &h2) { TH2D hnew = h1; hnew.Add(&h2,1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2D operator-(TH2D &h1, TH2D &h2) { TH2D hnew = h1; hnew.Add(&h2,-1); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2D operator*(TH2D &h1, TH2D &h2) { TH2D hnew = h1; hnew.Multiply(&h2); hnew.SetDirectory(0); return hnew; } //______________________________________________________________________________ TH2D operator/(TH2D &h1, TH2D &h2) { TH2D hnew = h1; hnew.Divide(&h2); hnew.SetDirectory(0); return hnew; }