// @(#)root/hist:$Name:  $:$Id: TProfile3D.cxx,v 1.5 2007/02/01 14:58:44 brun Exp $
// Author: Rene Brun   17/05/2006

/*************************************************************************
 * 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 "TProfile3D.h"
#include "THashList.h"
#include "TMath.h"
#include "THLimitsFinder.h"
#include "Riostream.h"
#include "TVirtualPad.h"
#include "TError.h"
#include "TClass.h"

const Int_t kNstat = 13;
Bool_t TProfile3D::fgApproximate = kFALSE;

ClassImp(TProfile3D)

//______________________________________________________________________________
//
//  Profile3D histograms are used to display the mean
//  value of T and its RMS for each cell in X,Y,Z.
//  Profile3D histograms are in many cases an
//  The inter-relation of three measured quantities X, Y, Z and T can always 
//  be visualized by a four-dimensional histogram or scatter-plot; 
//  its representation on the line-printer is not particularly
//  satisfactory, except for sparse data. If T is an unknown (but single-valued)
//  approximate function of X,Y,Z this function is displayed by a profile3D histogram with
//  much better precision than by a scatter-plot.
//
//  The following formulae show the cumulated contents (capital letters) and the values
//  displayed by the printing or plotting routines (small letters) of the elements for cell I, J.
//
//                                                                2
//      H(I,J,K)  =  sum T                      E(I,J,K)  =  sum T
//      l(I,J,K)  =  sum l                      L(I,J,K)  =  sum l
//      h(I,J,K)  =  H(I,J,K)/L(I,J,K)          s(I,J,K)  =  sqrt(E(I,J,K)/L(I,J,K)- h(I,J,K)**2)
//      e(I,J,K)  =  s(I,J,K)/sqrt(L(I,J,K))
//
//  In the special case where s(I,J,K) is zero (eg, case of 1 entry only in one cell)
//  e(I,J,K) is computed from the average of the s(I,J,K) for all cells.
//  This simple/crude approximation was suggested in order to keep the cell
//  during a fit operation.
//
//           Example of a profile3D histogram
//{
//  TCanvas *c1 = new TCanvas("c1","Profile histogram example",200,10,700,500);
//  hprof3d  = new TProfile3D("hprof3d","Profile of pt versus px, py and pz",40,-4,4,40,-4,4,40,0,20);
//  Double_t px, py, pz, pt;
//  TRandom3 r(0);
//  for ( Int_t i=0; i<25000; i++) {
//     r.Rannor(px,py);
//     pz = px*px + py*py;
//     pt = r.Landau(0,1);
//     hprof3d->Fill(px,py,pz,pt,1);
//  }
//  hprof3d->Draw();
//}
//
// NOTE: A TProfile3D is drawn as it was a simple TH3

//______________________________________________________________________________
TProfile3D::TProfile3D() : TH3D()
{
//*-*-*-*-*-*Default constructor for Profile3D histograms*-*-*-*-*-*-*-*-*
//*-*        ============================================
   fTsumwt = fTsumwt2 = 0;
   fScaling = kFALSE;
}

//______________________________________________________________________________
TProfile3D::~TProfile3D()
{
//*-*-*-*-*-*Default destructor for Profile3D histograms*-*-*-*-*-*-*-*-*
//*-*        ===========================================

}

//______________________________________________________________________________
TProfile3D::TProfile3D(const char *name,const char *title,Int_t nx,Double_t xlow,Double_t xup,Int_t ny,Double_t ylow,Double_t yup,Int_t nz, Double_t zlow,Double_t zup,Option_t *option)
    : TH3D(name,title,nx,xlow,xup,ny,ylow,yup,nz,zlow,zup)
{
//*-*-*-*-*-*Normal Constructor for Profile histograms*-*-*-*-*-*-*-*-*-*
//*-*        ==========================================
//
//  The first eleven parameters are similar to TH3D::TH3D.
//  All values of t are accepted at filling time.
//  To fill a profile3D histogram, one must use TProfile3D::Fill function.
//
//  Note that when filling the profile histogram the function Fill
//  checks if the variable t is betyween fTmin and fTmax.
//  If a minimum or maximum value is set for the T scale before filling,
//  then all values below tmin or above tmax will be discarded.
//  Setting the minimum or maximum value for the T scale before filling
//  has the same effect as calling the special TProfile3D constructor below
//  where tmin and tmax are specified.
//
//  H(I,J,K) is printed as the cell contents. The errors computed are s(I,J,K) if CHOPT='S'
//  (spread option), or e(I,J,K) if CHOPT=' ' (error on mean).
//
//        See TProfile3D::BuildOptions for explanation of parameters
//
//   see other constructors below with all possible combinations of
//   fix and variable bin size like in TH3D.

   BuildOptions(0,0,option);
   if (xlow >= xup || ylow >= yup || zlow >= zup) SetBuffer(fgBufferSize);
}

//______________________________________________________________________________
TProfile3D::TProfile3D(const char *name,const char *title,Int_t nx,const Double_t *xbins,Int_t ny,const Double_t *ybins,Int_t nz,const Double_t *zbins,Option_t *option)
    : TH3D(name,title,nx,xbins,ny,ybins,nz,zbins)
{
//  Create a 3-D Profile with variable bins in X , Y and Z

   BuildOptions(0,0,option);
}

//______________________________________________________________________________
void TProfile3D::BuildOptions(Double_t tmin, Double_t tmax, Option_t *option)
{
//*-*-*-*-*-*-*Set Profile3D histogram structure and options*-*-*-*-*-*-*-*-*
//*-*          =============================================
//
//    If a cell has N data points all with the same value T (especially
//    possible when dealing with integers), the spread in T for that cell
//    is zero, and the uncertainty assigned is also zero, and the cell is
//    ignored in making subsequent fits. If SQRT(T) was the correct error
//    in the case above, then SQRT(T)/SQRT(N) would be the correct error here.
//    In fact, any cell with non-zero number of entries N but with zero spread
//    should have an uncertainty SQRT(T)/SQRT(N).
//
//    Now, is SQRT(T)/SQRT(N) really the correct uncertainty?
//    that it is only in the case where the T variable is some sort
//    of counting statistics, following a Poisson distribution. This should
//    probably be set as the default case. However, T can be any variable
//    from an original NTUPLE, not necessarily distributed "Poissonly".
//    The computation of errors is based on the parameter option:
//    option:
//     ' '  (Default) Errors are Spread/SQRT(N) for Spread.ne.0. ,
//                      "     "  SQRT(T)/SQRT(N) for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//     's'            Errors are Spread  for Spread.ne.0. ,
//                      "     "  SQRT(T)  for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//     'i'            Errors are Spread/SQRT(N) for Spread.ne.0. ,
//                      "     "  1./SQRT(12.*N) for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//
//    The third case above corresponds to Integer T values for which the
//    uncertainty is +-0.5, with the assumption that the probability that T
//    takes any value between T-0.5 and T+0.5 is uniform (the same argument
//    goes for T uniformly distributed between T and T+1); this would be
//    useful if T is an ADC measurement, for example. Other, fancier options
//    would be possible, at the cost of adding one more parameter to the PROFILE2D
//    For example, if all T variables are distributed according to some
//    known Gaussian of standard deviation Sigma, then:
//     'G'            Errors are Spread/SQRT(N) for Spread.ne.0. ,
//                      "     "  Sigma/SQRT(N) for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//    For example, this would be useful when all T's are experimental quantities
//    measured with the same instrument with precision Sigma.
//
//

   SetErrorOption(option);

   fBinEntries.Set(fNcells);  //*-* create number of entries per cell array

   Sumw2();                   //*-* create sum of squares of weights array

   fTmin = tmin;
   fTmax = tmax;
   fScaling = kFALSE;
   fTsumwt  = fTsumwt2 = 0;
}

//______________________________________________________________________________
TProfile3D::TProfile3D(const TProfile3D &profile) : TH3D()
{
   //copy constructor
   ((TProfile3D&)profile).Copy(*this);
}


//______________________________________________________________________________
void TProfile3D::Add(TF1 *, Double_t , Option_t*)
{
   // Performs the operation: this = this + c1*f1

   Error("Add","Function not implemented for TProfile3D");
   return;
}


//______________________________________________________________________________
void TProfile3D::Add(const TH1 *h1, Double_t c1)
{
   // Performs the operation: this = this + c1*h1

   if (!h1) {
      Error("Add","Attempt to add a non-existing profile");
      return;
   }
   if (!h1->InheritsFrom("TProfile3D")) {
      Error("Add","Attempt to add a non-profile2D object");
      return;
   }
   TProfile3D *p1 = (TProfile3D*)h1;

//*-*- Check profile compatibility
   Int_t nx = GetNbinsX();
   if (nx != p1->GetNbinsX()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }
   Int_t ny = GetNbinsY();
   if (ny != p1->GetNbinsY()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }
   Int_t nz = GetNbinsZ();
   if (nz != p1->GetNbinsZ()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }

//*-*- Add statistics
   Double_t ac1 = TMath::Abs(c1);
   fEntries += ac1*p1->GetEntries();
   fTsumw   += ac1*p1->fTsumw;
   fTsumw2  += c1*c1*p1->fTsumw2;
   fTsumwx  += ac1*p1->fTsumwx;
   fTsumwx2 += ac1*p1->fTsumwx2;
   fTsumwy  += ac1*p1->fTsumwy;
   fTsumwy2 += ac1*p1->fTsumwy2;
   fTsumwxy += ac1*p1->fTsumwxy;
   fTsumwz  += ac1*p1->fTsumwz;
   fTsumwz2 += ac1*p1->fTsumwz2;
   fTsumwxz += ac1*p1->fTsumwxz;
   fTsumwyz += ac1*p1->fTsumwyz;
   fTsumwt  += ac1*p1->fTsumwt;
   fTsumwt2 += ac1*p1->fTsumwt2;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin,binx,biny,binz;
   Double_t *cu1 = p1->GetW();
   Double_t *er1 = p1->GetW2();
   Double_t *en1 = p1->GetB();
   for (binx =0;binx<=nx+1;binx++) {
      for (biny =0;biny<=ny+1;biny++) {
         for (binz =0;binz<=nz+1;binz++) {
            bin   = GetBin(binx,biny,binz);
            fArray[bin]             += c1*cu1[bin];
            //see http://savannah.cern.ch/bugs/?func=detailitem&item_id=14851
            fSumw2.fArray[bin]      += ac1*er1[bin];
            fBinEntries.fArray[bin] += ac1*en1[bin];
         }
      }
   }
}

//______________________________________________________________________________
void TProfile3D::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
{
//*-*-*-*-*Replace contents of this profile3D by the addition of h1 and h2*-*-*
//*-*      ===============================================================
//
//   this = c1*h1 + c2*h2
//

   if (!h1 || !h2) {
      Error("Add","Attempt to add a non-existing profile");
      return;
   }
   if (!h1->InheritsFrom("TProfile3D")) {
      Error("Add","Attempt to add a non-profile3D object");
      return;
   }
   TProfile3D *p1 = (TProfile3D*)h1;
   if (!h2->InheritsFrom("TProfile3D")) {
      Error("Add","Attempt to add a non-profile3D object");
      return;
   }
   TProfile3D *p2 = (TProfile3D*)h2;

//*-*- Check profile compatibility
   Int_t nx = GetNbinsX();
   if (nx != p1->GetNbinsX() || nx != p2->GetNbinsX()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }
   Int_t ny = GetNbinsY();
   if (ny != p1->GetNbinsY() || ny != p2->GetNbinsY()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }
   Int_t nz = GetNbinsZ();
   if (nz != p1->GetNbinsZ() || nz != p2->GetNbinsZ()) {
      Error("Add","Attempt to add profiles with different number of bins");
      return;
   }

//*-*- Add statistics
   Double_t ac1 = TMath::Abs(c1);
   Double_t ac2 = TMath::Abs(c2);
   fEntries = ac1*p1->GetEntries() + ac2*p2->GetEntries();
   fTsumw   = ac1*p1->fTsumw       + ac2*p2->fTsumw;
   fTsumw2  = c1*c1*p1->fTsumw2    + c2*c2*p2->fTsumw2;
   fTsumwx  = ac1*p1->fTsumwx      + ac2*p2->fTsumwx;
   fTsumwx2 = ac1*p1->fTsumwx2     + ac2*p2->fTsumwx2;
   fTsumwy  = ac1*p1->fTsumwy      + ac2*p2->fTsumwy;
   fTsumwy2 = ac1*p1->fTsumwy2     + ac2*p2->fTsumwy2;
   fTsumwxy = ac1*p1->fTsumwxy     + ac2*p2->fTsumwxy;
   fTsumwz  = ac1*p1->fTsumwz      + ac2*p2->fTsumwz;
   fTsumwz2 = ac1*p1->fTsumwz2     + ac2*p2->fTsumwz2;
   fTsumwxz = ac1*p1->fTsumwxz     + ac2*p2->fTsumwxz;
   fTsumwyz = ac1*p1->fTsumwyz     + ac2*p2->fTsumwyz;
   fTsumwt  = ac1*p1->fTsumwt      + ac2*p2->fTsumwt;
   fTsumwt2 = ac1*p1->fTsumwt2     + ac2*p2->fTsumwt2;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin,binx,biny,binz;
   Double_t *cu1 = p1->GetW();
   Double_t *cu2 = p2->GetW();
   Double_t *er1 = p1->GetW2();
   Double_t *er2 = p2->GetW2();
   Double_t *en1 = p1->GetB();
   Double_t *en2 = p2->GetB();
   for (binx =0;binx<=nx+1;binx++) {
      for (biny =0;biny<=ny+1;biny++) {
         for (binz =0;binz<=nz+1;binz++) {
            bin   = GetBin(binx,biny,binz);
            fArray[bin]  = c1*cu1[bin] +  c2*cu2[bin];
            if (fScaling) {
               //see http://savannah.cern.ch/bugs/?func=detailitem&item_id=14851
               fSumw2.fArray[bin]      = ac1*ac1*er1[bin] + ac2*ac2*er2[bin];
               fBinEntries.fArray[bin] = en1[bin];
            } else {
               fSumw2.fArray[bin]      = ac1*er1[bin] + ac2*er2[bin];
               fBinEntries.fArray[bin] = ac1*en1[bin] + ac2*en2[bin];
            }
         }
      }
   }
}


//______________________________________________________________________________
void TProfile3D::Approximate(Bool_t approx)
{
//     static function
// set the fgApproximate flag. When the flag is true, the function GetBinError
// will approximate the bin error with the average profile error on all bins
// in the following situation only
//  - the number of bins in the profile3D is less than 10404 (eg 100x100x100)
//  - the bin number of entries is small ( <5)
//  - the estimated bin error is extremely small compared to the bin content
//  (see TProfile3D::GetBinError)

   fgApproximate = approx;
}


//______________________________________________________________________________
Int_t TProfile3D::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;
      Double_t zmin = fBuffer[4];
      Double_t zmax = zmin;
      for (Int_t i=1;i<nbentries;i++) {
         Double_t x = fBuffer[5*i+2];
         if (x < xmin) xmin = x;
         if (x > xmax) xmax = x;
         Double_t y = fBuffer[5*i+3];
         if (y < ymin) ymin = y;
         if (y > ymax) ymax = y;
         Double_t z = fBuffer[5*i+4];
         if (z < zmin) zmin = z;
         if (z > zmax) zmax = z;
     }
      if (fXaxis.GetXmax() <= fXaxis.GetXmin() || fYaxis.GetXmax() <= fYaxis.GetXmin() || fZaxis.GetXmax() <= fZaxis.GetXmin()) {
         THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax,ymin,ymax,zmin,zmax);
      } 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");
         if (zmin <  fZaxis.GetXmin()) RebinAxis(zmin,"Z");
         if (zmax >= fZaxis.GetXmax()) RebinAxis(zmax,"Z");
         fBuffer = buffer;
         fBufferSize = keep;
      }
   }

   fBuffer = 0;
   for (Int_t i=0;i<nbentries;i++) {
      Fill(buffer[5*i+2],buffer[5*i+3],buffer[5*i+4],buffer[5*i+5],buffer[5*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 TProfile3D::BufferFill(Double_t x, Double_t y, Double_t z, Double_t t, Double_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
// fBuffer[4] = z of first entry
// fBuffer[5] = t 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 (5*nbentries+5 >= fBufferSize) {
      BufferEmpty(1);
      return Fill(x,y,z,t,w);
   }
   fBuffer[5*nbentries+1] = w;
   fBuffer[5*nbentries+2] = x;
   fBuffer[5*nbentries+3] = y;
   fBuffer[5*nbentries+4] = z;
   fBuffer[5*nbentries+5] = t;
   fBuffer[0] += 1;
   return -2;
}

//______________________________________________________________________________
void TProfile3D::Copy(TObject &obj) const
{
//*-*-*-*-*-*-*-*Copy a Profile3D histogram to a new profile2D histogram*-*-*-*
//*-*            =======================================================

   TH3D::Copy(((TProfile3D&)obj));
   fBinEntries.Copy(((TProfile3D&)obj).fBinEntries);
   ((TProfile3D&)obj).fTmin = fTmin;
   ((TProfile3D&)obj).fTmax = fTmax;
   ((TProfile3D&)obj).fScaling = fScaling;
   ((TProfile3D&)obj).fErrorMode = fErrorMode;
   ((TProfile3D&)obj).fTsumwt  = fTsumwt;
   ((TProfile3D&)obj).fTsumwt2 = fTsumwt2;
}


//______________________________________________________________________________
void TProfile3D::Divide(TF1 *, Double_t )
{
   // Performs the operation: this = this/(c1*f1)

   Error("Divide","Function not implemented for TProfile3D");
   return;
}

//______________________________________________________________________________
void TProfile3D::Divide(const TH1 *h1)
{
//*-*-*-*-*-*-*-*-*-*-*Divide this profile2D by h1*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ===========================
//
//   this = this/h1
//

   if (!h1) {
      Error("Divide","Attempt to divide a non-existing profile2D");
      return;
   }
   if (!h1->InheritsFrom("TProfile3D")) {
      Error("Divide","Attempt to divide a non-profile2D object");
      return;
   }
   TProfile3D *p1 = (TProfile3D*)h1;

//*-*- Check profile compatibility
   Int_t nx = GetNbinsX();
   if (nx != p1->GetNbinsX()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }
   Int_t ny = GetNbinsY();
   if (ny != p1->GetNbinsY()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }
   Int_t nz = GetNbinsZ();
   if (nz != p1->GetNbinsZ()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }

//*-*- Reset statistics
   fEntries = fTsumw   = fTsumw2 = fTsumwx = fTsumwx2 = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin,binx,biny,binz;
   Double_t *cu1 = p1->GetW();
   Double_t *er1 = p1->GetW2();
   Double_t *en1 = p1->GetB();
   Double_t c0,c1,w,u,x,y,z;
   for (binx =0;binx<=nx+1;binx++) {
      for (biny =0;biny<=ny+1;biny++) {
         for (binz =0;binz<=nz+1;binz++) {
            bin   = GetBin(binx,biny,binz);
            c0  = fArray[bin];
            c1  = cu1[bin];
            if (c1) w = c0/c1;
            else    w = 0;
            fArray[bin] = w;
            u = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            y = fYaxis.GetBinCenter(biny);
            z = fZaxis.GetBinCenter(binz);
            fEntries++;
            fTsumw   += u;
            fTsumw2  += u*u;
            fTsumwx  += u*x;
            fTsumwx2 += u*x*x;
            fTsumwy  += u*y;
            fTsumwy2 += u*y*y;
            fTsumwxy += u*x*y;
            fTsumwz  += u;
            fTsumwz2 += u*z;
            fTsumwxz += u*x*z;
            fTsumwyz += u*y*z;
            fTsumwt  += u;
            fTsumwt2 += u*u;
            Double_t e0 = fSumw2.fArray[bin];
            Double_t e1 = er1[bin];
            Double_t c12= c1*c1;
            if (!c1) fSumw2.fArray[bin] = 0;
            else     fSumw2.fArray[bin] = (e0*c1*c1 + e1*c0*c0)/(c12*c12);
            if (!en1[bin]) fBinEntries.fArray[bin] = 0;
            else           fBinEntries.fArray[bin] /= en1[bin];
         }
      }
   }
}


//______________________________________________________________________________
void TProfile3D::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
//*-*-*-*-*Replace contents of this profile2D by the division of h1 by h2*-*-*
//*-*      ==============================================================
//
//   this = c1*h1/(c2*h2)
//

   TString opt = option;
   opt.ToLower();
   Bool_t binomial = kFALSE;
   if (opt.Contains("b")) binomial = kTRUE;
   if (!h1 || !h2) {
      Error("Divide","Attempt to divide a non-existing profile2D");
      return;
   }
   if (!h1->InheritsFrom("TProfile3D")) {
      Error("Divide","Attempt to divide a non-profile2D object");
      return;
   }
   TProfile3D *p1 = (TProfile3D*)h1;
   if (!h2->InheritsFrom("TProfile3D")) {
      Error("Divide","Attempt to divide a non-profile2D object");
      return;
   }
   TProfile3D *p2 = (TProfile3D*)h2;

//*-*- Check histogram compatibility
   Int_t nx = GetNbinsX();
   if (nx != p1->GetNbinsX() || nx != p2->GetNbinsX()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }
   Int_t ny = GetNbinsY();
   if (ny != p1->GetNbinsY() || ny != p2->GetNbinsY()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }
   Int_t nz = GetNbinsZ();
   if (nz != p1->GetNbinsZ() || nz != p2->GetNbinsZ()) {
      Error("Divide","Attempt to divide profiles with different number of bins");
      return;
   }
   if (!c2) {
      Error("Divide","Coefficient of dividing profile cannot be zero");
      return;
   }

//*-*- Reset statistics
   fEntries = fTsumw   = fTsumw2 = fTsumwx = fTsumwx2 = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin,binx,biny,binz;
   Double_t *cu1 = p1->GetW();
   Double_t *cu2 = p2->GetW();
   Double_t *er1 = p1->GetW2();
   Double_t *er2 = p2->GetW2();
   Double_t *en1 = p1->GetB();
   Double_t *en2 = p2->GetB();
   Double_t b1,b2,w,u,x,y,z,ac1,ac2;
   ac1 = TMath::Abs(c1);
   ac2 = TMath::Abs(c2);
   for (binx =0;binx<=nx+1;binx++) {
      for (biny =0;biny<=ny+1;biny++) {
         for (binz =0;binz<=nz+1;binz++) {
            bin   = GetBin(binx,biny,binz);
            b1  = cu1[bin];
            b2  = cu2[bin];
            if (b2) w = c1*b1/(c2*b2);
            else    w = 0;
            fArray[bin] = w;
            u = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            y = fYaxis.GetBinCenter(biny);
            z = fZaxis.GetBinCenter(biny);
            fEntries++;
            fTsumw   += u;
            fTsumw2  += u*u;
            fTsumwx  += u*x;
            fTsumwx2 += u*x*x;
            fTsumwy  += u*y;
            fTsumwy2 += u*y*y;
            fTsumwxy += u*x*y;
            fTsumwz  += u*z;
            fTsumwz2 += u*z*z;
            fTsumwxz += u*x*z;
            fTsumwyz += u*y*z;
            fTsumwt  += u;
            fTsumwt2 += u*u;
            Double_t e1 = er1[bin];
            Double_t e2 = er2[bin];
          //Double_t b22= b2*b2*d2;
            Double_t b22= b2*b2*TMath::Abs(c2);
            if (!b2) fSumw2.fArray[bin] = 0;
            else {
               if (binomial) {
                  fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));
               } else {
                  fSumw2.fArray[bin] = ac1*ac2*(e1*b2*b2 + e2*b1*b1)/(b22*b22);
               }
            }
            if (!en2[bin]) fBinEntries.fArray[bin] = 0;
            else           fBinEntries.fArray[bin] = en1[bin]/en2[bin];
         }
      }
   }
}

//______________________________________________________________________________
TH1 *TProfile3D::DrawCopy(Option_t *option) const
{
//*-*-*-*-*-*-*-*Draw a copy of this profile3D histogram*-*-*-*-*-*-*-*-*-*-*
//*-*            =======================================
   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TProfile3D *newpf = new TProfile3D();
   Copy(*newpf);
   newpf->SetDirectory(0);
   newpf->SetBit(kCanDelete);
   newpf->AppendPad(option);
   return newpf;
}

//______________________________________________________________________________
Int_t TProfile3D::Fill(Double_t x, Double_t y, Double_t z, Double_t t)
{
//*-*-*-*-*-*-*-*-*-*-*Fill a Profile3D histogram (no weights)*-*-*-*-*-*-*-*
//*-*                  =======================================

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

   Int_t bin,binx,biny,binz;

   if (fTmin != fTmax) {
      if (t <fTmin || t> fTmax) return -1;
   }

   fEntries++;
   binx =fXaxis.FindBin(x);
   biny =fYaxis.FindBin(y);
   binz =fZaxis.FindBin(z);
   bin  = GetBin(binx,biny,binz);
   AddBinContent(bin, t);
   fSumw2.fArray[bin] += (Double_t)t*t;
   fBinEntries.fArray[bin] += 1;
   if (binx == 0 || binx > fXaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   if (biny == 0 || biny > fYaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   if (binz == 0 || binz > fZaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
//printf("x=%g, y=%g, z=%g, t=%g, binx=%d, biny=%d, binz=%d, bin=%d\n",x,y,z,t,binx,biny,binz,bin);
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   fTsumwy  += y;
   fTsumwy2 += y*y;
   fTsumwxy += x*y;
   fTsumwz  += z;
   fTsumwz2 += z*z;
   fTsumwxz += x*z;
   fTsumwyz += y*z;
   fTsumwt  += t;
   fTsumwt2 += t*t;
   return bin;
}

//______________________________________________________________________________
Int_t TProfile3D::Fill(Double_t x, Double_t y, Double_t z, Double_t t, Double_t w)
{
//*-*-*-*-*-*-*-*-*-*-*Fill a Profile3D histogram with weights*-*-*-*-*-*-*-*
//*-*                  =======================================

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

   Int_t bin,binx,biny,binz;

   if (fTmin != fTmax) {
      if (t <fTmin || z> fTmax) return -1;
   }

   Double_t u= (w > 0 ? w : -w);
   fEntries++;
   binx =fXaxis.FindBin(x);
   biny =fYaxis.FindBin(y);
   binz =fZaxis.FindBin(z);
   bin  = GetBin(binx,biny,binz);
   AddBinContent(bin, u*t);
   fSumw2.fArray[bin] += u*t*t;
   fBinEntries.fArray[bin] += u;
   if (binx == 0 || binx > fXaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   if (biny == 0 || biny > fYaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   if (binz == 0 || binz > fZaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   fTsumw   += u;
   fTsumw2  += u*u;
   fTsumwx  += u*x;
   fTsumwx2 += u*x*x;
   fTsumwy  += u*y;
   fTsumwy2 += u*y*y;
   fTsumwxy += u*x*y;
   fTsumwz  += u*z;
   fTsumwz2 += u*z*z;
   fTsumwxz += u*x*z;
   fTsumwyz += u*y*z;
   fTsumwt  += u*t;
   fTsumwt2 += u*t*t;
   return bin;
}

//______________________________________________________________________________
Double_t TProfile3D::GetBinContent(Int_t bin) const
{
//*-*-*-*-*-*-*Return bin content of a Profile3D histogram*-*-*-*-*-*-*-*-*
//*-*          ===========================================

   if (fBuffer) ((TProfile3D*)this)->BufferEmpty();

   if (bin < 0 || bin >= fNcells) return 0;
   if (fBinEntries.fArray[bin] == 0) return 0;
   if (!fArray) return 0;
   return fArray[bin]/fBinEntries.fArray[bin];
}

//______________________________________________________________________________
Double_t TProfile3D::GetBinEntries(Int_t bin) const
{
//*-*-*-*-*-*-*Return bin entries of a Profile3D histogram*-*-*-*-*-*-*-*-*
//*-*          ===========================================

   if (fBuffer) ((TProfile3D*)this)->BufferEmpty();

   if (bin < 0 || bin >= fNcells) return 0;
   return fBinEntries.fArray[bin];
}

//______________________________________________________________________________
Double_t TProfile3D::GetBinError(Int_t bin) const
{
//*-*-*-*-*-*-*Return bin error of a Profile3D histogram*-*-*-*-*-*-*-*-*
//
// Computing errors: A moving field
// =================================
// The computation of errors for a TProfile3D has evolved with the versions
// of ROOT. The difficulty is in computing errors for bins with low statistics.
// - prior to version 3.10, we had no special treatment of low statistic bins.
//   As a result, these bins had huge errors. The reason is that the
//   expression eprim2 is very close to 0 (rounding problems) or 0.
// - The algorithm is modified/protected for the case
//   when a TProfile3D is projected (ProjectionX). The previous algorithm
//   generated a N^2 problem when projecting a TProfile3D with a large number of
//   bins (eg 100000).
// - in version 3.10/02, a new static function TProfile::Approximate
//   is introduced to enable or disable (default) the approximation.
//   (see also comments in TProfile::GetBinError)

   if (fBuffer) ((TProfile3D*)this)->BufferEmpty();

   if (bin < 0 || bin >= fNcells) return 0;
   Double_t cont = fArray[bin];
   Double_t sum  = fBinEntries.fArray[bin];
   Double_t err2 = fSumw2.fArray[bin];
   if (sum == 0) return 0;
   Double_t eprim;
   Double_t contsum = cont/sum;
   Double_t eprim2  = TMath::Abs(err2/sum - contsum*contsum);
   eprim          = TMath::Sqrt(eprim2);
   Double_t test = 1;
   if (err2 != 0 && sum < 5) test = eprim2*sum/err2;
   //if (eprim <= 0) { //was in 3.10/01
   if (fgApproximate && fNcells <=1000404 && (test < 1.e-4 || eprim2 < 1e-6)) { //in 3.10/02
      Double_t scont, ssum, serr2;
      scont = ssum = serr2 = 0;
      for (Int_t i=1;i<fNcells;i++) {
         if (fSumw2.fArray[i] <= 0) continue; //added in 3.10/02
         scont += fArray[i];
         ssum  += fBinEntries.fArray[i];
         serr2 += fSumw2.fArray[i];
      }
      Double_t scontsum = scont/ssum;
      Double_t seprim2  = TMath::Abs(serr2/ssum - scontsum*scontsum);
      eprim           = TMath::Sqrt(seprim2);
   }
   if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
   else if (fErrorMode == kERRORSPREAD) return eprim;
   else if (fErrorMode == kERRORSPREADI) {
      if (eprim != 0) return eprim/TMath::Sqrt(sum);
      return 1/TMath::Sqrt(12*sum);
   }
   else if (fErrorMode == kERRORSPREADG) {
      return eprim/TMath::Sqrt(sum);
   }
   else return eprim;
}

//______________________________________________________________________________
Option_t *TProfile3D::GetErrorOption() const
{
//*-*-*-*-*-*-*-*-*-*Return option to compute profile2D errors*-*-*-*-*-*-*-*
//*-*                =========================================

   if (fErrorMode == kERRORSPREAD)  return "s";
   if (fErrorMode == kERRORSPREADI) return "i";
   if (fErrorMode == kERRORSPREADG) return "g";
   return "";
}

//______________________________________________________________________________
void TProfile3D::GetStats(Double_t *stats) const
{
   // fill the array stats from the contents of this profile
   // 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
   // stats[7] = sumwz
   // stats[8] = sumwz2
   // stats[9] = sumwxz
   // stats[10]= sumwyz
   // stats[11]= sumwt
   // stats[12]= sumwt2
   //
   // If no axis-subrange is specified (via TAxis::SetRange), the array stats
   // is simply a copy of the statistics quantities computed at filling time.
   // If a sub-range is specified, the function recomputes these quantities
   // from the bin contents in the current axis range.

   if (fBuffer) ((TProfile3D*)this)->BufferEmpty();

   // Loop on bins
   if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange) || fYaxis.TestBit(TAxis::kAxisRange)) {
      Int_t bin, binx, biny,binz;
      Double_t w;
      Double_t x,y,z;
      for (bin=0;bin<kNstat;bin++) stats[bin] = 0;
      if (!fBinEntries.fArray) return;
      for (binz=fZaxis.GetFirst();binz<=fZaxis.GetLast();biny++) {
         z = fZaxis.GetBinCenter(binz);
         for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
            y = fYaxis.GetBinCenter(biny);
            for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
               bin = GetBin(binx,biny,binz);
               w         = fBinEntries.fArray[bin];
               x         = fXaxis.GetBinCenter(binx);
               stats[0]  += w;
               stats[1]  += w*w;
               stats[2]  += w*x;
               stats[3]  += w*x*x;
               stats[4]  += w*y;
               stats[5]  += w*y*y;
               stats[6]  += w*x*y;
               stats[7]  += w*z;
               stats[8]  += w*z*z;
               stats[9]  += w*x*z;
               stats[10] += w*y*z;
               stats[11] += fArray[bin];
               stats[12] += fSumw2.fArray[bin];
            }
         }
      }
   } else {
      stats[0]  = fTsumw;
      stats[1]  = fTsumw2;
      stats[2]  = fTsumwx;
      stats[3]  = fTsumwx2;
      stats[4]  = fTsumwy;
      stats[5]  = fTsumwy2;
      stats[6]  = fTsumwxy;
      stats[7]  = fTsumwz;
      stats[8]  = fTsumwz2;
      stats[9]  = fTsumwxz;
      stats[10] = fTsumwyz;
      stats[11] = fTsumwt;
      stats[12] = fTsumwt2;
   }
}

//______________________________________________________________________________
Long64_t TProfile3D::Merge(TCollection *li)
{
   //Merge all histograms in the collection in 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 (!li) return 0;
   if (li->IsEmpty()) return (Int_t) GetEntries();

   TList inlist;
   TH1* hclone = (TH1*)Clone("FirstClone");
   R__ASSERT(hclone);
   BufferEmpty(1);         // To remove buffer.
   Reset();                // BufferEmpty sets limits so we can't use it later.
   SetEntries(0);
   inlist.Add(hclone);
   inlist.AddAll(li);

   TAxis newXAxis;
   TAxis newYAxis;
   TAxis newZAxis;
   Bool_t initialLimitsFound = kFALSE;
   Bool_t same = kTRUE;
   Bool_t allHaveLimits = kTRUE;

   TIter next(&inlist);
   while (TObject *o = next()) {
      TProfile3D* h = dynamic_cast<TProfile3D*> (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());
            newZAxis.Set(h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
                     h->GetZaxis()->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());
            }
            if (!RecomputeAxisLimits(newZAxis, *(h->GetZaxis()))) {
               Error("Merge", "Cannot merge histograms - limits are inconsistent:\n "
                     "first: (%d, %f, %f), second: (%d, %f, %f)",
                     newZAxis.GetNbins(), newZAxis.GetXmin(), newZAxis.GetXmax(),
                     h->GetZaxis()->GetNbins(), h->GetZaxis()->GetXmin(),
                     h->GetZaxis()->GetXmax());
            }
         }
      }
   }
   next.Reset();

   same = same && SameLimitsAndNBins(newXAxis, *GetXaxis())
               && SameLimitsAndNBins(newYAxis, *GetYaxis())
               && SameLimitsAndNBins(newZAxis, *GetZaxis());
   if (!same && initialLimitsFound)
      SetBins(newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
              newYAxis.GetNbins(), newYAxis.GetXmin(), newYAxis.GetXmax(),
              newZAxis.GetNbins(), newZAxis.GetXmin(), newZAxis.GetXmax());

   if (!allHaveLimits) {
      // fill this histogram with all the data from buffers of histograms without limits
      while (TProfile3D* h = dynamic_cast<TProfile3D*> (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[5*i + 2], h->fBuffer[5*i + 3],
                    h->fBuffer[5*i + 4], h->fBuffer[5*i + 5], h->fBuffer[5*i + 1]);
         }
      }
      if (!initialLimitsFound)
         return (Int_t) GetEntries();  // all histograms have been processed
      next.Reset();
   }

   //merge bin contents and errors
   Double_t stats[kNstat], totstats[kNstat];
   for (Int_t i=0;i<kNstat;i++) {totstats[i] = stats[i] = 0;}
   GetStats(totstats);
   Double_t nentries = GetEntries();
   Int_t binx, biny, binz, ix, iy, iz, nx, ny, nz, bin, ibin;
   Int_t nbix = fXaxis.GetNbins();
   Int_t nbiy = fYaxis.GetNbins();
   Bool_t canRebin=TestBit(kCanRebin);
   ResetBit(kCanRebin); // reset, otherwise setting the under/overflow will rebin

   while (TProfile3D* h=(TProfile3D*)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();
         nz = h->GetZaxis()->GetNbins();

         for (binz = 0; binz <= nz + 1; binz++) {
            iz = fZaxis.FindBin(h->GetZaxis()->GetBinCenter(binz));
            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 + (ny+2)*binz);
                  ibin = ix   +(nbix+2)*(iy + (nbiy+2)*iz);
                  if ((!same) && (binx == 0 || binx == nx + 1
                               || biny == 0 || biny == ny + 1 
                               || binz == 0 || binz == nz + 1)) {
                     if (h->GetW()[bin] != 0) {
                        Error("Merge", "Cannot merge histograms - the histograms have"
                                       " different limits and undeflows/overflows are present."
                                       " The initial histogram is now broken!");
                        return -1;
                     }
                  }
                  fArray[ibin]             += h->GetW()[bin];
                  fSumw2.fArray[ibin]      += h->GetW2()[bin];
                  fBinEntries.fArray[ibin] += h->GetB()[bin];
               }
            }
         }
         fEntries += h->GetEntries();
         fTsumw   += h->fTsumw;
         fTsumw2  += h->fTsumw2;
         fTsumwx  += h->fTsumwx;
         fTsumwx2 += h->fTsumwx2;
         fTsumwy  += h->fTsumwy;
         fTsumwy2 += h->fTsumwy2;
         fTsumwxy += h->fTsumwxy;
         fTsumwz  += h->fTsumwz;
         fTsumwz2 += h->fTsumwz2;
         fTsumwxz += h->fTsumwxz;
         fTsumwyz += h->fTsumwyz;
         fTsumwt  += h->fTsumwt;
         fTsumwt2 += h->fTsumwt2;
      }
   }
   if (canRebin) SetBit(kCanRebin);

   //copy merged stats
   PutStats(totstats);
   SetEntries(nentries);
   inlist.Remove(hclone);
   delete hclone;
   return (Long64_t)nentries;
}

//______________________________________________________________________________
void TProfile3D::Multiply(TF1 *, Double_t )
{
   // Performs the operation: this = this*c1*f1

   Error("Multiply","Function not implemented for TProfile3D");
   return;
}

//______________________________________________________________________________
void TProfile3D::Multiply(const TH1 *)
{
//*-*-*-*-*-*-*-*-*-*-*Multiply this profile2D by h1*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =============================
//
//   this = this*h1
//
   Error("Multiply","Multiplication of profile2D histograms not implemented");
}


//______________________________________________________________________________
void TProfile3D::Multiply(const TH1 *, const TH1 *, Double_t, Double_t, Option_t *)
{
//*-*-*-*-*Replace contents of this profile2D by multiplication of h1 by h2*-*
//*-*      ================================================================
//
//   this = (c1*h1)*(c2*h2)
//

   Error("Multiply","Multiplication of profile2D histograms not implemented");
}

//______________________________________________________________________________
TH3D *TProfile3D::ProjectionXYZ(const char *name, Option_t *option) const
{
//*-*-*-*-*Project this profile3D into a 3-D histogram along X,Y,Z*-*-*-*-*-*-*
//*-*      =====================================================
//
//   The projection is always of the type TH3D.
//
//   if option "E" is specified, the errors are computed. (default)
//   if option "B" is specified, the content of bin of the returned histogram
//      will be equal to the GetBinEntries(bin) of the profile,
//   if option "C=E" the bin contents of the projection are set to the
//       bin errors of the profile
//

   TString opt = option;
   opt.ToLower();
   Int_t nx = fXaxis.GetNbins();
   Int_t ny = fYaxis.GetNbins();
   Int_t nz = fZaxis.GetNbins();

   // Create the projection histogram
   char *pname = (char*)name;
   if (strcmp(name,"_px") == 0) {
      Int_t nch = strlen(GetName()) + 4;
      pname = new char[nch];
      sprintf(pname,"%s%s",GetName(),name);
   }
   TH3D *h1 = new TH3D(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),ny,fYaxis.GetXmin(),fYaxis.GetXmax(),nz,fZaxis.GetXmin(),fZaxis.GetXmax());
   Bool_t computeErrors = kFALSE;
   Bool_t cequalErrors  = kFALSE;
   Bool_t binEntries    = kFALSE;
   if (opt.Contains("b")) binEntries = kTRUE;
   if (opt.Contains("e")) computeErrors = kTRUE;
   if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
   if (computeErrors) h1->Sumw2();
   if (pname != name)  delete [] pname;

   // Fill the projected histogram
   Int_t bin,binx,biny,binz;
   Double_t cont,err;
   for (binx =0;binx<=nx+1;binx++) {
      for (biny =0;biny<=ny+1;biny++) {
         for (binz =0;binz<=nz+1;binz++) {
            bin = GetBin(binx,biny,binz);
            if (binEntries)    cont = GetBinEntries(bin);
            else               cont = GetBinContent(bin);
            err   = GetBinError(bin);
            if (cequalErrors)  h1->SetBinContent(binx,biny,binz, err);
            else               h1->SetBinContent(binx,biny,binz,cont);
            if (computeErrors) h1->SetBinError(binx,biny,binz,err);
         }
      }
   }
   h1->SetEntries(fEntries);
   return h1;
}

//______________________________________________________________________________
void TProfile3D::PutStats(Double_t *stats)
{
   // Replace current statistics with the values in array stats

   TH3::PutStats(stats);
   fTsumwt  = stats[11];
   fTsumwt2 = stats[12];
}

//______________________________________________________________________________
void TProfile3D::Reset(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*Reset contents of a Profile3D histogram*-*-*-*-*-*-*-*
//*-*                =======================================
   TH3D::Reset(option);
   fBinEntries.Reset();
   TString opt = option;
   opt.ToUpper();
   if (opt.Contains("ICE")) return;
   fTsumwt = fTsumwt2 = 0;
}

//______________________________________________________________________________
void TProfile3D::RebinAxis(Double_t x, const char* ax)
{
// Profile histogram is resized along ax such that x is in the axis range.
// The new axis limits are recomputed by doubling iteratively
// the current axis range until the specified value x is within the limits.
// The algorithm makes a copy of the histogram, then loops on all bins
// of the old histogram to fill the rebinned histogram.
// Takes into account errors (Sumw2) if any.
// The bit kCanRebin must be set before invoking this function.
//  Ex:  h->SetBit(TH1::kCanRebin);

   if (!TestBit(kCanRebin)) return;
   char achoice = toupper(ax[0]);
   TAxis *axis = &fXaxis;
   if (achoice == 'Y') axis = &fYaxis;
   if (achoice == 'Z') axis = &fZaxis;
   if (axis->GetXmin() >= axis->GetXmax()) return;
   if (axis->GetNbins() <= 0) return;

   Double_t xmin, xmax;
   if (!FindNewAxisLimits(axis, x, xmin, xmax))
      return;

   //save a copy of this histogram
   TProfile3D *hold = (TProfile3D*)Clone();
   hold->SetDirectory(0);
   //set new axis limits
   axis->SetLimits(xmin,xmax);

   Int_t  nbinsx = fXaxis.GetNbins();
   Int_t  nbinsy = fYaxis.GetNbins();
   Int_t  nbinsz = fZaxis.GetNbins();

   //now loop on all bins and refill
   Reset("ICE"); //reset only Integral, contents and Errors

   Double_t bx,by,bz;
   Int_t ix, iy, iz, binx, biny, binz;
   for (binz=1;binz<=nbinsz;binz++) {
      bz  = hold->GetZaxis()->GetBinCenter(binz);
      iz  = fZaxis.FindFixBin(bz);
      for (biny=1;biny<=nbinsy;biny++) {
         by  = hold->GetYaxis()->GetBinCenter(biny);
         iy  = fYaxis.FindFixBin(by);
         for (binx=1;binx<=nbinsx;binx++) {
            bx = hold->GetXaxis()->GetBinCenter(binx);
            ix  = fXaxis.FindFixBin(bx);

            Int_t sourceBin = hold->GetBin(binx,biny,binz);
            Int_t destinationBin = GetBin(ix,iy,iz);
            AddBinContent(destinationBin, hold->fArray[sourceBin]);
            fBinEntries.fArray[destinationBin] += hold->fBinEntries.fArray[sourceBin];
            fSumw2.fArray[destinationBin] += hold->fSumw2.fArray[sourceBin];
         }
      }
   }
   fTsumwt  = hold->fTsumwt;
   fTsumwt2 = hold->fTsumwt2;
   delete hold;
}

//______________________________________________________________________________
void TProfile3D::SavePrimitive(ostream &out, Option_t *option /*= ""*/)
{
   // Save primitive as a C++ statement(s) on output stream out

   //Note the following restrictions in the code generated:
   // - variable bin size not implemented
   // - SetErrorOption not implemented


   char quote = '"';
   out <<"   "<<endl;
   out <<"   "<<ClassName()<<" *";

   out << GetName() << " = new " << ClassName() << "(" << quote
       << GetName() << quote << "," << quote<< GetTitle() << quote
       << "," << GetXaxis()->GetNbins();
   out << "," << GetXaxis()->GetXmin()
       << "," << GetXaxis()->GetXmax();
   out << "," << GetYaxis()->GetNbins();
   out << "," << GetYaxis()->GetXmin()
       << "," << GetYaxis()->GetXmax();
   out << "," << GetZaxis()->GetNbins();
   out << "," << GetZaxis()->GetXmin()
       << "," << GetZaxis()->GetXmax();
   out << "," << fTmin
       << "," << fTmax;
   out << ");" << endl;


   // save bin entries
   Int_t bin;
   for (bin=0;bin<fNcells;bin++) {
      Double_t bi = GetBinEntries(bin);
      if (bi) {
         out<<"   "<<GetName()<<"->SetBinEntries("<<bin<<","<<bi<<");"<<endl;
      }
   }
   //save bin contents
   for (bin=0;bin<fNcells;bin++) {
      Double_t bc = fArray[bin];
      if (bc) {
         out<<"   "<<GetName()<<"->SetBinContent("<<bin<<","<<bc<<");"<<endl;
      }
   }
   // save bin errors
   if (fSumw2.fN) {
      for (bin=0;bin<fNcells;bin++) {
         Double_t be = TMath::Sqrt(fSumw2.fArray[bin]);
         if (be) {
            out<<"   "<<GetName()<<"->SetBinError("<<bin<<","<<be<<");"<<endl;
         }
      }
   }

   TH1::SavePrimitiveHelp(out, option);
}

