// @(#)root/hist:$Id$
// Author: Rene Brun   29/09/95

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

#include "TProfileHelper.h"

Bool_t TProfile::fgApproximate = kFALSE;

ClassImp(TProfile)

//______________________________________________________________________________
//
//  Profile histograms are used to display the mean
//  value of Y and its error for each bin in X. The displayed error is by default the
//  standard error on the mean (i.e. the standard deviation divided by the sqrt(n) )
//  Profile histograms are in many cases an
//  elegant replacement of two-dimensional histograms : the inter-relation of two
//  measured quantities X and Y can always be visualized by a two-dimensional
//  histogram or scatter-plot; its representation on the line-printer is not particularly
//  satisfactory, except for sparse data. If Y is an unknown (but single-valued)
//  approximate function of X, this function is displayed by a profile 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 bin J.
//
//                                                    2
//      H(J)  =  sum Y                  E(J)  =  sum Y
//      l(J)  =  sum l                  L(J)  =  sum l
//      h(J)  =  H(J)/L(J)                     mean of Y,
//      s(J)  =  sqrt(E(J)/L(J)- h(J)**2)      standard deviation of Y  (e.g. RMS)
//      e(J)  =  s(J)/sqrt(L(J))               standard error on the mean
//
//  The displayed bin content for bin J of a TProfile is always h(J). The corresponding bin error is by default
//  e(J). In case the option "s" is used (in the constructor or by calling TProfile::BuildOptions)
//  the displayed error is  s(J)
//
//  In the special case where s(J) is zero (eg, case of 1 entry only in one bin)
//  the bin error e(J) is computed from the average of the s(J) for all bins if
//  the static function TProfile::Approximate has been called.
//  This simple/crude approximation was suggested in order to keep the bin
//  during a fit operation. But note that this approximation is not the default behaviour.
//   See also TProfile::BuildOptions for other error options and more detailed explanations
//
//           Example of a profile histogram with its graphics output
//{
//  TCanvas *c1 = new TCanvas("c1","Profile histogram example",200,10,700,500);
//  hprof  = new TProfile("hprof","Profile of pz versus px",100,-4,4,0,20);
//  Float_t px, py, pz;
//  for ( Int_t i=0; i<25000; i++) {
//     gRandom->Rannor(px,py);
//     pz = px*px + py*py;
//     hprof->Fill(px,pz,1);
//  }
//  hprof->Draw();
//}
//Begin_Html
/*
 <img src="gif/profile.gif">
 */
//End_Html
//

//______________________________________________________________________________
TProfile::TProfile() : TH1D()
{
   //*-*-*-*-*-*Default constructor for Profile histograms*-*-*-*-*-*-*-*-*
   //*-*        ==========================================

   BuildOptions(0,0,"");
}

//______________________________________________________________________________
TProfile::~TProfile()
{
   //*-*-*-*-*-*Default destructor for Profile histograms*-*-*-*-*-*-*-*-*
   //*-*        =========================================

}

//______________________________________________________________________________
TProfile::TProfile(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup,Option_t *option)
: TH1D(name,title,nbins,xlow,xup)
{
   //*-*-*-*-*-*Normal Constructor for Profile histograms*-*-*-*-*-*-*-*-*-*
   //*-*        ==========================================
   //
   //  The first five parameters are similar to TH1D::TH1D.
   //  All values of y are accepted at filling time.
   //  To fill a profile histogram, one must use TProfile::Fill function.
   //
   //  Note that when filling the profile histogram the function Fill
   //  checks if the variable y is betyween fYmin and fYmax.
   //  If a minimum or maximum value is set for the Y scale before filling,
   //  then all values below ymin or above ymax will be discarded.
   //  Setting the minimum or maximum value for the Y scale before filling
   //  has the same effect as calling the special TProfile constructor below
   //  where ymin and ymax are specified.
   //
   //  H(J) is printed as the channel contents. The errors displayed are s(J) if CHOPT='S'
   //  (spread option), or e(J) if CHOPT=' ' (error on mean).
   //
   //        See TProfile::BuildOptions for explanation of parameters
   //
   // see also comments in the TH1 base class constructors

   BuildOptions(0,0,option);
}

//______________________________________________________________________________
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Float_t *xbins,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
   //*-*-*-*-*-*Constructor for Profile histograms with variable bin size*-*-*-*-*
   //*-*        =========================================================
   //
   //        See TProfile::BuildOptions for more explanations on errors
   //
   // see also comments in the TH1 base class constructors

   BuildOptions(0,0,option);
}

//______________________________________________________________________________
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Double_t *xbins,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
   //*-*-*-*-*-*Constructor for Profile histograms with variable bin size*-*-*-*-*
   //*-*        =========================================================
   //
   //        See TProfile::BuildOptions for more explanations on errors
   //
   // see also comments in the TH1 base class constructors

   BuildOptions(0,0,option);
}

//______________________________________________________________________________
TProfile::TProfile(const char *name,const char *title,Int_t nbins,const Double_t *xbins,Double_t ylow,Double_t yup,Option_t *option)
: TH1D(name,title,nbins,xbins)
{
   //*-*-*-*-*-*Constructor for Profile histograms with variable bin size*-*-*-*-*
   //*-*        =========================================================
   //
   //        See TProfile::BuildOptions for more explanations on errors
   //
   // see also comments in the TH1 base class constructors

   BuildOptions(ylow,yup,option);
}

//______________________________________________________________________________
TProfile::TProfile(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup,Double_t ylow,Double_t yup,Option_t *option)
: TH1D(name,title,nbins,xlow,xup)
{
   //*-*-*-*-*-*Constructor for Profile histograms with range in y*-*-*-*-*-*
   //*-*        ==================================================
   //  The first five parameters are similar to TH1D::TH1D.
   //  Only the values of Y between ylow and yup will be considered at filling time.
   //  ylow and yup will also be the maximum and minimum values
   //  on the y scale when drawing the profile.
   //
   //        See TProfile::BuildOptions for more explanations on errors
   //
   // see also comments in the TH1 base class constructors

   BuildOptions(ylow,yup,option);
}