//______________________________________________________________________________
void TProfile3D::Scale(Double_t c1)
{
//*-*-*-*-*Multiply this profile2D by a constant c1*-*-*-*-*-*-*-*-*
//*-*      ========================================
//
//   this = c1*this
//
// This function uses the services of TProfile3D::Add
//

   Double_t ent = fEntries;
   fScaling = kTRUE;
   Add(this,this,c1,0);
   fScaling = kFALSE;
   fEntries = ent;
}

//______________________________________________________________________________
void TProfile3D::SetBinEntries(Int_t bin, Double_t w)
{
//*-*-*-*-*-*-*-*-*Set the number of entries in bin*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ================================

   if (bin < 0 || bin >= fNcells) return;
   fBinEntries.fArray[bin] = w;
}

//______________________________________________________________________________
void TProfile3D::SetBins(Int_t nx, Double_t xmin, Double_t xmax, Int_t ny, Double_t ymin, Double_t ymax, Int_t nz, Double_t zmin, Double_t zmax)
{
//*-*-*-*-*-*-*-*-*Redefine  x axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ===========================

   fXaxis.Set(nx,xmin,xmax);
   fYaxis.Set(ny,ymin,ymax);
   fZaxis.Set(ny,zmin,zmax);
   fNcells = (nx+2)*(ny+2)*(nz+2);
   fBinEntries.Set(fNcells);
   fSumw2.Set(fNcells);
}