//______________________________________________________________________________
void TProfile::BuildOptions(Double_t ymin, Double_t ymax, Option_t *option)
{
   //*-*-*-*-*-*-*Set Profile histogram structure and options*-*-*-*-*-*-*-*-*
   //*-*          ===========================================
   //    ymin:  minimum value allowed for y
   //    ymax:  maximum value allowed for y
   //            if (ymin = ymax = 0) there are no limits on the allowed y values (ymin = -inf, ymax = +inf)
   //
   //    option:  this is the option for the computation of the y error of the profile ( TProfile::GetBinError )
   //             possible values for the options are:
   //
   //
   //     ' '  (Default) the bin errors are the standard error on the mean of Y  =  S(Y)/SQRT(N)
   //                    where S(Y) is the standard deviation (RMS) of the Y data in the bin
   //                    and N is the number of bin entries (from TProfile::GetBinEntries(ibin) )
   //                    (i.e the errors are the standard error on the bin content of the profile)
   //
   //     's'            Errors are the standard deviation of Y, S(Y)
   //
   //     'i'            Errors are S(Y)/SQRT(N) (standard error on the mean as in the default)
   //                    The only difference is only when the standard deviation in Y is zero.
   //                    In this  case the error a standard deviation = 1/SQRT(12) is assumed and the error is
   //                    1./SQRT(12*N).
   //                    This approximation assumes that the Y values are integer (e.g. ADC counts)
   //                    and have an implicit uncertainty of y +/- 0.5. With the assumption that the probability that y
   //                    takes any value between y-0.5 and y+0.5 is uniform, its standard error is 1/SQRT(12)
   //
   //     'g'            Errors are 1./SQRT(W) where W is the sum of the weights for the bin J
   //                    W is obtained as from TProfile::GetBinEntries(ibin)
   //                    This errors corresponds to the standard deviation of weighted mean where each
   //                    measurement Y is uncorrelated and has an error sigma, which is expressed in the
   //                    weight used to fill the Profile:  w = 1/sigma^2
   //                    The resulting  error in TProfile is then 1./SQRT( Sum(1./sigma^2) )

   //
   //    In the case of Profile filled weights and with TProfile::Sumw2() called,
   //    STD(Y) is the standard deviation of the weighted sample Y and N is in this case the
   //    number of effective entries (TProfile::GetBinEffectiveEntries(ibin) )
   //
   //    If a bin has N data points all with the same value Y (especially
   //    possible when dealing with integers), the spread in Y for that bin
   //    is zero, and the uncertainty assigned is also zero, and the bin is
   //    ignored in making subsequent fits.
   //    To avoid this problem one can use an approximation for the standard deviation S(Y),
   //    by using the average of all the S(Y) of the other Profile bins. To use this approximation
   //    one must call before TProfile::Approximate
   //    This approximayion applies only for the default and  the 's' options
   //

   SetErrorOption(option);

   // create extra profile data structire (bin entries/ y^2 and sum of weight square)
   TProfileHelper::BuildArray(this);

   fYmin = ymin;
   fYmax = ymax;
   fScaling = kFALSE;
   fTsumwy = fTsumwy2 = 0;

}

//______________________________________________________________________________
TProfile::TProfile(const TProfile &profile) : TH1D()
{
   // Copy constructor.

   ((TProfile&)profile).Copy(*this);
}


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

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


//______________________________________________________________________________
Bool_t TProfile::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 kFALSE;
   }
   if (!h1->InheritsFrom(TProfile::Class())) {
      Error("Add","Attempt to add a non-profile object");
      return kFALSE;
   }

   return TProfileHelper::Add(this, this, h1, 1, c1);
}

//______________________________________________________________________________
Bool_t TProfile::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
{
   //*-*-*-*-*Replace contents of this profile by the addition of h1 and h2*-*-*
   //*-*      =============================================================
   //
   //   this = c1*h1 + c2*h2
   //
   //   c1 and c2 are considered as weights applied to the two summed profiles.
   //   The operation acts therefore like merging the two profiles with a weight c1 and c2
   //

   if (!h1 || !h2) {
      Error("Add","Attempt to add a non-existing profile");
      return kFALSE;
   }
   if (!h1->InheritsFrom(TProfile::Class())) {
      Error("Add","Attempt to add a non-profile object");
      return kFALSE;
   }
   if (!h2->InheritsFrom(TProfile::Class())) {
      Error("Add","Attempt to add a non-profile object");
      return kFALSE;
   }
   return TProfileHelper::Add(this, h1, h2, c1, c2);
}


//______________________________________________________________________________
void TProfile::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 profile is less than 1002
   //  - the bin number of entries is small ( <5)
   //  - the estimated bin error is extremely small compared to the bin content
   //  (see TProfile::GetBinError)

   fgApproximate = approx;
}

//______________________________________________________________________________
Int_t TProfile::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("ICES"); // reset without deleting the functions
      fBuffer = buffer;
   }
   if (CanExtendAllAxes() || fXaxis.GetXmax() <= fXaxis.GetXmin()) {
      //find min, max of entries in buffer
      Double_t xmin = fBuffer[2];
      Double_t xmax = xmin;
      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;
      }
      if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
         THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax);
      } else {
         fBuffer = 0;
         Int_t keep = fBufferSize; fBufferSize = 0;
         if (xmin <  fXaxis.GetXmin()) ExtendAxis(xmin,&fXaxis);
         if (xmax >= fXaxis.GetXmax()) ExtendAxis(xmax,&fXaxis);
         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 TProfile::BufferFill(Double_t x, Double_t y, 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

   if (!fBuffer) return -2;
   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("ICES");  // reset without deleting the functions
         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 -2;
}

//______________________________________________________________________________
void TProfile::Copy(TObject &obj) const
{
   //*-*-*-*-*-*-*-*Copy a Profile histogram to a new profile histogram*-*-*-*-*
   //*-*            ===================================================

   try {
      TProfile & pobj = dynamic_cast<TProfile&>(obj);
      TH1D::Copy(pobj);
      fBinEntries.Copy(pobj.fBinEntries);
      fBinSumw2.Copy(pobj.fBinSumw2);
      for (int bin=0;bin<fNcells;bin++) {
         pobj.fArray[bin]        = fArray[bin];
         pobj.fSumw2.fArray[bin] = fSumw2.fArray[bin];
      }

      pobj.fYmin = fYmin;
      pobj.fYmax = fYmax;
      pobj.fScaling   = fScaling;
      pobj.fErrorMode = fErrorMode;
      pobj.fTsumwy    = fTsumwy;
      pobj.fTsumwy2   = fTsumwy2;

   } catch(...) {
      Fatal("Copy","Cannot copy a TProfile in a %s",obj.IsA()->GetName());
   }

}


//______________________________________________________________________________
Bool_t TProfile::Divide(TF1 *, Double_t )
{
   // Performs the operation: this = this/(c1*f1)
   // This function is not implemented for the TProfile

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

//______________________________________________________________________________
Bool_t TProfile::Divide(const TH1 *h1)
{
   //*-*-*-*-*-*-*-*-*-*-*Divide this profile by h1*-*-*-*-*-*-*-*-*-*-*-*-*
   //*-*                  =========================
   //
   //   this = this/h1
   // This function accepts to divide a TProfile by a histogram
   //
   // The function return kFALSE if the divide operation failed

   if (!h1) {
      Error("Divide","Attempt to divide a non-existing profile");
      return kFALSE;
   }
   if (!h1->InheritsFrom(TH1::Class())) {
      Error("Divide","Attempt to divide by a non-profile or non-histogram object");
      return kFALSE;
   }
   TProfile *p1 = (TProfile*)h1;

   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);


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

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

   //*-*- Loop on bins (including underflows/overflows)
   Int_t bin;
   Double_t *cu1=0, *er1=0, *en1=0;
   Double_t e0,e1,c12;
   if (h1->InheritsFrom(TProfile::Class())) {
      cu1 = p1->GetW();
      er1 = p1->GetW2();
      en1 = p1->GetB();
   }
   Double_t c0,c1,w,z,x;
   for (bin=0;bin<=nbinsx+1;bin++) {
      c0  = fArray[bin];
      if (cu1) c1  = cu1[bin];
      else     c1  = h1->GetBinContent(bin);
      if (c1) w = c0/c1;
      else    w = 0;
      fArray[bin] = w;
      z = TMath::Abs(w);
      x = fXaxis.GetBinCenter(bin);
      fEntries++;
      fTsumw   += z;
      fTsumw2  += z*z;
      fTsumwx  += z*x;
      fTsumwx2 += z*x*x;
      fTsumwy  += z*c1;
      fTsumwx2 += z*c1*c1;
      e0 = fSumw2.fArray[bin];
      if (er1) e1 = er1[bin];
      else    {e1 = h1->GetBinError(bin); e1*=e1;}
      c12= c1*c1;
      if (!c1) fSumw2.fArray[bin] = 0;
      else     fSumw2.fArray[bin] = (e0*c1*c1 + e1*c0*c0)/(c12*c12);
      if (!en1) continue;
      if (!en1[bin]) fBinEntries.fArray[bin] = 0;
      else           fBinEntries.fArray[bin] /= en1[bin];
   }
   // mantaining the correct sum of weights square is not supported when dividing
   // bin error resulting from division of profile needs to be checked
   if (fBinSumw2.fN) {
      Warning("Divide","Cannot preserve during the division of profiles the sum of bin weight square");
      fBinSumw2 = TArrayD();
   }

   return kTRUE;
}


//______________________________________________________________________________
Bool_t TProfile::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
   //*-*-*-*-*Replace contents of this profile by the division of h1 by h2*-*-*
   //*-*      ============================================================
   //
   //   this = c1*h1/(c2*h2)
   //
   // The function return kFALSE if the divide operation failed

   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 profile");
      return kFALSE;
   }
   if (!h1->InheritsFrom(TProfile::Class())) {
      Error("Divide","Attempt to divide a non-profile object");
      return kFALSE;
   }
   TProfile *p1 = (TProfile*)h1;
   if (!h2->InheritsFrom(TProfile::Class())) {
      Error("Divide","Attempt to divide by a non-profile object");
      return kFALSE;
   }
   TProfile *p2 = (TProfile*)h2;

   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);

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

   //THE ALGORITHM COMPUTING THE ERRORS IS WRONG. HELP REQUIRED
   printf("WARNING!!: The algorithm in TProfile::Divide computing the errors is not accurate\n");
   printf(" Instead of Divide(TProfile *h1, TProfile *h2, do:\n");
   printf("   TH1D *p1 = h1->ProjectionX();\n");
   printf("   TH1D *p2 = h2->ProjectionX();\n");
   printf("   p1->Divide(p2);\n");

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

   //*-*- Loop on bins (including underflows/overflows)
   Int_t bin;
   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,z,x,ac1,ac2;
   //d1 = c1*c1;
   //d2 = c2*c2;
   ac1 = TMath::Abs(c1);
   ac2 = TMath::Abs(c2);
   for (bin=0;bin<=nbinsx+1;bin++) {
      b1  = cu1[bin];
      b2  = cu2[bin];
      if (b2) w = c1*b1/(c2*b2);
      else    w = 0;
      fArray[bin] = w;
      z = TMath::Abs(w);
      x = fXaxis.GetBinCenter(bin);
      fEntries++;
      fTsumw   += z;
      fTsumw2  += z*z;
      fTsumwx  += z*x;
      fTsumwx2 += z*x*x;
      //fTsumwy  += z*x;
      //fTsumwy2 += z*x*x;
      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)/b2);
         } else {
            //fSumw2.fArray[bin] = d1*d2*(e1*b2*b2 + e2*b1*b1)/(b22*b22);
            fSumw2.fArray[bin] = ac1*ac2*(e1*b2*b2 + e2*b1*b1)/(b22*b22);
         }
      }
      if (en2[bin]) fBinEntries.fArray[bin] = en1[bin]/en2[bin];
      else          fBinEntries.fArray[bin] = 0;
   }

   // mantaining the correct sum of weights square is not supported when dividing
   // bin error resulting from division of profile needs to be checked
   if (fBinSumw2.fN) {
      Warning("Divide","Cannot preserve during the division of profiles the sum of bin weight square");
      fBinSumw2 = TArrayD();
   }

   return kTRUE;
}