//______________________________________________________________________________
void TProfile3D::SetBuffer(Int_t buffersize, Option_t *)
{
// set the buffer size in units of 8 bytes (double)

   if (fBuffer) {
      BufferEmpty();
      delete [] fBuffer;
      fBuffer = 0;
   }
   if (buffersize <= 0) {
      fBufferSize = 0;
      return;
   }
   if (buffersize < 100) buffersize = 100;
   fBufferSize = 1 + 5*buffersize; 
   fBuffer = new Double_t[fBufferSize];
   memset(fBuffer,0,8*fBufferSize);
}

//______________________________________________________________________________
void TProfile3D::SetErrorOption(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*Set option to compute profile2D errors*-*-*-*-*-*-*-*
//*-*                =======================================
//
//    The computation of errors is based on the parameter option:
//    option:
//     ' '  (Default) Errors are Spread/SQRT(N) for Spread.ne.0. ,
//                      "     "  SQRT(T)/SQRT(N) for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//     's'            Errors are Spread  for Spread.ne.0. ,
//                      "     "  SQRT(T)  for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//     'i'            Errors are Spread/SQRT(N) for Spread.ne.0. ,
//                      "     "  1./SQRT(12.*N) for Spread.eq.0,N.gt.0 ,
//                      "     "  0.  for N.eq.0
//   See TProfile3D::BuildOptions for explanation of all options

   TString opt = option;
   opt.ToLower();
   fErrorMode = kERRORMEAN;
   if (opt.Contains("s")) fErrorMode = kERRORSPREAD;
   if (opt.Contains("i")) fErrorMode = kERRORSPREADI;
   if (opt.Contains("g")) fErrorMode = kERRORSPREADG;
}