//______________________________________________________________________________
Int_t TProfile::Fill(Double_t x, Double_t y)
{
   //*-*-*-*-*-*-*-*-*-*-*Fill a Profile histogram (no weights)*-*-*-*-*-*-*-*
   //*-*                  =====================================

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

   Int_t bin;
   if (fYmin != fYmax) {
      if (y <fYmin || y> fYmax || TMath::IsNaN(y) ) return -1;
   }

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

//______________________________________________________________________________
Int_t TProfile::Fill(const char *namex, Double_t y)
{
   // Fill a Profile histogram (no weights)
   //
   Int_t bin;
   if (fYmin != fYmax) {
      if (y <fYmin || y> fYmax || TMath::IsNaN(y) ) return -1;
   }

   fEntries++;
   bin =fXaxis.FindBin(namex);
   AddBinContent(bin, y);
   fSumw2.fArray[bin] += (Double_t)y*y;
   fBinEntries.fArray[bin] += 1;
   if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += 1;
   if (bin == 0 || bin > fXaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   Double_t x = fXaxis.GetBinCenter(bin);
   fTsumw++;
   fTsumw2++;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   fTsumwy  += y;
   fTsumwy2 += y*y;
   return bin;
}

//______________________________________________________________________________
Int_t TProfile::Fill(Double_t x, Double_t y, Double_t w)
{
   //*-*-*-*-*-*-*-*-*-*-*Fill a Profile histogram with weights*-*-*-*-*-*-*-*
   //*-*                  =====================================

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

   Int_t bin;
   if (fYmin != fYmax) {
      if (y <fYmin || y> fYmax || TMath::IsNaN(y) ) return -1;
   }

   Double_t u= w;
   fEntries++;
   bin =fXaxis.FindBin(x);
   AddBinContent(bin, u*y);
   fSumw2.fArray[bin] += u*y*y;
   if (!fBinSumw2.fN && u != 1.0 && !TestBit(TH1::kIsNotW))  Sumw2();  // must be called before accumulating the entries
   if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += u*u;
   fBinEntries.fArray[bin] += u;
   if (bin == 0 || bin > fXaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   fTsumw   += u;
   fTsumw2  += u*u;
   fTsumwx  += u*x;
   fTsumwx2 += u*x*x;
   fTsumwy  += u*y;
   fTsumwy2 += u*y*y;
   return bin;
}

//______________________________________________________________________________
Int_t TProfile::Fill(const char *namex, Double_t y, Double_t w)
{
   // Fill a Profile histogram with weights
   //
   Int_t bin;

   if (fYmin != fYmax) {
      if (y <fYmin || y> fYmax || TMath::IsNaN(y) ) return -1;
   }

   Double_t u= w; // (w > 0 ? w : -w);
   fEntries++;
   bin =fXaxis.FindBin(namex);
   AddBinContent(bin, u*y);
   fSumw2.fArray[bin] += u*y*y;
   if (!fBinSumw2.fN && u != 1.0 && !TestBit(TH1::kIsNotW))  Sumw2();  // must be called before accumulating the entries
   if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += u*u;
   fBinEntries.fArray[bin] += u;
   if (bin == 0 || bin > fXaxis.GetNbins()) {
      if (!fgStatOverflows) return -1;
   }
   Double_t x = fXaxis.GetBinCenter(bin);
   fTsumw   += u;
   fTsumw2  += u*u;
   fTsumwx  += u*x;
   fTsumwx2 += u*x*x;
   fTsumwy  += u*y;
   fTsumwy2 += u*y*y;
   return bin;
}


//______________________________________________________________________________
void TProfile::FillN(Int_t ntimes, const Double_t *x, const Double_t *y, const Double_t *w, Int_t stride)
{
   //*-*-*-*-*-*-*-*-*-*-*Fill a Profile histogram with weights*-*-*-*-*-*-*-*
   //*-*                  =====================================
   Int_t bin,i;
   ntimes *= stride;
   Int_t ifirst = 0; 
   //If a buffer is activated, fill buffer
   // (note that this function must not be called from TH2::BufferEmpty)
   if (fBuffer) {
      for (i=0;i<ntimes;i+=stride) {
         if (!fBuffer) break; // buffer can be deleted in BufferFill when is empty
         if (w) BufferFill(x[i],y[i],w[i]);
         else BufferFill(x[i], y[i], 1.);
      }
      // fill the remaining entries if the buffer has been deleted 
      if (i < ntimes && fBuffer==0) 
         ifirst = i;  // start from i 
      else
         return;
   }

   for (i=ifirst;i<ntimes;i+=stride) {
      if (fYmin != fYmax) {
         if (y[i] <fYmin || y[i]> fYmax || TMath::IsNaN(y[i])) continue;
      }

      Double_t u = (w) ? w[i] : 1; // (w[i] > 0 ? w[i] : -w[i]);
      fEntries++;
      bin =fXaxis.FindBin(x[i]);
      AddBinContent(bin, u*y[i]);
      fSumw2.fArray[bin] += u*y[i]*y[i];
      if (!fBinSumw2.fN && u != 1.0 && !TestBit(TH1::kIsNotW))  Sumw2();  // must be called before accumulating the entries
      if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += u*u;
      fBinEntries.fArray[bin] += u;
      if (bin == 0 || bin > fXaxis.GetNbins()) {
         if (!fgStatOverflows) continue;
      }
      fTsumw   += u;
      fTsumw2  += u*u;
      fTsumwx  += u*x[i];
      fTsumwx2 += u*x[i]*x[i];
      fTsumwy  += u*y[i];
      fTsumwy2 += u*y[i]*y[i];
   }
}

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

   if (fBuffer) ((TProfile*)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 TProfile::GetBinEntries(Int_t bin) const
{
   //*-*-*-*-*-*-*Return bin entries of a Profile histogram*-*-*-*-*-*-*-*-*-*
   //*-*          =========================================

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

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

//______________________________________________________________________________
Double_t TProfile::GetBinEffectiveEntries(Int_t bin) const
{
   //            Return bin effective entries for a weighted filled Profile histogram.
   //            In case of an unweighted profile, it is equivalent to the number of entries per bin
   //            The effective entries is defined as the square of the sum of the weights divided by the
   //            sum of the weights square.
   //            TProfile::Sumw2() must be called before filling the profile with weights.
   //            Only by calling this method the  sum of the square of the weights per bin is stored.
   //
   //*-*          =========================================

   return TProfileHelper::GetBinEffectiveEntries((TProfile*)this, bin);
}

//______________________________________________________________________________
Double_t TProfile::GetBinError(Int_t bin) const
{
   // *-*-*-*-*-*-*Return bin error of a Profile histogram*-*-*-*-*-*-*-*-*-*
   // *-*          =======================================
   //
   // Computing errors: A moving field
   // =================================
   // The computation of errors for a TProfile has evolved with the versions
   // of ROOT. The difficulty is in computing errors for bins with low statistics.
   // - prior to version 3.00, 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.
   // - in version 3.00 (18 Dec 2000), the algorithm is protected for values of
   //   eprim2 very small and the bin errors set to the average bin errors, following
   //   recommendations from a group of users.
   // - in version 3.01 (19 Apr 2001), it is realized that the algorithm above
   //   should be applied only to low statistic bins.
   // - in version 3.02 (26 Sep 2001), the same group of users recommend instead
   //   to take two times the average error on all bins for these low
   //   statistics bins giving a very small value for eprim2.
   // - in version 3.04 (Nov 2002), the algorithm is modified/protected for the case
   //   when a TProfile is projected (ProjectionX). The previous algorithm
   //   generated a N^2 problem when projecting a TProfile with a large number of
   //   bins (eg 100000).
   // - in version 3.05/06, a new static function TProfile::Approximate
   //   is introduced to enable or disable (default) the approximation.
   //
   // Ideas for improvements of this algorithm are welcome. No suggestions
   // received since our call for advice to roottalk in Jul 2002.
   // see for instance: http://root.cern.ch/root/roottalk/roottalk02/2916.html

   return TProfileHelper::GetBinError((TProfile*)this, bin);
}

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

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

//______________________________________________________________________________
void TProfile::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
   //
   // 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) ((TProfile*)this)->BufferEmpty();

   // Loop on bins
   Int_t bin, binx;
   if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange)) {
      for (bin=0;bin<6;bin++) stats[bin] = 0;
      if (!fBinEntries.fArray) return;
      Int_t firstBinX = fXaxis.GetFirst();
      Int_t lastBinX  = fXaxis.GetLast();
      // include underflow/overflow if TH1::StatOverflows(kTRUE) in case no range is set on the axis
      if (fgStatOverflows && !fXaxis.TestBit(TAxis::kAxisRange)) {
         if (firstBinX == 1) firstBinX = 0;
         if (lastBinX ==  fXaxis.GetNbins() ) lastBinX += 1;
      }
      for (binx = firstBinX; binx <= lastBinX; binx++) {
         Double_t w   = fBinEntries.fArray[binx];
         Double_t w2  = (fBinSumw2.fN ? fBinSumw2.fArray[binx] : w);
         Double_t x   = fXaxis.GetBinCenter(binx);
         stats[0] += w;
         stats[1] += w2;
         stats[2] += w*x;
         stats[3] += w*x*x;
         stats[4] += fArray[binx];
         stats[5] += fSumw2.fArray[binx];
      }
   } else {
      if (fTsumwy == 0 && fTsumwy2 == 0) {
         //this case may happen when processing TProfiles with version <=3
         TProfile *p = (TProfile*)this; // chheting with const
         for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
            p->fTsumwy  += fArray[binx];
            p->fTsumwy2 += fSumw2.fArray[binx];
         }
      }
      stats[0] = fTsumw;
      stats[1] = fTsumw2;
      stats[2] = fTsumwx;
      stats[3] = fTsumwx2;
      stats[4] = fTsumwy;
      stats[5] = fTsumwy2;
   }
}

//___________________________________________________________________________
void TProfile::LabelsDeflate(Option_t *option)
{
   // Reduce the number of bins for this axis to the number of bins having a label.

   TProfileHelper::LabelsDeflate(this, option);
}

//___________________________________________________________________________
void TProfile::LabelsInflate(Option_t *options)
{
   // Double the number of bins for axis.
   // Refill histogram
   // This function is called by TAxis::FindBin(const char *label)

   TProfileHelper::LabelsInflate(this, options);
}

//___________________________________________________________________________
void TProfile::LabelsOption(Option_t *option, Option_t * /*ax */)
{
   //  Set option(s) to draw axis with labels
   //  option = "a" sort by alphabetic order
   //         = ">" sort by decreasing values
   //         = "<" sort by increasing values
   //         = "h" draw labels horizonthal
   //         = "v" draw labels vertical
   //         = "u" draw labels up (end of label right adjusted)
   //         = "d" draw labels down (start of label left adjusted)

   THashList *labels = fXaxis.GetLabels();
   if (!labels) {
      Warning("LabelsOption","Cannot sort. No labels");
      return;
   }
   TString opt = option;
   opt.ToLower();
   if (opt.Contains("h")) {
      fXaxis.SetBit(TAxis::kLabelsHori);
      fXaxis.ResetBit(TAxis::kLabelsVert);
      fXaxis.ResetBit(TAxis::kLabelsDown);
      fXaxis.ResetBit(TAxis::kLabelsUp);
   }
   if (opt.Contains("v")) {
      fXaxis.SetBit(TAxis::kLabelsVert);
      fXaxis.ResetBit(TAxis::kLabelsHori);
      fXaxis.ResetBit(TAxis::kLabelsDown);
      fXaxis.ResetBit(TAxis::kLabelsUp);
   }
   if (opt.Contains("u")) {
      fXaxis.SetBit(TAxis::kLabelsUp);
      fXaxis.ResetBit(TAxis::kLabelsVert);
      fXaxis.ResetBit(TAxis::kLabelsDown);
      fXaxis.ResetBit(TAxis::kLabelsHori);
   }
   if (opt.Contains("d")) {
      fXaxis.SetBit(TAxis::kLabelsDown);
      fXaxis.ResetBit(TAxis::kLabelsVert);
      fXaxis.ResetBit(TAxis::kLabelsHori);
      fXaxis.ResetBit(TAxis::kLabelsUp);
   }
   Int_t sort = -1;
   if (opt.Contains("a")) sort = 0;
   if (opt.Contains(">")) sort = 1;
   if (opt.Contains("<")) sort = 2;
   if (sort < 0) return;

   Int_t n = TMath::Min(fXaxis.GetNbins(), labels->GetSize());
   Int_t *a = new Int_t[n+2];
   Int_t i,j;
   Double_t *cont   = new Double_t[n+2];
   Double_t *sumw   = new Double_t[n+2];
   Double_t *errors = new Double_t[n+2];
   Double_t *ent    = new Double_t[n+2];
   THashList *labold = new THashList(labels->GetSize(),1);
   TIter nextold(labels);
   TObject *obj;
   while ((obj=nextold())) {
      labold->Add(obj);
   }
   labels->Clear();
   if (sort > 0) {
      //---sort by values of bins
      for (i=1;i<=n;i++) {
         sumw[i-1]   = fArray[i];
         errors[i-1] = fSumw2.fArray[i];
         ent[i-1]    = fBinEntries.fArray[i];
         if (fBinEntries.fArray[i] == 0) cont[i-1] = 0;
         else cont[i-1] = fArray[i]/fBinEntries.fArray[i];
      }
      if (sort ==1) TMath::Sort(n,cont,a,kTRUE);  //sort by decreasing values
      else          TMath::Sort(n,cont,a,kFALSE); //sort by increasing values
      for (i=1;i<=n;i++) {
         fArray[i] = sumw[a[i-1]];
         fSumw2.fArray[i] = errors[a[i-1]];
         fBinEntries.fArray[i] = ent[a[i-1]];
      }
      for (i=1;i<=n;i++) {
         obj = labold->At(a[i-1]);
         labels->Add(obj);
         obj->SetUniqueID(i);
      }
   } else {
      //---alphabetic sort
      const UInt_t kUsed = 1<<18;
      TObject *objk=0;
      a[0] = 0;
      a[n+1] = n+1;
      for (i=1;i<=n;i++) {
         const char *label = "zzzzzzzzzzzz";
         for (j=1;j<=n;j++) {
            obj = labold->At(j-1);
            if (!obj) continue;
            if (obj->TestBit(kUsed)) continue;
            //use strcasecmp for case non-sensitive sort (may be an option)
            if (strcmp(label,obj->GetName()) < 0) continue;
            objk = obj;
            a[i] = j;
            label = obj->GetName();
         }
         if (objk) {
            objk->SetUniqueID(i);
            labels->Add(objk);
            objk->SetBit(kUsed);
         }
      }
      for (i=1;i<=n;i++) {
         obj = labels->At(i-1);
         if (!obj) continue;
         obj->ResetBit(kUsed);
      }

      for (i=1;i<=n;i++) {
         sumw[i]   = fArray[a[i]];
         errors[i] = fSumw2.fArray[a[i]];
         ent[i]    = fBinEntries.fArray[a[i]];
      }
      for (i=1;i<=n;i++) {
         fArray[i] = sumw[i];
         fSumw2.fArray[i] = errors[i];
         fBinEntries.fArray[i] = ent[i];
      }
   }
   delete labold;
   if (a)      delete [] a;
   if (sumw)   delete [] sumw;
   if (cont)   delete [] cont;
   if (errors) delete [] errors;
   if (ent)    delete [] ent;
}

//______________________________________________________________________________
Long64_t TProfile::Merge(TCollection *li)
{
   //Merge all histograms in the collection in this histogram.
   //This function computes the min/max for the x axis,
   //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 axis x 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.

   return TProfileHelper::Merge(this, li);
}


//______________________________________________________________________________
Bool_t TProfile::Multiply(TF1 *f1, Double_t c1)
{
   // Performs the operation: this = this*c1*f1
   //
   // The function return kFALSE if the Multiply operation failed


   if (!f1) {
      Error("Multiply","Attempt to multiply by a null function");
      return kFALSE;
   }

   Int_t nbinsx = GetNbinsX();

   //*-*- Add statistics
   Double_t xx[1], cf1, ac1 = TMath::Abs(c1);
   Double_t s1[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = 0;}
   PutStats(s1);

   SetMinimum();
   SetMaximum();

   //*-*- Loop on bins (including underflows/overflows)
   Int_t bin;
   for (bin=0;bin<=nbinsx+1;bin++) {
      xx[0] = fXaxis.GetBinCenter(bin);
      if (!f1->IsInside(xx)) continue;
      TF1::RejectPoint(kFALSE);
      cf1 = f1->EvalPar(xx);
      if (TF1::RejectedPoint()) continue;
      fArray[bin]             *= c1*cf1;
      //see http://savannah.cern.ch/bugs/?func=detailitem&item_id=14851
      //fSumw2.fArray[bin]      *= c1*c1*cf1*cf1;
      fSumw2.fArray[bin]      *= ac1*cf1*cf1;
      //fBinEntries.fArray[bin] *= ac1*TMath::Abs(cf1);
   }
   return kTRUE;
}

//______________________________________________________________________________
Bool_t TProfile::Multiply(const TH1 *)
{
   //*-*-*-*-*-*-*-*-*-*-*Multiply this profile by h1*-*-*-*-*-*-*-*-*-*-*-*-*
   //*-*                  =============================
   //
   //   this = this*h1
   //
   Error("Multiply","Multiplication of profile histograms not implemented");
   return kFALSE;
}


//______________________________________________________________________________
Bool_t TProfile::Multiply(const TH1 *, const TH1 *, Double_t, Double_t, Option_t *)
{
   //*-*-*-*-*Replace contents of this profile by multiplication of h1 by h2*-*
   //*-*      ================================================================
   //
   //   this = (c1*h1)*(c2*h2)
   //
   Error("Multiply","Multiplication of profile histograms not implemented");
   return kFALSE;
}

//______________________________________________________________________________
TH1D *TProfile::ProjectionX(const char *name, Option_t *option) const
{
   //*-*-*-*-*Project this profile into a 1-D histogram along X*-*-*-*-*-*-*
   //*-*      =================================================
   //
   //   The projection is always of the type TH1D.
   //
   //   if option "E" is specified the errors of the projected histogram are computed and set
   //      to be equal to the errors of the profile.
   //      Option "E" is defined as the default one in the header file.
   //   if option "" is specified the histogram errors are simply the sqrt of its content
   //   if option "B" is specified, the content of bin of the returned histogram
   //      will be equal to the GetBinEntries(bin) of the profile,
   //      otherwise (default) it will be equal to GetBinContent(bin)
   //   if option "C=E" the bin contents of the projection are set to the
   //       bin errors of the profile
   //   if option "W" is specified the bin content of the projected histogram  is set to the
   //       product of the bin content of the profile and the entries.
   //       With this option the returned histogram will be equivalent to the one obtained by
   //       filling directly a TH1D using the 2-nd value as a weight.
   //       This makes sense only for profile filled with weights =1. If not, the error of the
   //        projected histogram obtained with this option will not be correct.


   TString opt = option;
   opt.ToLower();
   Int_t nx = fXaxis.GetNbins();

   // Create the projection histogram
   TString pname = name;
   if (pname == "_px") {
      pname = GetName();
      pname.Append("_px");
   }
   TH1D *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);
   }
   Bool_t computeErrors = kFALSE;
   Bool_t cequalErrors  = kFALSE;
   Bool_t binEntries    = kFALSE;
   Bool_t binWeight     = kFALSE;
   if (opt.Contains("b")) binEntries = kTRUE;
   if (opt.Contains("e")) computeErrors = kTRUE;
   if (opt.Contains("w")) binWeight = kTRUE;
   if (opt.Contains("c=e")) {cequalErrors = kTRUE; computeErrors=kFALSE;}
   if (computeErrors || binWeight || (binEntries && fBinSumw2.fN) ) h1->Sumw2();

   // Fill the projected histogram
   Double_t cont;
   for (Int_t bin =0;bin<=nx+1;bin++) {

      if (binEntries)         cont = GetBinEntries(bin);
      else if (cequalErrors)  cont = GetBinError(bin);
      else if (binWeight)     cont = fArray[bin];  // bin content * bin entries
      else                    cont = GetBinContent(bin);    // default case

      h1->SetBinContent(bin ,cont);

      // if option E projected histogram errors are same as profile
      if (computeErrors ) h1->SetBinError(bin , GetBinError(bin) );
      // in case of option W bin error is deduced from bin sum of z**2 values of profile
      // this is correct only if the profile is filled with weights =1
      if (binWeight) h1->GetSumw2()->fArray[bin] = fSumw2.fArray[bin];
      // in case of bin entries and profile is weighted, we need to set also the bin error
      if (binEntries && fBinSumw2.fN ) {
         R__ASSERT(  h1->GetSumw2() );
         h1->GetSumw2()->fArray[bin] =  fBinSumw2.fArray[bin];
      }

   }

   // Copy the axis attributes and the axis labels if needed.
   h1->GetXaxis()->ImportAttributes(this->GetXaxis());
   h1->GetYaxis()->ImportAttributes(this->GetYaxis());
   THashList* labels=this->GetXaxis()->GetLabels();
   if (labels) {
      TIter iL(labels);
      TObjString* lb;
      Int_t i = 1;
      while ((lb=(TObjString*)iL())) {
         h1->GetXaxis()->SetBinLabel(i,lb->String().Data());
         i++;
      }
   }

   h1->SetEntries(fEntries);
   return h1;
}

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

   fTsumw   = stats[0];
   fTsumw2  = stats[1];
   fTsumwx  = stats[2];
   fTsumwx2 = stats[3];
   fTsumwy  = stats[4];
   fTsumwy2 = stats[5];
}

//______________________________________________________________________________
TH1 *TProfile::Rebin(Int_t ngroup, const char*newname, const Double_t *xbins)
{
   //*-*-*-*-*Rebin this profile grouping ngroup bins together*-*-*-*-*-*-*-*-*
   //*-*      ================================================
   //  -case 1  xbins=0
   //   if newname is not blank a new temporary profile hnew is created.
   //   else the current profile is modified (default)
   //   The parameter ngroup indicates how many bins of this have to me merged
   //   into one bin of hnew
   //   If the original profile has errors stored (via Sumw2), the resulting
   //   profile has new errors correctly calculated.
   //
   //   examples: if hp is an existing TProfile histogram with 100 bins
   //     hp->Rebin();  //merges two bins in one in hp: previous contents of hp are lost
   //     hp->Rebin(5); //merges five bins in one in hp
   //     TProfile *hnew = hp->Rebin(5,"hnew"); // creates a new profile hnew
   //                                       //merging 5 bins of hp in one bin
   //
   //   NOTE:  If ngroup is not an exact divider of the number of bins,
   //          the top limit of the rebinned profile is changed
   //          to the upper edge of the bin=newbins*ngroup and the corresponding
   //          bins are added to the overflow bin.
   //          Statistics will be recomputed from the new bin contents.
   //
   //  -case 2  xbins!=0
   //   a new profile is created (you should specify newname).
   //   The parameter ngroup is the number of variable size bins in the created profile
   //   The array xbins must contain ngroup+1 elements that represent the low-edge
   //   of the bins.
   //   The data of the old bins are added to the new bin which contains the bin center
   //   of the old bins. It is possible that information from the old binning are attached
   //   to the under-/overflow bins of the new binning.
   //
   //   examples: if hp is an existing TProfile with 100 bins
   //     Double_t xbins[25] = {...} array of low-edges (xbins[25] is the upper edge of last bin
   //     hp->Rebin(24,"hpnew",xbins);  //creates a new variable bin size profile hpnew

   Int_t nbins    = fXaxis.GetNbins();
   Double_t xmin  = fXaxis.GetXmin();
   Double_t xmax  = fXaxis.GetXmax();
   if ((ngroup <= 0) || (ngroup > nbins)) {
      Error("Rebin", "Illegal value of ngroup=%d",ngroup);
      return 0;
   }
   if (!newname && xbins) {
      Error("Rebin","if xbins is specified, newname must be given");
      return 0;
   }

   Int_t newbins = nbins/ngroup;
   if (!xbins) {
      Int_t nbg = nbins/ngroup;
      if (nbg*ngroup != nbins) {
         Warning("Rebin", "ngroup=%d must be an exact divider of nbins=%d",ngroup,nbins);
      }
   }
   else {
      // in the case of xbins given (rebinning in variable bins) ngroup is the new number of bins.
      //  and number of grouped bins is not constant.
      // when looping for setting the contents for the new histogram we
      // need to loop on all bins of original histogram. Set then ngroup=nbins
      newbins = ngroup;
      ngroup = nbins;
   }

   // Save old bin contents into a new array
   Double_t *oldBins   = new Double_t[nbins+2];
   Double_t *oldCount  = new Double_t[nbins+2];
   Double_t *oldErrors = new Double_t[nbins+2];
   Double_t *oldBinw2  = (fBinSumw2.fN ? new Double_t[nbins+2] : 0  );
   Int_t bin, i;
   Double_t *cu1 = GetW();
   Double_t *er1 = GetW2();
   Double_t *en1 = GetB();
   Double_t *ew1 = GetB2();

   for (bin=0;bin<=nbins+1;bin++) {
      oldBins[bin]   = cu1[bin];
      oldCount[bin]  = en1[bin];
      oldErrors[bin] = er1[bin];
      if (ew1 && fBinSumw2.fN) oldBinw2[bin]  = ew1[bin];
   }

   // create a clone of the old histogram if newname is specified
   TProfile *hnew = this;
   if ((newname && strlen(newname) > 0) || xbins) {
      hnew = (TProfile*)Clone(newname);
   }

   // in case of ngroup not an excat divider of nbins,
   // top limit is changed (see NOTE in method comment)
   if(!xbins && (newbins*ngroup != nbins)) {
      xmax = fXaxis.GetBinUpEdge(newbins*ngroup);
      hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin
   }

   // set correctly the axis and resizes the bin arrays
   if(!xbins && (fXaxis.GetXbins()->GetSize() > 0)){
      // for rebinning of variable bins in a constant group
      Double_t *bins = new Double_t[newbins+1];
      for(i = 0; i <= newbins; ++i) bins[i] = fXaxis.GetBinLowEdge(1+i*ngroup);
      hnew->SetBins(newbins,bins); //this also changes the bin array's
      delete [] bins;
   } else if (xbins) {
      // when rebinning in variable bins
      hnew->SetBins(newbins,xbins);
   } else {
      hnew->SetBins(newbins,xmin,xmax);
   }

   // merge bin contents ignoring now underflow/overflows
   if (fBinSumw2.fN) hnew->Sumw2();

   // Start merging only once the new lowest edge is reached
   Int_t startbin = 1;
   const Double_t newxmin = hnew->GetXaxis()->GetBinLowEdge(1);
   while( fXaxis.GetBinCenter(startbin) < newxmin && startbin <= nbins ) {
      startbin++;
   }

   Double_t *cu2 = hnew->GetW();
   Double_t *er2 = hnew->GetW2();
   Double_t *en2 = hnew->GetB();
   Double_t *ew2 = hnew->GetB2();
   Int_t oldbin = startbin;
   Double_t binContent, binCount, binError, binSumw2;
   for (bin = 1;bin<=newbins;bin++) {
      binContent = 0;
      binCount   = 0;
      binError   = 0;
      binSumw2   = 0;

      //for xbins != 0: ngroup == nbins
      Int_t imax = ngroup;
      Double_t xbinmax = hnew->GetXaxis()->GetBinUpEdge(bin);
      for (i=0;i<ngroup;i++) {
         if((hnew == this && (oldbin+i > nbins)) ||
            (hnew != this && (fXaxis.GetBinCenter(oldbin+i) > xbinmax)))
         {
            imax = i;
            break;
         }

         binContent += oldBins[oldbin+i];
         binCount   += oldCount[oldbin+i];
         binError   += oldErrors[oldbin+i];
         if (fBinSumw2.fN) binSumw2 += oldBinw2[oldbin+i];
      }

      cu2[bin] = binContent;
      er2[bin] = binError;
      en2[bin] = binCount;
      if (fBinSumw2.fN) ew2[bin] = binSumw2;
      oldbin += imax;
   }
   // set bin statistics for underflow bin
   binContent = 0;
   binCount   = 0;
   binError   = 0;
   binSumw2   = 0;
   for(i=0;i<startbin;i++)
   {
      binContent += oldBins[i];
      binCount   += oldCount[i];
      binError   += oldErrors[i];
      if (fBinSumw2.fN) binSumw2 += oldBinw2[i];
   }
   hnew->fArray[0] = binContent;
   hnew->fBinEntries[0] = binCount;
   hnew->fSumw2[0] = binError;
   if ( fBinSumw2.fN ) hnew->fBinSumw2[0] = binSumw2;

   // set bin statistics for overflow bin
   binContent = 0;
   binCount   = 0;
   binError   = 0;
   binSumw2   = 0;
   for(i=oldbin;i<=nbins+1;i++)
   {
      binContent += oldBins[i];
      binCount   += oldCount[i];
      binError   += oldErrors[i];
      if (fBinSumw2.fN) binSumw2 += oldBinw2[i];
   }
   hnew->fArray[newbins+1] = binContent;
   hnew->fBinEntries[newbins+1] = binCount;
   hnew->fSumw2[newbins+1] = binError;
   if ( fBinSumw2.fN ) hnew->fBinSumw2[newbins+1] = binSumw2;


   delete [] oldBins;
   delete [] oldCount;
   delete [] oldErrors;
   if (oldBinw2) delete [] oldBinw2;
   return hnew;
}

//______________________________________________________________________________
void TProfile::ExtendAxis(Double_t x, TAxis *axis)
{
// Profile histogram is resized along x axis 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 extended histogram.
// Takes into account errors (Sumw2) if any.
// The axis must be extendable before invoking this function.
// Ex: h->GetXaxis()->SetCanExtend(kTRUE)

   TProfile*  hold = TProfileHelper::ExtendAxis(this, x, axis);
   if ( hold ) {
      fTsumwy  = hold->fTsumwy;
      fTsumwy2 = hold->fTsumwy2;

      delete hold;
   }
}

//______________________________________________________________________________
void TProfile::Reset(Option_t *option)
{
   //*-*-*-*-*-*-*-*-*-*Reset contents of a Profile histogram*-*-*-*-*-*-*-*-*
   //*-*                =====================================
   TH1D::Reset(option);
   fBinEntries.Reset();
   fBinSumw2.Reset();
   TString opt = option;
   opt.ToUpper();
   if (opt.Contains("ICE") && !opt.Contains("S")) return;
   fTsumwy  = 0;
   fTsumwy2 = 0;
}

//______________________________________________________________________________
void TProfile::SavePrimitive(std::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

   Bool_t nonEqiX = kFALSE;
   Int_t i;
   // Check if the profile has equidistant X bins or not.  If not, we
   // create an array holding the bins.
   if (GetXaxis()->GetXbins()->fN && GetXaxis()->GetXbins()->fArray) {
      nonEqiX = kTRUE;
      out << "   Double_t xAxis[" << GetXaxis()->GetXbins()->fN
      << "] = {";
      for (i = 0; i < GetXaxis()->GetXbins()->fN; i++) {
         if (i != 0) out << ", ";
         out << GetXaxis()->GetXbins()->fArray[i];
      }
      out << "}; " << std::endl;
   }

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

   //histogram pointer has by default teh histogram name.
   //however, in case histogram has no directory, it is safer to add a incremental suffix
   static Int_t hcounter = 0;
   TString histName = GetName();
   if (!fDirectory) {
      hcounter++;
      histName += "__";
      histName += hcounter;
   }
   const char *hname = histName.Data();

   out<<hname<<" = new "<<ClassName()<<"("<<quote<<GetName()<<quote<<","<<quote<<GetTitle()<<quote
   <<","<<GetXaxis()->GetNbins();
   if (nonEqiX)
      out << ", xAxis";
   else
      out << "," << GetXaxis()->GetXmin()
          << "," << GetXaxis()->GetXmax()
          <<","<<quote<<GetErrorOption()<<quote<<");"<<std::endl;

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

   TH1::SavePrimitiveHelp(out, hname, option);
}

//______________________________________________________________________________
void TProfile::Scale(Double_t c1, Option_t * option)
{
   // *-*-*-*-*Multiply this profile by a constant c1*-*-*-*-*-*-*-*-*
   // *-*      ======================================
   //
   //   this = c1*this
   //
   // This function uses the services of TProfile::Add
   //

   TProfileHelper::Scale(this, c1, option);
}

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

   TProfileHelper::SetBinEntries(this, bin, w);
}

//______________________________________________________________________________
void TProfile::SetBins(Int_t nx, Double_t xmin, Double_t xmax)
{
   //*-*-*-*-*-*-*-*-*Redefine  x axis parameters*-*-*-*-*-*-*-*-*-*-*-*
   //*-*              ===========================

   fXaxis.Set(nx,xmin,xmax);
   fNcells = nx+2;
   SetBinsLength(fNcells);
}

//______________________________________________________________________________
void TProfile::SetBins(Int_t nx, const Double_t *xbins)
{
   //*-*-*-*-*-*-*-*-*Redefine  x axis parameters*-*-*-*-*-*-*-*-*-*-*-*
   //*-*              ===========================

   fXaxis.Set(nx,xbins);
   fNcells = nx+2;
   SetBinsLength(fNcells);
}

//______________________________________________________________________________
void TProfile::SetBinsLength(Int_t n)
{
   // Set total number of bins including under/overflow
   // Reallocate bin contents array
   TH1D::SetBinsLength(n);
   TProfileHelper::BuildArray(this);
}

//______________________________________________________________________________
void TProfile::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 + 3*buffersize;
   fBuffer = new Double_t[fBufferSize];
   memset(fBuffer,0,sizeof(Double_t)*fBufferSize);
}

//______________________________________________________________________________
void TProfile::SetErrorOption(Option_t *option)
{
   //*-*-*-*-*-*-*-*-*-*Set option to compute profile errors*-*-*-*-*-*-*-*-*
   //*-*                =====================================
   //
   //    The computation of the bin errors is based on the parameter option:
   //    option:
   //     ' '  (Default) The bin errors are the standard error on the mean of the bin profiled values (Y),
   //                    i.e. the standard error of the bin contents.
   //                    Note that if TProfile::Approximate()  is called, an approximation is used when
   //                    the spread in Y is 0 and the number of bin entries  is > 0
   //
   //     's'            The bin errors are the standard deviations of the Y bin values
   //                    Note that if TProfile::Approximate()  is called, an approximation is used when
   //                    the spread in Y is 0 and the number of bin entries is > 0
   //
   //     'i'            Errors are as in default case (standard errors of the bin contents)
   //                    The only difference is for the case when the spread in Y is zero.
   //                    In this case for N > 0 the error is  1./SQRT(12.*N)
   //
   //     'g'            Errors are 1./SQRT(W)  for W not equal to 0 and 0 for W = 0.
   //                    W is the sum in the bin of the weights of the profile.
   //                    This option is for combining measurements y +/- dy,
   //                    and  the profile is filled with values y and weights w = 1/dy**2
   //
   //   See TProfile::BuildOptions for a detailed explanation of all options

   TProfileHelper::SetErrorOption(this, option);
}

//______________________________________________________________________________
void TProfile::Streamer(TBuffer &R__b)
{
   // Stream an object of class TProfile.

   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) {
         R__b.ReadClassBuffer(TProfile::Class(), this, R__v, R__s, R__c);
         return;
      }
      //====process old versions before automatic schema evolution
      TH1D::Streamer(R__b);
      fBinEntries.Streamer(R__b);
      Int_t errorMode;
      R__b >> errorMode;
      fErrorMode = (EErrorType)errorMode;
      if (R__v < 2) {
         Float_t ymin,ymax;
         R__b >> ymin; fYmin = ymin;
         R__b >> ymax; fYmax = ymax;
      } else {
         R__b >> fYmin;
         R__b >> fYmax;
      }
      R__b.CheckByteCount(R__s, R__c, TProfile::IsA());
      //====end of old versions

   } else {
      R__b.WriteClassBuffer(TProfile::Class(),this);
   }
}
//______________________________________________________________________________
void TProfile::Sumw2(Bool_t flag)
{
   // Create/delete structure to store sum of squares of weights per bin  *-*-*-*-*-*-*-*
   //   This is needed to compute  the correct statistical quantities
   //    of a profile filled with weights
   //
   //
   //  This function is automatically called when the histogram is created
   //  if the static function TH1::SetDefaultSumw2 has been called before.
   // If flag is false the structure is deleted

   TProfileHelper::Sumw2(this, flag);
}