// @(#)root/hist:$Name:  $:$Id: TH1.cxx,v 1.25 2000/09/13 16:49:53 brun Exp $
// Author: Rene Brun   26/12/94

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

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include <fstream.h>
#include <iostream.h>
#include <float.h>

#include "TROOT.h"
#include "TH1.h"
#include "TH2.h"
#include "TF2.h"
#include "TF3.h"
#include "TVirtualPad.h"
#include "Foption.h"
#include "TMath.h"
#include "TRandom.h"
#include "TVirtualFitter.h"
#include "TProfile.h"
#include "TStyle.h"

//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*The H I S T O G R A M   Classes*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ===============================
//*-*
//*-*  ROOT supports the following histogram types:
//*-*
//*-*   1-D histograms:
//*-*      TH1C : histograms with one byte per channel.   Maximum bin content = 255
//*-*      TH1S : histograms with one short per channel.  Maximum bin content = 65535
//*-*      TH1F : histograms with one float per channel.  Maximum precision 7 digits
//*-*      TH1D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-*   2-D histograms:
//*-*      TH2C : histograms with one byte per channel.   Maximum bin content = 255
//*-*      TH2S : histograms with one short per channel.  Maximum bin content = 65535
//*-*      TH2F : histograms with one float per channel.  Maximum precision 7 digits
//*-*      TH2D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-*   3-D histograms:
//*-*      TH3C : histograms with one byte per channel.   Maximum bin content = 255
//*-*      TH3S : histograms with one short per channel.  Maximum bin content = 65535
//*-*      TH3F : histograms with one float per channel.  Maximum precision 7 digits
//*-*      TH3D : histograms with one double per channel. Maximum precision 14 digits
//*-*
//*-*   Profile histograms: See classes  TProfile and TProfile2D
//*-*   Profile histograms are used to display the mean value of Y and its RMS
//*-*   for each bin in X. 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; 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.
///*-*
//*-*- All histogram classes are derived from the base class TH1
//*-*
//*-*                             TH1
//*-*                              ^
//*-*                              |
//*-*                              |
//*-*                              |
//*-*      -----------------------------------------------------
//*-*             |                |       |      |      |     |
//*-*             |                |      TH1C   TH1S   TH1F  TH1D
//*-*             |                |                           |
//*-*             |                |                           |
//*-*             |               TH2                      TProfile
//*-*             |                |
//*-*             |                |
//*-*             |                -----------------------------
//*-*             |                        |      |      |     |
//*-*             |                       TH2C   TH2S   TH2F  TH2D
//*-*             |                                            |
//*-*            TH3                                           |
//*-*             |                                        TProfile2D
//*-*             |
//*-*             ------------------------------
//*-*                     |      |      |      |
//*-*                    TH3C   TH3S   TH3F   TH3D
//*-*
//*-*   The TH*C classes also inherit from the array class TArrayC.
//*-*   The TH*S classes also inherit from the array class TArrayS.
//*-*   The TH*F classes also inherit from the array class TArrayF.
//*-*   The TH*D classes also inherit from the array class TArrayD.
//*-*
//*-*  Creating histograms
//*-*  ===================
//*-*  Histograms are created by invoking one of the constructors, eg
//*-*    TH1F *h1 = new TH1F("h1","h1 title",100,0,4.4);
//*-*    TH2F *h2 = new TH2F("h2","h2 title",40,0,4,30,-3,3);
//*-*  histograms may also be created by:
//*-*    - calling the Clone function, see below
//*-*    - making a projection from a 2-D or 3-D histogram, see below
//*-*    - reading an histogram from a file
//*-*  When an histogram is created, a reference to it is automatically added
//*-*  to the list of in-memory objects for the current file or directory.
//*-*  This default behaviour can be changed by:
//*-*    h->SetDirectory(0);         // for the current histogram h
//*-*    TH1::AddDirectory(kFALSE);  // sets a global switch disabling the reference
//*-*  When the histogram is deleted, the reference to it is removed from
//*-*  the list of objects in memory.
//*-*  When a file is closed, all histograms in memory associated with this file
//*-*  are automatically deleted.
//*-*
//*-*   Fix or variable bin size
//*-*   ========================
//*-*
//*-*  All histogram types support either fix or variable bin sizes.
//*-*  2-D histograms may have fix size bins along X and variable size bins
//*-*  along Y or vice-versa. The functions to fill, manipulate, draw or access
//*-*  histograms are identical in both cases.
//*-*  Each histogram always contains 3 objects TAxis: fXaxis, fYaxis and fZaxis
//*-*  To access the axis parameters, do:
//*-*     TAxis *xaxis = h->GetXaxis(); etc.
//*-*     Double_t binCenter = xaxis->GetBinCenter(bin), etc.
//*-*  See class TAxis for a description of all the access functions.
//*-*  The axis range is always stored internally in double precision.
//*-*
//*-*   Convention for numbering bins
//*-*   =============================
//*-*   For all histogram types: nbins, xlow, xup
//*-*     bin = 0;       underflow bin
//*-*     bin = 1;       first bin with low-edge xlow INCLUDED
//*-*     bin = nbins;   last bin with upper-edge xup EXCLUDED
//*-*     bin = nbins+1; overflow bin
//*-*   In case of 2-D or 3-D histograms, a "global bin" number is defined.
//*-*   For example, assuming a 3-D histogram with binx,biny,binz, the function
//*-*     Int_t bin = h->GetBin(binx,biny,binz);
//*-*   returns a global/linearized bin number. This
//*-*   to access the bin information independently of the dimension.
//*-*
//*-*  Filling histograms
//*-*  ==================
//*-*  An histogram is typically filled with statements like:
//*-*    h1->Fill(x);
//*-*    h1->Fill(x,w); //fill with weight
//*-*    h2->Fill(x,y)
//*-*    h2->Fill(x,y,w)
//*-*    h3->Fill(x,y,z)
//*-*    h3->Fill(x,y,z,w)
//*-*  The Fill functions compute the bin number corresponding to the given
//*-*  x,y or z argument and increment this bin by the given weight.
//*-*  The Fill functions return the bin number for 1-D histograms or global
//*-*  bin number for 2-D and 3-D histograms.
//*-*  If TH1::Sumw2 has been called before filling, the sum of squares of
//*-*  weights is also stored.
//*-*  One can also increment directly a bin number via TH1::AddBinContent
//*-*  or replace the existing content via TH1::SetBinContent.
//*-*  To access the bin content of a given bin, do:
//*-*    Double_t binContent = h->GetBinContent(bin);
//*-*
//*-*  By default, the bin number is computed using the current axis ranges.
//*-*  If the automatic binning option has been set via
//*-*         h->SetBit(TH1::kCanRebin);
//*-*  then, the Fill Function will automatically extend the axis range to
//*-*  accomodate the new value specified in the Fill argument. The method
//*-*  used is to double the bin size until the new value fits in the range,
//*-*  merging bins two by two. This automatic binning options is extensively
//*-*  used by the TTree::Draw function when histogramming Tree variables
//*-*  with an unknown range.
//*-*  This automatic binning option is supported for 1-d, 2-D and 3-D histograms.
//*-*
//*-*  During filling, some statistics parameters are incremented to compute
//*-*  the mean value and Root Mean Square with the maximum precision.
//*-*
//*-*  In case of histograms of type TH1C, TH1S, TH2C, TH2S, TH3C, TH3S
//*-*  a check is made that the bin contents do not exceed the maximum positive
//*-*  capacity (127 or 65535). Histograms of all types may have positive
//*-*  or/and negative bin contents.
//*-*
//*-*  Rebinning
//*-*  =========
//*-*  At any time, an histogram can be rebinned via TH1::Rebin. This function
//*-*  returns a new histogram with the rebinned contents.
//*-*  If bin errors were stored, they are recomputed during the rebinning.
//*-*
//*-*  Associated errors
//*-*  =================
//*-*  By default, for each bin, the sum of weights is computed at fill time.
//*-*  One can also call TH1::Sumw2 to force the storage and computation
//*-*  of the sum of the square of weights per bin.
//*-*  If Sumw2 has been called, the error per bin is computed as the
//*-*  sqrt(sum of squares of weights), otherwise the error is set equal
//*-*  to the sqrt(bin content).
//*-*  To return the error for a given bin number, do:
//*-*     Double_t error = h->GetBinError(bin);
//*-*
//*-*  Associated functions
//*-*  ====================
//*-*  One or more object (typically a TF1*) can be added to the list
//*-*  of functions (fFunctions) associated to each histogram.
//*-*  When TF1::Fit is invoked, the fitted function is added to this list.
//*-*  Given an histogram h, one can retrieve an associated function
//*-*  with:  TF1 *myfunc = h->GetFunction("myfunc");
//*-*
//*-*  Operations on histograms
//*-*  ========================
//*-*
//*-*  Many types of operations are supported on histograms or between histograms
//*-*  - Addition of an histogram to the current histogram
//*-*  - Additions of two histograms with coefficients and storage into the current
//*-*    histogram
//*-*  - Multiplications and Divisions are supported in the same way as additions.
//*-*  - The Add, Divide and Multiply functions also exist to add,divide or multiply
//*-*    an histogram by a function.
//*-*  If an histogram has associated error bars (TH1::Sumw2 has been called),
//*-*  the resulting error bars are also computed assuming independent histograms.
//*-*  In case of divisions, Binomial errors are also supported.
//*-*
//*-*
//*-*  Fitting histograms
//*-*  ==================
//*-*
//*-*  Histograms (1-D,2-D,3-D and Profiles) can be fitted with a user
//*-*  specified function via TH1::Fit. When an histogram is fitted, the
//*-*  resulting function with its parameters is added to the list of functions
//*-*  of this histogram. If the histogram is made persistent, the list of
//*-*  associated functions is also persistent. Given a pointer (see above)
//*-*  to an associated function myfunc, one can retrieve the function/fit
//*-*  parameters with calls such as:
//*-*    Double_t chi2 = myfunc->GetChisquare();
//*-*    Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
//*-*    Double_t err0 = myfunc->GetParError(0);  //error on first parameter
//*-*
//*-*
//*-*  Projections of histograms
//*-*  ========================
//*-*  One can:
//*-*   - make a 1-D projection of a 2-D histogram or Profile
//*-*     see functions TH2::ProjectionX,Y, TH2::ProfileX,Y, TProfile::ProjectionX
//*-*   - make a 1-D, 2-D or profile out of a 3-D histogram
//*-*     see functions TH3::ProjectionZ, TH3::Project3D.
//*-*
//*-*  One can fit these projections via:
//*-*   TH2::FitSlicesX,Y, TH3::FitSlicesZ.
//*-*
//*-*  Random Numbers and histograms
//*-*  =============================
//*-*  TH1::FillRandom can be used to randomly fill an histogram using
//*-*                 the contents of an existing TF1 function or another
//*-*                 TH1 histogram (for all dimensions).
//*-*  For example the following two statements create and fill an histogram
//*-*  10000 times with a default gaussian distribution of mean 0 and sigma 1:
//*-*    TH1F h1("h1","histo from a gaussian",100,-3,3);
//*-*    h1.FillRandom("gaus",10000);
//*-*  TH1::GetRandom can be used to return a random number distributed
//*-*                 according the contents of an histogram.
//*-*
//*-*  Making a copy of an histogram
//*-*  =============================
//*-*  Like for any other Root object derived from TObject, one can use
//*-*  the Clone function. This makes an identical copy of the original histogram
//*-*  including all associated errors and functions, eg:
//*-*    TH1F *hnew = (TH1F*)h->Clone();
//*-*    hnew->SetName("hnew"); //recommended, otherwise you get 2 histograms
//*-*                           //with the same name
//*-*
//*-*  Normalizing histograms
//*-*  ======================
//*-*  One can scale an histogram such that the bins integral is equal to
//*-*  to the normalization parameter via TH1::Scale(Double_t norm);
//*-*
//*-*  Drawing histograms
//*-*  ==================
//*-*  Histograms are drawn via the THistPainter class. Each histogram has
//*-*  a pointer to its own painter (to be usable in a multithreaded program).
//*-*  Many drawing options are supported.
//*-*  See THistPainter::Paint for more details.
//*-*  The same histogram can be drawn with different options in different pads.
//*-*  When an histogram drawn in a pad is deleted, the histogram is
//*-*  automatically removed from the pad or pads where it was drawn.
//*-*  If an histogram is drawn in a pad, then filled again, the new status
//*-*  of the histogram will be automatically shown in the pad next time
//*-*  the pad is updated. One does not need to redraw the histogram.
//*-*  To draw the current version of an histogram in a pad, one can use
//*-*     h->DrawCopy();
//*-*  This makes a clone (see Clone below) of the histogram. Once the clone
//*-*  is drawn, the original histogram may be modified or deleted without
//*-*  affecting the aspect of the clone.
//*-*
//*-*  One can use TH1::SetMaximum and TH1::SetMinimum to force a particular
//*-*  value for the maximum or the minimum scale on the plot.
//*-*
//*-*  TH1::UseCurrentStyle can be used to change all histogram graphics
//*-*  attributes to correspond to the current selected style.
//*-*  This function must be called for each histogram.
//*-*  In case one reads and draws many histograms from a file, one can force
//*-*  the histograms to inherit automatically the current graphics style
//*-*  by calling before gROOT->ForceStyle();
//*-*
//*-*
//*-*  Setting Drawing histogram contour levels (2-D hists only)
//*-*  =========================================================
//*-*  By default contours are automatically generated at equidistant
//*-*  intervals. A default value of 20 levels is used. This can be modified
//*-*  via TH1::SetContour or TH1::SetContourLevel.
//*-*  the contours level info is used by the drawing options "cont", "surf",
//*-*  and "lego".
//*-*
//*-*  Setting histogram graphics attributes
//*-*  =====================================
//*-*  The histogram classes inherit from the attribute classes:
//*-*    TAttLine, TAttFill, TAttMarker and TAttText.
//*-*  See the member functions of these classes for the list of options.
//*-*
//*-*  Giving titles to the X, Y and Z axis
//*-*  =================================
//*-*    h->GetXaxis()->SetTitle("X axis title");
//*-*    h->GetYaxis()->SetTitle("Y axis title");
//*-*  The histogram title and the axis titles can be any TLatex string.
//*-*  The titles are part of the persistent histogram.
//*-*
//*-*  Saving/Reading histograms to/from a Root file
//*-*  ================================
//*-*  The following statements create a Root file and store an histogram
//*-*  on the file. Because TH1 derives from TNamed, the key identifier on
//*-*  the file is the histogram name:
//*-*     TFile f("histos.root","new");
//*-*     TH1F h1("hgaus","histo from a gaussian",100,-3,3);
//*-*     h1.FillRandom("gaus",10000);
//*-*     h1->Write();
//*-*  To Read this histogram in another Root session, do:
//*-*     TFile f("histos.root");
//*-*     TH1F *h = (TH1F*)f.Get("hgaus");
//*-*  One can save all histograms in memory to the file by:
//*-*  file->Write();
//*-*
//*-*  Miscelaneous operations
//*-*  =======================
//*-*
//*-*   TH1::KolmogorovTest: Statistical test of compatibility in shape between
//*-*                        two histograms.
//*-*   TH1::Smooth smooths the bin contents of a 1-d histogram
//*-*   TH1::Integral returns the integral of bin contents in a given bin range
//*-*   TH1::GetMean(int axis) returns the mean value along axis
//*-*   TH1::GetRMS(int axis)  returns the Root Mean Square along axis
//*-*   TH1::GetEntries returns the number of entries
//*-*   TH1::Reset() resets the bin contents and errors of an histogram.
//*-*
//
/*

*/
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

Bool_t TH1::fgAddDirectory = kTRUE;

Foption_t Foption;

TAxis *xaxis=0;
TAxis *yaxis=0;
TAxis *zaxis=0;

TF1 *gF1=0;

TVirtualFitter *hFitter=0;

static Int_t xfirst,xlast,yfirst,ylast,zfirst,zlast;
static Axis_t xmin, xmax, ymin, ymax, zmin, zmax, binwidx, binwidy, binwidz;

extern void H1InitGaus();
extern void H1InitExpo();
extern void H1InitPolynom();
extern void H1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void H1FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a);
extern void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail);
extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);


ClassImp(TH1)

//______________________________________________________________________________
TH1::TH1(): TNamed(), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*-*-*-*-*Histogram default constructor*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =============================
   fDirectory = 0;
   fFunctions = new TList;
   fIntegral  = 0;
   fPainter   = 0;
   fXaxis.SetName("xaxis");
   fYaxis.SetName("yaxis");
   fZaxis.SetName("zaxis");
   fXaxis.SetParent(this);
   fYaxis.SetParent(this);
   fZaxis.SetParent(this);
}

//______________________________________________________________________________
TH1::~TH1()
{
//*-*-*-*-*-*-*-*-*-*-*Histogram default destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ============================

   if (!TestBit(kNotDeleted)) return;
   if (fIntegral) {delete [] fIntegral; fIntegral = 0;}
   if (fFunctions) { fFunctions->Delete(); delete fFunctions; }
   if (fDirectory) {
      if (!fDirectory->TestBit(TDirectory::kCloseDirectory))
         fDirectory->GetList()->Remove(this);
   }
   delete fPainter;
   fDirectory = 0;
   fFunctions = 0;
}

//______________________________________________________________________________
TH1::TH1(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*-*-*Normal constructor for fix bin size histograms*-*-*-*-*-*-*
//*-*              ==============================================
//
//     Creates the main histogram structure:
//        name   : name of histogram (avoid blanks)
//        title  : histogram title
//        nbins  : number of bins
//        xlow   : low edge of first bin
//        xup    : upper edge of last bin (not included in last bin)
//
//      When an histogram is created, it is automatically added to the list
//      of special objects in the current directory.
//      To find the pointer to this histogram in the current directory
//      by its name, do:
//      TH1F *h1 = (TH1F*)gDirectory->FindObject(name);
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   Build();
   if (nbins <= 0) nbins = 1;
   fXaxis.Set(nbins,xlow,xup);
   fNcells = fXaxis.GetNbins()+2;
}

//______________________________________________________________________________
TH1::TH1(const char *name,const char *title,Int_t nbins,Float_t *xbins)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*Normal constructor for variable bin size histograms*-*-*-*-*-*-*
//*-*          ===================================================
//
//  Creates the main histogram structure:
//     name   : name of histogram (avoid blanks)
//     title  : histogram title
//     nbins  : number of bins
//     xbins  : array of low-edges for each bin
//              This is an array of size nbins+1
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Build();
   if (nbins <= 0) nbins = 1;
   if (xbins) fXaxis.Set(nbins,xbins);
   else       fXaxis.Set(nbins,0,1);
   fNcells = fXaxis.GetNbins()+2;
}

//______________________________________________________________________________
TH1::TH1(const char *name,const char *title,Int_t nbins,Double_t *xbins)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
//*-*-*-*-*-*-*Normal constructor for variable bin size histograms*-*-*-*-*-*-*
//*-*          ===================================================
//
//  Creates the main histogram structure:
//     name   : name of histogram (avoid blanks)
//     title  : histogram title
//     nbins  : number of bins
//     xbins  : array of low-edges for each bin
//              This is an array of size nbins+1
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Build();
   if (nbins <= 0) nbins = 1;
   if (xbins) fXaxis.Set(nbins,xbins);
   else       fXaxis.Set(nbins,0,1);
   fNcells = fXaxis.GetNbins()+2;
}

//______________________________________________________________________________
void TH1::Browse(TBrowser *)
{
    Draw();
    gPad->Update();
}


//______________________________________________________________________________
void TH1::Build()
{
//*-*-*-*-*-*-*-*-*-*Creates histogram basic data structure*-*-*-*-*-*-*-*-*-*
//*-*                ======================================

   fPainter       = 0;
   fIntegral      = 0;
   fEntries       = 0;
   fNormFactor    = 0;
   fTsumw         = fTsumw2=fTsumwx=fTsumwx2=0;
   fMaximum       = -1111;
   fMinimum       = -1111;
   fXaxis.SetName("xaxis");
   fYaxis.SetName("yaxis");
   fZaxis.SetName("zaxis");
   fYaxis.Set(1,0.,1.);
   fZaxis.Set(1,0.,1.);
   fXaxis.SetParent(this);
   fYaxis.SetParent(this);
   fZaxis.SetParent(this);

   UseCurrentStyle();

   if (fgAddDirectory && gDirectory) {
      TH1 *hold = (TH1*)gDirectory->GetList()->FindObject(GetName());
      if (hold) {
         Warning("Build","Replacing existing histogram: %s",GetName());
         gDirectory->GetList()->Remove(hold);
       //  delete hold;
      }
      gDirectory->Append(this);
   }
   fDirectory = gDirectory;
   fFunctions = new TList;
}

//______________________________________________________________________________
void TH1::Add(TF1 *f1, Double_t c1)
{
   // Performs the operation: this = this + c1*f1
   // if errors are defined (see TH1::Sumw2), errors are also recalculated.

   if (!f1) {
      Error("Add","Attempt to add a non-existing function");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-*- Add statistics
   Stat_t s1[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = 0;}
   PutStats(s1);

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu;
   Double_t xx[3];
   Double_t *params = 0;
   f1->InitArgs(xx,params);
   for (binz=0;binz<=nbinsz+1;binz++) {
      xx[2] = fZaxis.GetBinCenter(binz);
      for (biny=0;biny<=nbinsy+1;biny++) {
         xx[1] = fYaxis.GetBinCenter(biny);
         for (binx=0;binx<=nbinsx+1;binx++) {
            xx[0] = fXaxis.GetBinCenter(binx);
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            cu  = c1*f1->EvalPar(xx);
            Double_t error1 = GetBinError(bin);
            AddBinContent(bin,cu);
            if (fSumw2.fN) {
               fSumw2.fArray[bin] = c1*c1*error1*error1;
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Add(TH1 *h1, Double_t c1)
{
   // Performs the operation: this = this + c1*h1
   // if errors are defined (see TH1::Sumw2), errors are also recalculated.
   // Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
   // if not already set.

   if (!h1) {
      Error("Add","Attempt to add a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
      Error("Add","Attempt to add histograms with different number of bins");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Add","Attempt to add histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 has Sumw2 set
   if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();

//*-*- Add statistics
   Double_t ac1 = TMath::Abs(c1);
   fEntries += ac1*h1->GetEntries();
   Stat_t s1[10], s2[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = s2[i] = 0;}
   GetStats(s1);
   h1->GetStats(s2);
   for (i=0;i<10;i++) s1[i] += ac1*s2[i];
   PutStats(s1);

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            cu  = c1*h1->GetBinContent(bin);
            AddBinContent(bin,cu);
            if (fSumw2.fN) {
               Double_t error1 = h1->GetBinError(bin);
               fSumw2.fArray[bin] += c1*c1*error1*error1;
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Add(TH1 *h1, TH1 *h2, Double_t c1, Double_t c2)
{
//*-*-*-*-*Replace contents of this histogram by the addition of h1 and h2*-*-*
//*-*      ===============================================================
//
//   this = c1*h1 + c2*h2
//   if errors are defined (see TH1::Sumw2), errors are also recalculated
//   Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
//   if not already set.
//

   if (!h1 || !h2) {
      Error("Add","Attempt to add a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
    || nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
      Error("Add","Attempt to add histograms with different number of bins");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Add","Attempt to add histograms with different axis limits");
   }
   if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
       Warning("Add","Attempt to add histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 or h2 have Sumw2 set
   if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();

//*-*- Add statistics
   Double_t ac1 = TMath::Abs(c1);
   Double_t ac2 = TMath::Abs(c2);
   fEntries = ac1*h1->GetEntries() + ac2*h2->GetEntries();
   Stat_t s1[10], s2[10], s3[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = s2[i] = s3[i] = 0;}
   GetStats(s1);
   h1->GetStats(s2);
   for (i=0;i<10;i++) s3[i] = ac1*s1[i] +ac2*s2[i];
   PutStats(s3);

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            cu  = c1*h1->GetBinContent(bin)+ c2*h2->GetBinContent(bin);
            SetBinContent(bin,cu);
            if (fSumw2.fN) {
               Double_t error1 = h1->GetBinError(bin);
               Double_t error2 = h2->GetBinError(bin);
               fSumw2.fArray[bin] = c1*c1*error1*error1 + c2*c2*error2*error2;
            }
         }
      }
   }
}


//______________________________________________________________________________
void TH1::AddBinContent(Int_t)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                ==========================
   AbstractMethod("AddBinContent");
}

//______________________________________________________________________________
void TH1::AddBinContent(Int_t, Stat_t)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by a weight w*-*-*-*-*-*-*-*-*-*-*
//*-*                ===================================
   AbstractMethod("AddBinContent");
}

//______________________________________________________________________________
void TH1::AddDirectory(Bool_t add)
{
// Sets the flag controlling the automatic add of histograms in memory
//
// By default (fAddDirectory = kTRUE), histograms are automatically added
// to the list of objects in memory.
// Note that one histogram can be removed from its support directory
// by calling h->SetDirectory(0) or h->SetDirectory(dir) to add it
// to the list of objects in the directory dir.
//
//  NOTE that this is a static function. To call it, use;
//     TH1::AddDirectory

   fgAddDirectory = add;
}

//______________________________________________________________________________
Double_t TH1::ComputeIntegral()
{
//  Compute integral (cumulative sum of bins)
//  The result stored in fIntegral is used by the GetRandom functions.
//  This function is automatically called by GetRandom when the fIntegral
//  array does not exist or when the number of entries in the histogram
//  has changed since the previous call to GetRandom.
//  The resulting integral is normalized to 1

   Int_t bin, binx, biny, binz, ibin;

// delete previously computed integral (if any)
   if (fIntegral) delete [] fIntegral;

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

   fIntegral = new Double_t[nbins+2];
   ibin = 0;
   fIntegral[ibin] = 0;
   Double_t dx,dy,dz;
   for (binz=1;binz<=nbinsz;binz++) {
      dz = fZaxis.GetBinWidth(binz);
      for (biny=1;biny<=nbinsy;biny++) {
         dy = fYaxis.GetBinWidth(biny);
         for (binx=1;binx<=nbinsx;binx++) {
            dx = fXaxis.GetBinWidth(binx);
            ibin++;
            bin  = GetBin(binx, biny, binz);
            fIntegral[ibin] = fIntegral[ibin-1] + GetBinContent(bin)*dx*dy*dz;
         }
      }
   }

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

//______________________________________________________________________________
void TH1::Copy(TObject &obj)
{
//*-*-*-*-*-*-*Copy this histogram structure to newth1*-*-*-*-*-*-*-*-*-*-*-*
//*-*          =======================================
//
// Note that this function does not copy the list of associated functions.
// Use TObJect::Clone to make a full copy of an histogram.

   TNamed::Copy(obj);
   ((TH1&)obj).fDimension = fDimension;
   ((TH1&)obj).fNormFactor= fNormFactor;
   ((TH1&)obj).fEntries   = fEntries;
   ((TH1&)obj).fNcells    = fNcells;
   ((TH1&)obj).fBarOffset = fBarOffset;
   ((TH1&)obj).fBarWidth  = fBarWidth;
   ((TH1&)obj).fTsumw     = fTsumw;
   ((TH1&)obj).fTsumw2    = fTsumw2;
   ((TH1&)obj).fTsumwx    = fTsumwx;
   ((TH1&)obj).fTsumwx2   = fTsumwx2;
   ((TH1&)obj).fMaximum   = fMaximum;
   ((TH1&)obj).fMinimum   = fMinimum;
   ((TH1&)obj).fOption    = fOption;
   TAttLine::Copy(((TH1&)obj));
   TAttFill::Copy(((TH1&)obj));
   TAttMarker::Copy(((TH1&)obj));
   fXaxis.Copy(((TH1&)obj).fXaxis);
   fYaxis.Copy(((TH1&)obj).fYaxis);
   fZaxis.Copy(((TH1&)obj).fZaxis);
   ((TH1&)obj).fXaxis.SetParent(this);
   ((TH1&)obj).fYaxis.SetParent(this);
   ((TH1&)obj).fZaxis.SetParent(this);
   fContour.Copy(((TH1&)obj).fContour);
   fSumw2.Copy(((TH1&)obj).fSumw2);
//   fFunctions->Copy(((TH1&)obj).fFunctions);
   gDirectory->Append(&obj);
//   ((TH1&)obj).AppendDirectory();
   ((TH1&)obj).fDirectory = gDirectory;
}

//______________________________________________________________________________
Int_t TH1::DistancetoPrimitive(Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Compute distance from point px,py to a line*-*-*-*-*-*
//*-*                  ===========================================
//*-*  Compute the closest distance of approach from point px,py to elements
//*-*  of an histogram.
//*-*  The distance is computed in pixels units.
//*-*
//*-*  Algorithm:
//*-*  Currently, this simple model computes the distance from the mouse
//*-*  to the histogram contour only.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (!fPainter) return 9999;
   return fPainter->DistancetoPrimitive(px,py);
}

//______________________________________________________________________________
void TH1::Divide(TF1 *f1, Double_t c1)
{
   // Performs the operation: this = this/(c1*f1)
   // if errors are defined (see TH1::Sumw2), errors are also recalculated.

   if (!f1) {
      Error("Add","Attempt to divide by a non-existing function");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-*- Add statistics
   Stat_t s1[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = 0;}
   PutStats(s1);

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu,w;
   Double_t xx[3];
   Double_t *params = 0;
   f1->InitArgs(xx,params);
   for (binz=0;binz<=nbinsz+1;binz++) {
      xx[2] = fZaxis.GetBinCenter(binz);
      for (biny=0;biny<=nbinsy+1;biny++) {
         xx[1] = fYaxis.GetBinCenter(biny);
         for (binx=0;binx<=nbinsx+1;binx++) {
            xx[0] = fXaxis.GetBinCenter(binx);
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            Double_t error1 = GetBinError(bin);
            cu  = c1*f1->EvalPar(xx);
            if (cu) w = GetBinContent(bin)/cu;
            else    w = 0;
            SetBinContent(bin,w);
            if (fSumw2.fN) {
               fSumw2.fArray[bin] = c1*c1*error1*error1;
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Divide(TH1 *h1)
{
//*-*-*-*-*-*-*-*-*-*-*Divide this histogram by h1*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ===========================
//
//   this = this/h1
//   if errors are defined (see TH1::Sumw2), errors are also recalculated.
//   Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
//   if not already set.
//   The resulting errors are calculated assuming uncorrelated histograms.
//   See the other TH1::Divide that gives the possibility to optionaly
//   compute Binomial errors.

   if (!h1) {
      Error("Divide","Attempt to divide by a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
      Error("Divide","Attempt to divide histograms with different number of bins");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Divide","Attempt to divide histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 has Sumw2 set
   if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();

//*-*- Reset statistics
   fEntries = fTsumw = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t c0,c1,w;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = GetBin(binx,biny,binz);
            c0  = GetBinContent(bin);
            c1  = h1->GetBinContent(bin);
            if (c1) w = c0/c1;
            else    w = 0;
            SetBinContent(bin,w);
            fEntries++;
            if (fSumw2.fN) {
               Double_t e0 = GetBinError(bin);
               Double_t e1 = h1->GetBinError(bin);
               Double_t c12= c1*c1;
               if (!c1) { fSumw2.fArray[bin] = 0; continue;}
               fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
            }
         }
      }
   }
   Stat_t s[10];
   GetStats(s);
   PutStats(s);
}


//______________________________________________________________________________
void TH1::Divide(TH1 *h1, TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
//*-*-*-*-*Replace contents of this histogram by the division of h1 by h2*-*-*
//*-*      ==============================================================
//
//   this = c1*h1/(c2*h2)
//
//   if errors are defined (see TH1::Sumw2), errors are also recalculated
//   Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
//   if not already set.
//   The resulting errors are calculated assuming uncorrelated histograms.
//   However, if option ="B" is specified, Binomial errors are computed.
//

   TString opt = option;
   opt.ToLower();
   Bool_t binomial = kFALSE;
   if (opt.Contains("b")) binomial = kTRUE;
   if (!h1 || !h2) {
      Error("Divide","Attempt to divide by a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
    || nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
      Error("Divide","Attempt to divide histograms with different number of bins");
      return;
   }
   if (!c2) {
      Error("Divide","Coefficient of dividing histogram cannot be zero");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Divide","Attempt to divide histograms with different axis limits");
   }
   if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
       Warning("Divide","Attempt to divide histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 or h2 have Sumw2 set
   if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();

//*-*- Reset statistics
   fEntries = fTsumw = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t b1,b2,w,d1,d2;
   d1 = c1*c1;
   d2 = c2*c2;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            b1  = h1->GetBinContent(bin);
            b2  = h2->GetBinContent(bin);
            if (b2) w = c1*b1/(c2*b2);
            else    w = 0;
            SetBinContent(bin,w);
            fEntries++;
            if (fSumw2.fN) {
               Double_t e1 = h1->GetBinError(bin);
               Double_t e2 = h2->GetBinError(bin);
               Double_t b22= b2*b2*d2;
               if (!b2) { fSumw2.fArray[bin] = 0; continue;}
               if (binomial) {
                  fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));
               } else {
                  fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22);
               }
            }
         }
      }
   }
   Stat_t s[10];
   GetStats(s);
   PutStats(s);
}

//______________________________________________________________________________
void TH1::Draw(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this histogram with options*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ================================
//*-*
//*-*  Histograms are drawn via the THistPainter class. Each histogram has
//*-*  a pointer to its own painter (to be usable in a multithreaded program).
//*-*  The same histogram can be drawn with different options in different pads.
//*-*  When an histogram drawn in a pad is deleted, the histogram is
//*-*  automatically removed from the pad or pads where it was drawn.
//*-*  If an histogram is drawn in a pad, then filled again, the new status
//*-*  of the histogram will be automatically shown in the pad next time
//*-*  the pad is updated. One does not need to redraw the histogram.
//*-*  To draw the current version of an histogram in a pad, one can use
//*-*     h->DrawCopy();
//*-*  This makes a clone of the histogram. Once the clone is drawn, the original
//*-*  histogram may be modified or deleted without affecting the aspect of the
//*-*  clone.
//*-*  By default, TH1::Draw clears the current pad.
//*-*
//*-*  One can use TH1::SetMaximum and TH1::SetMinimum to force a particular
//*-*  value for the maximum or the minimum scale on the plot.
//*-*
//*-*  TH1::UseCurrentStyle can be used to change all histogram graphics
//*-*  attributes to correspond to the current selected style.
//*-*  This function must be called for each histogram.
//*-*  In case one reads and draws many histograms from a file, one can force
//*-*  the histograms to inherit automatically the current graphics style
//*-*  by calling before gROOT->ForceStyle();
//*-*
//*-*  See THistPainter::Paint for a description of all the drawing options
//*-*  =======================
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) {
      //the following statement is necessary in case one attempts to draw
      //a temporary histogram already in the current pad
      if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
      gPad->Clear();
   }
   AppendPad(opt.Data());
}

//______________________________________________________________________________
TH1 *TH1::DrawCopy(Option_t *)
{
//*-*-*-*-*-*-*Copy this histogram and Draw in the current pad*-*-*-*-*-*-*-*
//*-*          ===============================================
//*-*
//*-*  Once the histogram is drawn into the pad, any further modification
//*-*  using graphics input will be made on the copy of the histogram,
//*-*  and not to the original object.
//*-*
//*-*  See Draw for the list of options
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   AbstractMethod("DrawCopy");
   return 0;
}

//______________________________________________________________________________
void TH1::DrawPanel()
{
//*-*-*-*-*-*-*Display a panel with all histogram drawing options*-*-*-*-*-*
//*-*          ==================================================
//*-*
//*-*   See class TDrawPanelHist for example

   if (fPainter) fPainter->DrawPanel();
}

//______________________________________________________________________________
void TH1::Eval(TF1 *f1, Option_t *option)
{
//*-*-*-*-*Evaluate function f1 at the center of bins of this histogram-*-*-*-*
//*-*      ============================================================
//*-*
//*-*  If option "R" is specified, the function is evaluated only
//*-*  for the bins included in the function range.
//*-*  If option "A" is specified, the value of the function is added to the
//*-*  existing bin contents
//*-*  If option "S" is specified, the value of the function is used to
//*-*  generate an integer value, distributed according to the Poisson
//*-*  distribution, with f1 as the mean.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Double_t x[3];
   Int_t stat,add,bin,binx,biny,binz,nbinsx, nbinsy, nbinsz;
   if (!f1) return;
   Double_t fu;
   Double_t e=0;
   TString opt = option;
   opt.ToLower();
   if (opt.Contains("a")) add  = 1;
   else                   add  = 0;
   if (opt.Contains("s")) stat = 1;
   else                   stat = 0;
   nbinsx  = fXaxis.GetNbins();
   nbinsy  = fYaxis.GetNbins();
   nbinsz  = fZaxis.GetNbins();
   if (!add) Reset();

   for (binz=1;binz<=nbinsz;binz++) {
      x[2]  = fZaxis.GetBinCenter(binz);
      for (biny=1;biny<=nbinsy;biny++) {
         x[1]  = fYaxis.GetBinCenter(biny);
         for (binx=1;binx<=nbinsx;binx++) {
            bin = GetBin(binx,biny,binz);
            x[0]  = fXaxis.GetBinCenter(binx);
            fu = f1->Eval(x[0],x[1],x[2]);
            if (stat) fu = gRandom->Poisson(fu);
            if (fSumw2.fN) e = fSumw2.fArray[bin];
            AddBinContent(bin,fu);
            if (fSumw2.fN) fSumw2.fArray[bin] = e+ fu*fu;
         }
      }
   }
}

//______________________________________________________________________________
void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
//*-*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*
//*-*                  =========================================
//*-*  This member function is called when a histogram is clicked with the locator
//*-*
//*-*  If Left button clicked on the bin top value, then the content of this bin
//*-*  is modified according to the new position of the mouse when it is released.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (fPainter) fPainter->ExecuteEvent(event, px, py);
}

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

   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(x);
   AddBinContent(bin);
   if (fSumw2.fN) ++fSumw2.fArray[bin];
   if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   return bin;
}

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

   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(x);
   AddBinContent(bin, w);
   if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
   if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
   Stat_t z= (w > 0 ? w : -w);
   fTsumw   += z;
   fTsumw2  += z*z;
   fTsumwx  += z*x;
   fTsumwx2 += z*x*x;
   return bin;
}

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

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

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

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

   Double_t *integral = new Double_t[nbinsx+1];
   ibin = 0;
   integral[ibin] = 0;
   for (binx=1;binx<=nbinsx;binx++) {
      xv[0] = fXaxis.GetBinCenter(binx);
//      xv[0] = fXaxis.GetBinUpEdge(binx);
      ibin++;
      integral[ibin] = integral[ibin-1] + f1->Eval(xv[0])*fXaxis.GetBinWidth(binx);
   }

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

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

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

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

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

   Int_t loop;
   Axis_t x;
   for (loop=0;loop<ntimes;loop++) {
      x = h->GetRandom();
      Fill(x);
   }
}


//______________________________________________________________________________
Int_t TH1::FindBin(Axis_t x, Axis_t y, Axis_t z)
{
//*-*-*-*-*-*Return Global bin number corresponding to x,y,z*-*-*-*-*-*-*
//*-*        ===============================================
//*-*
//*-*   2-D and 3-D histograms are represented with a one dimensional
//*-*   structure.
//*-*   This has the advantage that all existing functions, such as
//*-*     GetBinContent, GetBinError, GetBinFunction work for all dimensions.
//*-*  See also TH1::GetBin
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (GetDimension() < 2) {
      return fXaxis.FindBin(x);
   }
   if (GetDimension() < 3) {
      Int_t nx   = fXaxis.GetNbins()+2;
      Int_t binx = fXaxis.FindBin(x);
      Int_t biny = fYaxis.FindBin(y);
      return  binx + nx*biny;
   }
   if (GetDimension() < 4) {
      Int_t nx   = fXaxis.GetNbins()+2;
      Int_t ny   = fYaxis.GetNbins()+2;
      Int_t binx = fXaxis.FindBin(x);
      Int_t biny = fYaxis.FindBin(y);
      Int_t binz = fZaxis.FindBin(z);
      return  binx + nx*(biny +ny*binz);
   }
   return -1;
}

//______________________________________________________________________________
void TH1::Fit(const char *fname ,Option_t *option ,Option_t *goption, Axis_t xxmin, Axis_t xxmax)
{
//*-*-*-*-*-*-*-*-*-*-*Fit histogram with function fname*-*-*-*-*-*-*-*-*-*-*
//*-*                  =================================
//*-*
//*-*   fname is the name of an already predefined function created by TF1 or TF2
//*-*   Predefined functions such as gaus, expo and poln are automatically
//*-*   created by ROOT.
//*-*
//*-*   The list of fit options is given in parameter option.
//*-*      option = "W"  Set all errors to 1
//*-*             = "I" Use integral of function in bin instead of value at bin center
//*-*             = "L" Use Loglikelihood method (default is chisquare method)
//*-*             = "U" Use a User specified fitting algorithm (via SetFCN)
//*-*             = "Q" Quiet mode (minimum printing)
//*-*             = "V" Verbose mode (default is between Q and V)
//*-*             = "E" Perform better Errors estimation using Minos technique
//*-*             = "M" More. Improve fit results
//*-*             = "R" Use the Range specified in the function range
//*-*             = "N" Do not store the graphics function, do not draw
//*-*             = "0" Do not plot the result of the fit. By default the fitted function
//*-*                   is drawn unless the option"N" above is specified.
//*-*             = "+" Add this new fitted function to the list of fitted functions
//*-*                   (by default, any previous function is deleted)
//*-*
//*-*   When the fit is drawn (by default), the parameter goption may be used
//*-*   to specify a list of graphics options. See TH1::Draw for a complete
//*-*   list of these options.
//*-*
//*-*   In order to use the Range option, one must first create a function
//*-*   with the expression to be fitted. For example, if your histogram
//*-*   has a defined range between -4 and 4 and you want to fit a gaussian
//*-*   only in the interval 1 to 3, you can do:
//*-*        TF1 *f1 = new TF1("f1","gaus",1,3);
//*-*        histo->Fit("f1","R");
//*-*
//*-*   You can specify boundary limits for some or all parameters via
//*-*        f1->SetParLimits(p_number, parmin, parmax);
//*-*   if parmin>=parmax, the parameter is fixed
//*-*   Note that you are not forced to fix the limits for all parameters.
//*-*   For example, if you fit a function with 6 parameters, you can do:
//*-*     func->SetParameters(0,3.1,1.e-6,0.1,-8,100);
//*-*     func->SetParLimits(4,-10,-4);
//*-*     func->SetParLimits(5, 1,1);
//*-*   With this setup, parameters 0->3 can vary freely
//*-*   Parameter 4 has boundaries [-10,-4] with initial value -8
//*-*   Parameter 5 is fixed to 100.
//*-*
//*-*   Note that option "I" gives better results but is slower.
//*-*
//*-*
//*-*   Changing the fitting function
//*-*   =============================
//*-*  By default the fitting function H1FitChisquare is used.
//*-*  To specify a User defined fitting function, specify option "U" and
//*-*  call the following functions:
//*-*    TVirtualFitter::Fitter(myhist)->SetFCN(MyFittingFunction)
//*-*  where MyFittingFunction is of type:
//*-*  extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   Int_t i, npar,nvpar,nparx;
   Double_t par, we, al, bl;
   Double_t eplus,eminus,eparab,globcc,amin,edm,errdef,werr;
   Double_t params[100], arglist[100];
   TF1 *fnew1;
   TF2 *fnew2;
   TF3 *fnew3;

   xfirst  = fXaxis.GetFirst();
   xlast   = fXaxis.GetLast();
   binwidx = fXaxis.GetBinWidth(xlast);
   xmin    = fXaxis.GetBinLowEdge(xfirst);
   xmax    = fXaxis.GetBinLowEdge(xlast) +binwidx;
   yfirst  = fYaxis.GetFirst();
   ylast   = fYaxis.GetLast();
   binwidy = fYaxis.GetBinWidth(ylast);
   ymin    = fYaxis.GetBinLowEdge(yfirst);
   ymax    = fYaxis.GetBinLowEdge(ylast) +binwidy;
   zfirst  = fZaxis.GetFirst();
   zlast   = fZaxis.GetLast();
   binwidz = fZaxis.GetBinWidth(zlast);
   zmin    = fZaxis.GetBinLowEdge(zfirst);
   zmax    = fZaxis.GetBinLowEdge(zlast) +binwidz;
   xaxis   = &fXaxis;
   yaxis   = &fYaxis;
   zaxis   = &fZaxis;

//*-*- Check if Minuit is initialized and create special functions
   hFitter = TVirtualFitter::Fitter(this);
   hFitter->Clear();

//*-*- Get pointer to the function by searching in the list of functions in ROOT
   gF1 = (TF1*)gROOT->GetFunction(fname);
   if (!gF1) { Error("Fit", "Unknown function: %s",fname); return; }
   npar = gF1->GetNpar();
   if (npar <=0) { Error("Fit", "Illegal number of parameters = %d",npar); return; }

//*-*- Check that function has same dimension as histogram
   if (gF1->GetNdim() == 1 && GetDimension() > 1) {
      Error("Fit", "Function %s is not 2-D",fname); return; }
   if (gF1->GetNdim() == 2 && GetDimension() < 2) {
      Error("Fit", "Function %s is not 1-D",fname); return; }
   if (xxmin != xxmax) gF1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);

//*-*- Decode list of options into Foption
   if (!FitOptionsMake(option)) return;
   if (xxmin != xxmax) {
      gF1->SetRange(xxmin,ymin,zmin,xxmax,ymax,zmax);
      Foption.Range = 1;
   }
//*-*- Is a Fit range specified?
   Double_t fxmin, fymin, fzmin, fxmax, fymax, fzmax;
   if (Foption.Range) {
      gF1->GetRange(fxmin, fymin, fzmin, fxmax, fymax, fzmax);
      if (fxmin > xmin) xmin = fxmin;
      if (fymin > ymin) ymin = fymin;
      if (fzmin > zmin) zmin = fzmin;
      if (fxmax < xmax) xmax = fxmax;
      if (fymax < ymax) ymax = fymax;
      if (fzmax < zmax) zmax = fzmax;
      xfirst = fXaxis.FindFixBin(xmin); if (xfirst < 1) xfirst = 1;
      xlast  = fXaxis.FindFixBin(xmax); if (xlast > fXaxis.GetLast()) xlast = fXaxis.GetLast();
      yfirst = fYaxis.FindFixBin(ymin); if (yfirst < 1) yfirst = 1;
      ylast  = fYaxis.FindFixBin(ymax); if (ylast > fYaxis.GetLast()) ylast = fYaxis.GetLast();
      zfirst = fZaxis.FindFixBin(zmin); if (zfirst < 1) zfirst = 1;
      zlast  = fZaxis.FindFixBin(zmax); if (zlast > fZaxis.GetLast()) zlast = fZaxis.GetLast();
   } else {
      gF1->SetRange(xmin,ymin,zmin,xmax,ymax,zmax);
   }

//*-*- If case of a predefined function, then compute initial values of parameters
   Int_t special = gF1->GetNumber();
   if      (special == 100)      H1InitGaus();
   else if (special == 400)      H1InitGaus();
   else if (special == 200)      H1InitExpo();
   else if (special == 299+npar) H1InitPolynom();

//*-*- Some initialisations
   if (!Foption.Verbose) {
      arglist[0] = -1;
      hFitter->ExecuteCommand("SET PRINT", arglist,1);
      arglist[0] = 0;
      hFitter->ExecuteCommand("SET NOW",   arglist,0);
   }

//*-*- Set error criterion for chisquare or likelihood methods
//*-*-  MINUIT ERRDEF should not be set to 0.5 in case of loglikelihood fit.
//*-*-  because the FCN is already multiplied by 2 in H1FitLikelihood
//*-*-  if Hoption.User is specified, assume that the user has already set
//*-*-  his minimization function via SetFCN.
   arglist[0] = 1;
   if (Foption.Like) {
      hFitter->SetFCN(H1FitLikelihood);
   } else {
      if (!Foption.User) hFitter->SetFCN(H1FitChisquare);
   }
   hFitter->ExecuteCommand("SET ERR",arglist,1);

//*-*- Transfer names and initial values of parameters to Minuit
   Int_t nfixed = 0;
   for (i=0;i<npar;i++) {
      par = gF1->GetParameter(i);
      gF1->GetParLimits(i,al,bl);
      if (al*bl != 0 && al >= bl) {
         al = bl = 0;
         arglist[nfixed] = i+1;
         nfixed++;
      }
      we = 0.1*TMath::Abs(bl-al);
      if (we == 0) we = 0.3*TMath::Abs(par);
//      if (we == 0) we = 1000*TMath::Abs(par);
      if (we == 0) we = binwidx;
      hFitter->SetParameter(i,gF1->GetParName(i),par,we,al,bl);
   }
   if(nfixed > 0)hFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto

//*-*- Set Gradient
   if (Foption.Gradient) {
      if (Foption.Gradient == 1) arglist[0] = 1;
      else                       arglist[0] = 0;
      hFitter->ExecuteCommand("SET GRAD",arglist,1);
   }

//*-*- Reset Print level
   if (Foption.Verbose) {
      arglist[0] = 0; hFitter->ExecuteCommand("SET PRINT", arglist,1);
   }

//*-*- Compute sum of squares of errors in the bin range
   Double_t ey, sumw2=0;
   for (i=xfirst;i<=xlast;i++) {
      ey = GetBinError(i);
      sumw2 += ey*ey;
   }

//*-*- Perform minimization
   arglist[0] = TVirtualFitter::GetMaxIterations();
   arglist[1] = sumw2*TVirtualFitter::GetPrecision();
   hFitter->ExecuteCommand("MIGRAD",arglist,2);
   if (Foption.More) {
      hFitter->ExecuteCommand("IMPROVE",arglist,0);
   }
   if (Foption.Errors) {
      hFitter->ExecuteCommand("HESSE",arglist,0);
      hFitter->ExecuteCommand("MINOS",arglist,0);
   }

//*-*- Get return status
   char parName[50];
   for (i=0;i<npar;i++) {
      hFitter->GetParameter(i,parName, par,we,al,bl);
      if (Foption.Errors) werr = we;
      else {
         hFitter->GetErrors(i,eplus,eminus,eparab,globcc);
         if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);
         else                         werr = we;
      }
      params[i] = par;
      gF1->SetParameter(i,par);
      gF1->SetParError(i,werr);
   }
   hFitter->GetStats(amin,edm,errdef,nvpar,nparx);

//*-*- Print final values of parameters.
   if (!Foption.Quiet) {
      if (Foption.Errors) hFitter->PrintResults(4,amin);
      else                hFitter->PrintResults(3,amin);
   }

//*-*  If Log Likelihood, compute an equivalent chisquare
   if (Foption.Like) H1FitChisquare(npar, params, amin, params, 1);

   gF1->SetChisquare(amin);

//*-*- Store fitted function in histogram functions list and draw
   if (!Foption.Nostore) {
      if (!Foption.Plus) {
         TIter next(fFunctions, kIterBackward);
         TObject *obj;
         while ((obj = next())) {
            if (obj->InheritsFrom(TF1::Class())) {
               fFunctions->Remove(obj);
               delete obj;
            }
         }
      }
      if (GetDimension() < 2) {
         fnew1 = new TF1();
         gF1->Copy(*fnew1);
         fFunctions->Add(fnew1);
         fnew1->SetParent(this);
         fnew1->Save(xmin,xmax);
      } else if (GetDimension() < 3) {
         fnew2 = new TF2();
         gF1->Copy(*fnew2);
         fFunctions->Add(fnew2);
         fnew2->SetParent(this);
      } else {
         fnew3 = new TF3();
         gF1->Copy(*fnew3);
         fFunctions->Add(fnew3);
         fnew3->SetParent(this);
      }
      if (TestBit(kCanDelete)) return;
      if (!Foption.Nograph && GetDimension() < 3) Draw(goption);
   }
}

//______________________________________________________________________________
void TH1::FitPanel()
{
//*-*-*-*-*-*-*Display a panel with all histogram fit options*-*-*-*-*-*
//*-*          ==============================================
//*-*
//*-*   See class TFitPanel for example

   if (fPainter) fPainter->FitPanel();
}

//______________________________________________________________________________
char *TH1::GetObjectInfo(Int_t px, Int_t py)
{
//   Redefines TObject::GetObjectInfo.
//   Displays the histogram info (bin number, contents, integral up to bin
//   corresponding to cursor position px,py
//
   return fPainter->GetObjectInfo(px,py);
}

//______________________________________________________________________________
Int_t TH1::FitOptionsMake(Option_t *choptin)
{
//*-*-*-*-*-*-*-*-*Decode string choptin and fill Foption structure*-*-*-*-*-*
//*-*              ================================================

   Foption.Quiet   = 0;
   Foption.Verbose = 0;
   Foption.Bound   = 0;
   Foption.Like    = 0;
   Foption.User    = 0;
   Foption.W1      = 0;
   Foption.Errors  = 0;
   Foption.More    = 0;
   Foption.Range   = 0;
   Foption.Gradient= 0;
   Foption.Nograph = 0;
   Foption.Nostore = 0;
   Foption.Plus    = 0;
   Foption.Integral= 0;

   Int_t nch = strlen(choptin);
   if (!nch) return 1;

   char chopt[32];
   strcpy(chopt,choptin);

   for (Int_t i=0;i<nch;i++) chopt[i] = toupper(choptin[i]);

   if (strstr(chopt,"Q")) Foption.Quiet   = 1;
   if (strstr(chopt,"V")){Foption.Verbose = 1; Foption.Quiet = 0;}
   if (strstr(chopt,"L")) Foption.Like    = 1;
   if (strstr(chopt,"W")) Foption.W1      = 1;
   if (strstr(chopt,"E")) Foption.Errors  = 1;
   if (strstr(chopt,"M")) Foption.More    = 1;
   if (strstr(chopt,"R")) Foption.Range   = 1;
   if (strstr(chopt,"G")) Foption.Gradient= 1;
   if (strstr(chopt,"N")) Foption.Nostore = 1;
   if (strstr(chopt,"0")) Foption.Nograph = 1;
   if (strstr(chopt,"+")) Foption.Plus    = 1;
   if (strstr(chopt,"I")) Foption.Integral= 1;
   if (strstr(chopt,"U")){Foption.User    = 1; Foption.Like = 0;}
   return 1;
}

//______________________________________________________________________________
void H1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//*-*-*-*-*-*Minimization function for H1s using a Chisquare method*-*-*-*-*-*
//*-*        ======================================================

   Double_t cu,eu,fu,fsum;
   Double_t dersum[100], grad[100];
   Double_t x[3];
   Int_t bin,binx,biny,binz,k;
   Axis_t binlow, binup, binsize;

   Int_t npfits = 0;

   npar = gF1->GetNpar();
   if (flag == 2) for (k=0;k<npar;k++) dersum[k] = gin[k] = 0;

   TH1 *hfit = (TH1*)hFitter->GetObjectFit();
   gF1->InitArgs(x,u);
   f = 0;
   for (binz=zfirst;binz<=zlast;binz++) {
      x[2]  = zaxis->GetBinCenter(binz);
      for (biny=yfirst;biny<=ylast;biny++) {
         x[1]  = yaxis->GetBinCenter(biny);
         for (binx=xfirst;binx<=xlast;binx++) {
            bin = hfit->GetBin(binx,biny,binz);
            cu  = hfit->GetBinContent(bin);
            x[0]  = xaxis->GetBinCenter(binx);
            if (Foption.Integral) {
               binlow  = xaxis->GetBinLowEdge(binx);
               binsize = xaxis->GetBinWidth(binx);
               binup   = binlow + binsize;
               fu      = gF1->Integral(binlow,binup,u)/binsize;
            } else {
               fu = gF1->EvalPar(x,u);
            }
            if (Foption.W1) {
               if (cu == 0) continue;
               eu = 1;
            } else {
               eu  = hfit->GetBinError(bin);
               if (eu <= 0) continue;
            }
            if (flag == 2) {
               for (k=0;k<npar;k++) dersum[k] += 1; //should be the derivative
            }
            npfits++;
            if (flag == 2) {
               for (k=0;k<npar;k++) grad[k] += dersum[k]*(fu-cu)/eu; dersum[k] = 0;
            }
            fsum = (cu-fu)/eu;
            f += fsum*fsum;
         }
      }
   }
//printf("f=%g, npfits=%dn",f,npfits);
   gF1->SetNumberFitPoints(npfits);
}

//______________________________________________________________________________
void H1FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//*-*-*-*-*-*Minimization function for H1s using a Likelihood method*-*-*-*-*-*
//*-*        =======================================================
//     Basically, it forms the likelihood by determining the Poisson
//     probability that given a number of entries in a particualar bin,
//     the fit would predict it's value.  This is then done for each bin,
//     and the sum of the logs is taken as the likelihood.

   Double_t cu,fu,fobs,fsub;
   Double_t dersum[100];
   Double_t x[3];
   Int_t bin,binx,biny,binz,k,icu;
   Axis_t binlow, binup, binsize;

   Int_t npfits = 0;

   npar = gF1->GetNpar();
   if (flag == 2) for (k=0;k<npar;k++) dersum[k] = gin[k] = 0;

   TH1 *hfit = (TH1*)hFitter->GetObjectFit();
   gF1->InitArgs(x,u);
   f = 0;
   for (binz=zfirst;binz<=zlast;binz++) {
      x[2]  = zaxis->GetBinCenter(binz);
      for (biny=yfirst;biny<=ylast;biny++) {
         x[1]  = yaxis->GetBinCenter(biny);
         for (binx=xfirst;binx<=xlast;binx++) {
            bin = hfit->GetBin(binx,biny,binz);
            cu  = hfit->GetBinContent(bin);
            x[0]  = xaxis->GetBinCenter(binx);
            if (Foption.Integral) {
               binlow  = xaxis->GetBinLowEdge(binx);
               binsize = xaxis->GetBinWidth(binx);
               binup   = binlow + binsize;
               fu      = gF1->Integral(binlow,binup,u)/binsize;
            } else {
               fu = gF1->EvalPar(x,u);
            }
            npfits++;
            if (flag == 2) {
               for (k=0;k<npar;k++) {
                  dersum[k] += 1; //should be the derivative
                  //grad[k] += dersum[k]*(fu-cu)/eu; dersum[k] = 0;
               }
            }
            if (fu < 1.e-9) fu = 1.e-9;
            icu  = Int_t(cu);
            fsub = -fu +icu*TMath::Log(fu);
            fobs = hFitter->GetSumLog(icu);

            fsub -= fobs;
            f -= fsub;
         }
      }
   }
   f *= 2;
   gF1->SetNumberFitPoints(npfits);
}

//______________________________________________________________________________
void H1InitGaus()
{
//*-*-*-*-*-*Compute Initial values of parameters for a gaussian*-*-*-*-*-*-*
//*-*        ===================================================

   Double_t allcha, sumx, sumx2, x, val, rms, mean;
   Int_t bin;
   static Double_t sqrtpi = 2.506628;

//*-*- Compute mean value and RMS of the histogram in the given range
   TH1 *curHist = (TH1*)hFitter->GetObjectFit();
   Double_t valmax = curHist->GetBinContent(xfirst);
   allcha = sumx = sumx2 = 0;
   for (bin=xfirst;bin<=xlast;bin++) {
      x       = curHist->GetBinCenter(bin);
      val     = TMath::Abs(curHist->GetBinContent(bin));
      if (val > valmax) valmax = val;
      sumx   += val*x;
      sumx2  += val*x*x;
      allcha += val;
   }
   if (allcha == 0) return;
   mean = sumx/allcha;
   rms  = TMath::Sqrt(sumx2/allcha - mean*mean);
   if (rms == 0) rms = binwidx*(xlast-xfirst+1)/4;
   //if the distribution is really gaussian, the best approximation
   //is binwidx*allcha/(sqrtpi*rms)
   //However, in case of non-gaussian tails, this underestimates
   //the normalisation constant. In this case the maximum value
   //is a better approximation.
   //We take the average of both quantities
   Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*rms));

   //In case the mean value is outside the histo limits and
   //the RMS is bigger than the range, we take
   //  mean = center of bins
   //  rms  = half range
   Axis_t xmin = curHist->GetXaxis()->GetXmin();
   Axis_t xmax = curHist->GetXaxis()->GetXmax();
   if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
      mean = 0.5*(xmax+xmin);
      rms  = 0.5*(xmax-xmin);
   }
   gF1->SetParameter(0,constant);
   gF1->SetParameter(1,mean);
   gF1->SetParameter(2,rms);
   gF1->SetParLimits(2,0,10*rms);
}

//______________________________________________________________________________
void H1InitExpo()
{
//*-*-*-*-*-*Compute Initial values of parameters for an exponential*-*-*-*-*
//*-*        =======================================================

   Double_t constant, slope;
   Int_t ifail;
   Int_t nchanx = xlast - xfirst + 1;

   H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);

   gF1->SetParameter(0,constant);
   gF1->SetParameter(1,slope);

}

//______________________________________________________________________________
void H1InitPolynom()
{
//*-*-*-*-*-*Compute Initial values of parameters for a polynom*-*-*-*-*-*-*
//*-*        ===================================================

   Double_t fitpar[25];

   Int_t nchanx = xlast - xfirst + 1;
   Int_t npar   = gF1->GetNpar();

   if (nchanx <=1 || npar == 1) {
      TH1 *curHist = (TH1*)hFitter->GetObjectFit();
      fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
   } else {
      H1LeastSquareFit( nchanx, npar, fitpar);
   }
   for (Int_t i=0;i<npar;i++) gF1->SetParameter(i, fitpar[i]);
}

//______________________________________________________________________________
void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a)
{
//*-*-*-*-*-*-*-*Least squares lpolynomial fitting without weights*-*-*-*-*-*-*
//*-*            =================================================
//*-*
//*-*  n   number of points to fit
//*-*  m   number of parameters
//*-*  a   array of parameters
//*-*
//*-*   based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
//*-*   (E.Keil.  revised by B.Schorr, 23.10.1981.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    const Double_t zero = 0.;
    const Double_t one = 1.;
    const Int_t idim = 20;

    Double_t  b[400]	/* was [20][20] */;
    Int_t i, k, l, ifail;
    Double_t power;
    Double_t da[20], xk, yk;

    if (m <= 2) {
       H1LeastSquareLinearFit(n, a[0], a[1], ifail);
       return;
    }
    if (m > idim || m > n) return;
    b[0]  = Double_t(n);
    da[0] = zero;
    for (l = 2; l <= m; ++l) {
	b[l-1]           = zero;
	b[m + l*20 - 21] = zero;
	da[l-1]          = zero;
    }
    TH1 *curHist = (TH1*)hFitter->GetObjectFit();
    for (k = xfirst; k <= xlast; ++k) {
	xk     = curHist->GetBinCenter(k);
	yk     = curHist->GetBinContent(k);
	power  = one;
	da[0] += yk;
	for (l = 2; l <= m; ++l) {
	    power   *= xk;
	    b[l-1]  += power;
	    da[l-1] += power*yk;
	}
	for (l = 2; l <= m; ++l) {
	    power            *= xk;
	    b[m + l*20 - 21] += power;
	}
    }
    for (i = 3; i <= m; ++i) {
	for (k = i; k <= m; ++k) {
	    b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
	}
    }
    H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);

    for (i=0; i<m; ++i) a[i] = da[i];

}

//______________________________________________________________________________
void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
{
//*-*-*-*-*-*-*-*-*-*Least square linear fit without weights*-*-*-*-*-*-*-*-*
//*-*                =======================================
//*-*
//*-*   extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
//*-*   (added to LSQ by B. Schorr, 15.02.1982.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

    Double_t xbar, ybar, x2bar;
    Int_t i, n;
    Double_t xybar;
    Double_t fn, xk, yk;
    Double_t det;

    n     = TMath::Abs(ndata);
    ifail = -2;
    xbar  = ybar = x2bar = xybar = 0;
    TH1 *curHist = (TH1*)hFitter->GetObjectFit();
    for (i = xfirst; i <= xlast; ++i) {
	xk = curHist->GetBinCenter(i);
	yk = curHist->GetBinContent(i);
	if (ndata < 0) {
	    if (yk <= 0) yk = 1e-9;
	    yk = TMath::Log(yk);
	}
	xbar  += xk;
	ybar  += yk;
	x2bar += xk*xk;
	xybar += xk*yk;
    }
    fn    = Double_t(n);
    det   = fn*x2bar - xbar*xbar;
    ifail = -1;
    if (det <= 0) {
       a0 = ybar/fn;
       a1 = 0;
       return;
    }
    ifail = 0;
    a0 = (x2bar*ybar - xbar*xybar) / det;
    a1 = (fn*xybar - xbar*ybar) / det;

}

//______________________________________________________________________________
void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
{
//*-*-*-*-*-*-*-*Extracted from CERN Program library routine DSEQN*-*-*-*-*-*
//*-*            =================================================
//*-*
//*-*        : Translated to C++ by Rene Brun
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    Int_t a_dim1, a_offset, b_dim1, b_offset;
    Int_t nmjp1, i, j, l;
    Int_t im1, jp1, nm1, nmi;
    Double_t s1, s21, s22;
    const Double_t one = 1.;

    /* Parameter adjustments */
    b_dim1 = idim;
    b_offset = b_dim1 + 1;
    b -= b_offset;
    a_dim1 = idim;
    a_offset = a_dim1 + 1;
    a -= a_offset;

    if (idim < n) return;

    ifail = 0;
    for (j = 1; j <= n; ++j) {
	if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
	a[j + j*a_dim1] = one / a[j + j*a_dim1];
	if (j == n) continue;
	jp1 = j + 1;
	for (l = jp1; l <= n; ++l) {
	    a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
	    s1 = -a[l + (j+1)*a_dim1];
	    for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
	    a[l + (j+1)*a_dim1] = -s1;
	}
    }
    if (k <= 0) return;

    for (l = 1; l <= k; ++l) {
	b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
    }
    if (n == 1) return;
    for (l = 1; l <= k; ++l) {
	for (i = 2; i <= n; ++i) {
	    im1 = i - 1;
	    s21 = -b[i + l*b_dim1];
	    for (j = 1; j <= im1; ++j) {
		s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
	    }
	    b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
	}
	nm1 = n - 1;
	for (i = 1; i <= nm1; ++i) {
	    nmi = n - i;
	    s22 = -b[nmi + l*b_dim1];
	    for (j = 1; j <= i; ++j) {
		nmjp1 = n - j + 1;
		s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
	    }
	    b[nmi + l*b_dim1] = -s22;
	}
    }
}


//______________________________________________________________________________
Int_t TH1::GetBin(Int_t binx, Int_t biny, Int_t binz)
{
//*-*-*-*-*-*Return Global bin number corresponding to binx,y,z*-*-*-*-*-*-*
//*-*        ==================================================
//*-*
//*-*   2-D and 3-D histograms are represented with a one dimensional
//*-*   structure.
//*-*   This has the advantage that all existing functions, such as
//*-*     GetBinContent, GetBinError, GetBinFunction work for all dimensions.
//*-*
//*-*  In case of a TH1x, returns binx directly.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Int_t nx, ny, nz;
   if (GetDimension() < 2) {
      nx  = fXaxis.GetNbins()+2;
      if (binx < 0)   binx = 0;
      if (binx >= nx) binx = nx-1;
      return binx;
   }
   if (GetDimension() < 3) {
      nx  = fXaxis.GetNbins()+2;
      if (binx < 0)   binx = 0;
      if (binx >= nx) binx = nx-1;
      ny  = fYaxis.GetNbins()+2;
      if (biny < 0)   biny = 0;
      if (biny >= ny) biny = ny-1;
      return  binx + nx*biny;
   }
   if (GetDimension() < 4) {
      nx  = fXaxis.GetNbins()+2;
      if (binx < 0)   binx = 0;
      if (binx >= nx) binx = nx-1;
      ny  = fYaxis.GetNbins()+2;
      if (biny < 0)   biny = 0;
      if (biny >= ny) biny = ny-1;
      nz  = fZaxis.GetNbins()+2;
      if (binz < 0)   binz = 0;
      if (binz >= nz) binz = nz-1;
      return  binx + nx*(biny +ny*binz);
   }
   return -1;
}

//______________________________________________________________________________
Axis_t TH1::GetRandom()
{
// return a random number distributed according the histogram bin contents.
// This function checks if the bins integral exists. If not, the integral
// is evaluated, normalized to one.
// The integral is automatically recomputed if the number of entries
// is not the same then when the integral was computed.
// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.

   if (fDimension > 1) {
      Error("GetRandom","Function only valid for 1-d histograms");
      return 0;
   }
   Int_t nbinsx = GetNbinsX();
   Double_t integral;
   if (fIntegral) {
      if (fIntegral[nbinsx+1] != fEntries) integral = ComputeIntegral();
   } else {
      integral = ComputeIntegral();
      if (integral == 0 || fIntegral == 0) return 0;
   }
   Double_t r1 = gRandom->Rndm();
   Int_t ibin = TMath::BinarySearch(nbinsx,&fIntegral[0],r1);
   return GetBinLowEdge(ibin+1)
      +GetBinWidth(ibin+1)*(fIntegral[ibin+1]-r1)/(fIntegral[ibin+1] - fIntegral[ibin]);
}

//______________________________________________________________________________
Stat_t TH1::GetBinContent(Int_t)
{
//*-*-*-*-*-*-*Return content of global bin number bin*-*-*-*-*-*-*-*-*-*-*-*
//*-*          =======================================
   AbstractMethod("GetBinContent");
   return 0;
}

//______________________________________________________________________________
void TH1::Multiply(TF1 *f1, Double_t c1)
{
   // Performs the operation: this = this*c1*f1
   // if errors are defined (see TH1::Sumw2), errors are also recalculated.

   if (!f1) {
      Error("Add","Attempt to multiply by a non-existing function");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-*- Add statistics
   Stat_t s1[10];
   Int_t i;
   for (i=0;i<10;i++) {s1[i] = 0;}
   PutStats(s1);

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu,w;
   Double_t xx[3];
   Double_t *params = 0;
   f1->InitArgs(xx,params);
   for (binz=0;binz<=nbinsz+1;binz++) {
      xx[2] = fZaxis.GetBinCenter(binz);
      for (biny=0;biny<=nbinsy+1;biny++) {
         xx[1] = fYaxis.GetBinCenter(biny);
         for (binx=0;binx<=nbinsx+1;binx++) {
            xx[0] = fXaxis.GetBinCenter(binx);
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            Double_t error1 = GetBinError(bin);
            cu  = f1->EvalPar(xx);
            w = GetBinContent(bin)*c1*cu;
            SetBinContent(bin,w);
            if (fSumw2.fN) {
               fSumw2.fArray[bin] = c1*c1*error1*error1;
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Multiply(TH1 *h1)
{
//*-*-*-*-*-*-*-*-*-*-*Multiply this histogram by h1*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  =============================
//
//   this = this*h1
//
//   If errors of this are available (TH1::Sumw2), errors are recalculated.
//   Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
//   if not already set.

   if (!h1) {
      Error("Multiply","Attempt to multiply by a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
      Error("Multiply","Attempt to multiply histograms with different number of bins");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Multiply","Attempt to multiply histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 has Sumw2 set
   if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();

//*-*- Reset statistics
   fEntries = fTsumw = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t c0,c1,w;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = GetBin(binx,biny,binz);
            c0  = GetBinContent(bin);
            c1  = h1->GetBinContent(bin);
            w   = c0*c1;
            SetBinContent(bin,w);
            fEntries++;
            if (fSumw2.fN) {
               Double_t e0 = GetBinError(bin);
               Double_t e1 = h1->GetBinError(bin);
               fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0);
            }
         }
      }
   }
   Stat_t s[10];
   GetStats(s);
   PutStats(s);
}


//______________________________________________________________________________
void TH1::Multiply(TH1 *h1, TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
//*-*-*-*-*Replace contents of this histogram by multiplication of h1 by h2*-*
//*-*      ================================================================
//
//   this = (c1*h1)*(c2*h2)
//
//   If errors of this are available (TH1::Sumw2), errors are recalculated.
//   Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
//   if not already set.
//

   TString opt = option;
   opt.ToLower();
//   Bool_t binomial = kFALSE;
//   if (opt.Contains("b")) binomial = kTRUE;
   if (!h1 || !h2) {
      Error("Multiply","Attempt to multiply by a non-existing histogram");
      return;
   }

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
//*-*- Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
    || nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
      Error("Multiply","Attempt to multiply histograms with different number of bins");
      return;
   }
//*-*- Issue a Warning if histogram limits are different
   if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Multiply","Attempt to multiply histograms with different axis limits");
   }
   if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
       GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
       GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
       GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
       GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
       GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
       Warning("Multiply","Attempt to multiply histograms with different axis limits");
   }
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;

//*-* Create Sumw2 if h1 or h2 have Sumw2 set
   if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();

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

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t b1,b2,w,d1,d2;
   d1 = c1*c1;
   d2 = c2*c2;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            b1  = h1->GetBinContent(bin);
            b2  = h2->GetBinContent(bin);
            w   = (c1*b1)*(c2*b2);
            SetBinContent(bin,w);
            fEntries++;
            if (fSumw2.fN) {
               Double_t e1 = h1->GetBinError(bin);
               Double_t e2 = h2->GetBinError(bin);
               fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1);
            }
         }
      }
   }
   Stat_t s[10];
   GetStats(s);
   PutStats(s);
}

//______________________________________________________________________________
void TH1::Paint(Option_t *option)
{
//*-*-*-*-*-*-*-*-*Control routine to paint any kind of histograms*-*-*-*-*-*-*
//*-*              ===============================================
//
//  This function is automatically called by TCanvas::Update.
//  (see TH1::Draw for the list of options)

   if (!fPainter) fPainter = TVirtualHistPainter::HistPainter(this);
   if (fPainter) {
      if (strlen(option) > 0) fPainter->Paint(option);
      else                    fPainter->Paint(fOption.Data());
   }
}

//______________________________________________________________________________
TH1 *TH1::Rebin(Int_t ngroup, const char*newname)
{
//*-*-*-*-*Rebin this histogram grouping ngroup bins together*-*-*-*-*-*-*-*-*
//*-*      ==================================================
//   if newname is not blank a new temporary histogram hnew is created.
//   else the current histogram is modified (default)
//   The parameter ngroup indicates how many bins of this have to me merged
//   into one bin of hnew
//   If the original histogram has errors stored (via Sumw2), the resulting
//   histograms has new errors correctly calculated.
//
//   examples: if h1 is an existing TH1F histogram with 100 bins
//     h1->Rebin();  //merges two bins in one in h1: previous contents of h1 are lost
//     h1->Rebin(5); //merges five bins in one in h1
//     TH1F *hnew = h1->Rebin(5,"hnew"); // creates a new histogram hnew
//                                       //merging 5 bins of h1 in one bin
//   NOTE1: This function is currently implemented only for 1-D histograms
//
//   NOTE2: If ngroup is not an exact divider of the number of bins,
//          the top limit of the rebinned histogram 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.

   Int_t nbins   = fXaxis.GetNbins();
   Axis_t xmin  = fXaxis.GetXmin();
   Axis_t xmax  = fXaxis.GetXmax();
   if ((ngroup <= 0) || (ngroup > nbins)) {
      Error("Rebin", "Illegal value of ngroup=%d",ngroup);
      return 0;
   }
   if (fDimension > 1 || InheritsFrom("TProfile")) {
      Error("Rebin", "Operation valid on 1-D histograms only");
      return 0;
   }
   Int_t newbins = nbins/ngroup;

   // Save old bin contents into a new array
   Double_t *oldBins = new Double_t[nbins];
   Int_t bin, i;
   for (bin=0;bin<nbins;bin++) oldBins[bin] = GetBinContent(bin+1);
   Double_t *oldErrors = 0;
   if (fSumw2.fN != 0) {
      oldErrors = new Double_t[nbins];
      for (bin=0;bin<nbins;bin++) oldErrors[bin] = GetBinError(bin+1);
   }

   // create a clone of the old histogram if newname is specified
   TH1 *hnew = this;
   if (strlen(newname) > 0) {
      hnew = (TH1*)Clone();
      hnew->SetName(newname);
   }

   // change axis specs and rebuild bin contents array
   if(newbins*ngroup != nbins) {
      xmax = fXaxis.GetBinUpEdge(newbins*ngroup);
      hnew->fTsumw = 0; //stats must be reset because top bins will be moved to overflow bin
   }
   hnew->SetBins(newbins,xmin,xmax); //this also changes errors array (if any)

   // copy merged bin contents (ignore under/overflows)
   Int_t oldbin = 0;
   Double_t binContent, binError;
   for (bin = 0;bin<=newbins;bin++) {
      binContent = 0;
      binError   = 0;
      for (i=0;i<ngroup;i++) {
         if (oldbin+i >= nbins) break;
         binContent += oldBins[oldbin+i];
         if (oldErrors) binError += oldErrors[oldbin+i]*oldErrors[oldbin+i];
      }
      hnew->SetBinContent(bin+1,binContent);
      if (oldErrors) hnew->SetBinError(bin+1,TMath::Sqrt(binError));
      oldbin += ngroup;
   }

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

//______________________________________________________________________________
void TH1::RebinAxis(Axis_t x, const char *ax)
{
// 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 algorithm works for 1-d, 2-d and 3-d histograms.
// 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;
   Axis_t cxmin = axis->GetXmin();
   Axis_t cxmax = axis->GetXmax();
   Int_t  nbinsx = fXaxis.GetNbins();
   Int_t  nbinsy = fYaxis.GetNbins();
   Int_t  nbinsz = fZaxis.GetNbins();
   Axis_t range = cxmax-cxmin;

    //recompute new axis limits by doubling the current range
   Int_t bin;
   if (x < cxmin) {
      while (1) {
         range *= 2;
         if (x < cxmax-range) continue;
         xmin = cxmin - range/4;
         xmax = xmin + range;
         if (x < xmin) {
            xmin = cxmax - range;
            xmax = cxmax;
         }
         if ( xmin < 0 && x >= 0) {
            xmin = 0;
            xmax = range;
         }
         if (xmax >= 0 && cxmax <= 0) {
            xmax = 0;
            xmin = xmax - range;
         }
         break;
      }
   } else {
      while (1) {
         range *= 2;
         if (x >= cxmin+range) continue;
         xmax = cxmax + range/4;
         xmin = xmax - range;
         if (x >= xmax) {
            xmax = cxmin + range;
            xmin = cxmin;
         }
         if ( xmax > 0 && x < 0) {
            xmax = 0;
            xmin = -range;
         }
         if (xmin < 0 && cxmin >= 0) {
            xmin = 0;
            xmax = xmin + range;
         }
         break;
      }
   }

   //save a copy of this histogram
   TH1 *hold = (TH1*)Clone();
   hold->SetDirectory(0);

   //set new axis limits
   axis->SetLimits(xmin,xmax);

   //now loop on all bins and refill
   Double_t err,cu;
   Axis_t bx,by,bz;
   Int_t errors = GetSumw2N();
   Int_t ix,iy,iz,ibin,binx,biny,binz;
   Reset("ICE"); //reset only Integral, contents and Errors
   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);
            bin = hold->GetBin(binx,biny,binz);
            ibin= GetBin(ix,iy,iz);
            cu  = hold->GetBinContent(bin);
            AddBinContent(ibin,cu);
            if (errors) {
               err = hold->GetBinError(bin);
               fSumw2.fArray[ibin] += err*err;
            }
         }
      }
   }
   delete hold;
}

//______________________________________________________________________________
void TH1::Scale(Double_t c1)
{
//*-*-*-*-*Multiply this histogram by a constant c1*-*-*-*-*-*-*-*-*
//*-*      ========================================
//
//   this = c1*this
//
// Note that both contents and errors(if any) are scaled.
// This function uses the services of TH1::Add
//

   Double_t ent = fEntries;
   Add(this,this,c1,0);
   fEntries = ent;

   //if contours set, must also scale contours
   Int_t ncontours = GetContour();
   if (ncontours == 0) return;
   Double_t *levels = fContour.GetArray();
   for (Int_t i=0;i<ncontours;i++) {
      levels[i] *= c1;
   }
}
// -------------------------------------------------------------------------
void  TH1::SmoothArray(Int_t NN, Double_t *XX, Int_t ntimes)
{
// smooth array XX, translation of Hbook routine hsmoof.F
// based on algorithm 353QH twice presented by J. Friedman
// in Proc.of the 1974 CERN School of Computing, Norway, 11-24 August, 1974.

   Int_t ii, jj, ik, jk, kk, nn1, nn2;
   Double_t hh[5];
   Double_t *YY = new Double_t[NN];
   Double_t *ZZ = new Double_t[NN];
   Double_t *RR = new Double_t[NN];

   for (Int_t pass=0;pass<ntimes;pass++) {
      // first copy original data into temp array

      for (ii = 0; ii < NN; ii++) {
         YY[ii] = XX[ii];
      }

//  do 353 i.e. running median 3, 5, and 3 in a single loop
      for  (kk = 0; kk < 3; kk++)  {
         ik = 0;
         if  (kk == 1)  ik = 1;
         nn1 = ik + 1;
         nn2 = NN - ik - 1;
       // do all elements beside the first and last point for median 3
       //  and first two and last 2 for median 5
         for  (ii = nn1; ii < nn2; ii++)  {
            for  (jj = 0; jj < 3; jj++)   {
               hh[jj] = YY[ii + jj - 1];
            }
            ZZ[ii] = TH1::SmoothMedian(3 + 2*ik, hh);
         }

         if  (kk == 0)  {   // first median 3
// first point
            hh[0] = 3*YY[1] - 2*YY[2];
            hh[1] = YY[0];
            hh[2] = YY[2];
            ZZ[0] = TH1::SmoothMedian(3, hh);
// last point
            hh[0] = YY[NN - 2];
            hh[1] = YY[NN - 1];
            hh[2] = 3*YY[NN - 2] - 2*YY[NN - 3];
            ZZ[NN - 1] = TH1::SmoothMedian(3, hh);
         }
         if  (kk == 1)  {   //  median 5
  //  first point remains the same
            ZZ[0] = YY[0];
            for  (ii = 0; ii < 3; ii++) {
               hh[ii] = YY[ii];
            }
            ZZ[1] = TH1::SmoothMedian(3, hh);
// last two points
            for  (ii = 0; ii < 3; ii++) {
               hh[ii] = YY[NN - 3 + ii];
            }
            ZZ[NN - 2] = TH1::SmoothMedian(3, hh);
            ZZ[NN - 1] = YY[NN - 1];
         }
      }

// quadratic interpolation for flat segments
      nn2 = nn2 - 1;
      for (ii = nn1 + 1; ii < nn2; ii++) {
         if  (ZZ[ii - 1] != ZZ[ii]) continue;
         if  (ZZ[ii] != ZZ[ii + 1]) continue;
         hh[0] = ZZ[ii - 2] - ZZ[ii];
         hh[1] = ZZ[ii + 2] - ZZ[ii];
         if  (hh[0] * hh[1] < 0) continue;
         jk = 0;
         if  ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
         YY[ii] = 0.5*ZZ[ii - 2*jk] + ZZ[ii - jk]/0.75 + ZZ[ii + 2*jk] /6.;
         YY[ii + jk] = 0.5*(ZZ[ii + 2*jk] - ZZ[ii - 2*jk]) + ZZ[ii - jk];
      }

// running means
      for  (ii = 1; ii < NN - 1; ii++) {
         RR[ii] = 0.25*YY[ii - 1] + 0.5*YY[ii] + 0.25*YY[ii + 1];
      }
      RR[0] = YY[0];
      RR[NN - 1] = YY[NN - 1];

// now do the same for residuals

      for  (ii = 0; ii < NN; ii++)  {
         YY[ii] = XX[ii] - RR[ii];
      }

//  do 353 i.e. running median 3, 5, and 3 in a single loop
      for  (kk = 0; kk < 3; kk++)  {
         ik = 0;
         if  (kk == 1)  ik = 1;
         nn1 = ik + 1;
         nn2 = NN - ik - 1;
       // do all elements beside the first and last point for median 3
       //  and first two and last 2 for median 5
         for  (ii = nn1; ii < nn2; ii++)  {
            for  (jj = 0; jj < 3; jj++) {
               hh[jj] = YY[ii + jj - 1];
            }
            ZZ[ii] = TH1::SmoothMedian(3 + 2*ik, hh);
         }

         if  (kk == 0)  {   // first median 3
// first point
            hh[0] = 3*YY[1] - 2*YY[2];
            hh[1] = YY[0];
            hh[2] = YY[2];
            ZZ[0] = TH1::SmoothMedian(3, hh);
// last point
            hh[0] = YY[NN - 2];
            hh[1] = YY[NN - 1];
            hh[2] = 3*YY[NN - 2] - 2*YY[NN - 3];
            ZZ[NN - 1] = TH1::SmoothMedian(3, hh);
         }
         if  (kk == 1)  {   //  median 5
//  first point remains the same
            ZZ[0] = YY[0];
            for  (ii = 0; ii < 3; ii++) {
               hh[ii] = YY[ii];
            }
            ZZ[1] = TH1::SmoothMedian(3, hh);
// last two points
            for  (ii = 0; ii < 3; ii++) {
               hh[ii] = YY[NN - 3 + ii];
            }
            ZZ[NN - 2] = TH1::SmoothMedian(3, hh);
            ZZ[NN - 1] = YY[NN - 1];
         }
      }

// quadratic interpolation for flat segments
      nn2 = nn2 - 1;
      for (ii = nn1 + 1; ii < nn2; ii++) {
         if  (ZZ[ii - 1] != ZZ[ii]) continue;
         if  (ZZ[ii] != ZZ[ii + 1]) continue;
         hh[0] = ZZ[ii - 2] - ZZ[ii];
         hh[1] = ZZ[ii + 2] - ZZ[ii];
         if  (hh[0] * hh[1] < 0) continue;
         jk = 0;
         if  ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
         YY[ii] = 0.5*ZZ[ii - 2*jk] + ZZ[ii - jk]/0.75 + ZZ[ii + 2*jk]/6.;
         YY[ii + jk] = 0.5*(ZZ[ii + 2*jk] - ZZ[ii - 2*jk]) + ZZ[ii - jk];
      }

// running means
      for  (ii = 1; ii < NN - 1; ii++) {
         ZZ[ii] = 0.25*YY[ii - 1] + 0.5*YY[ii] + 0.25*YY[ii + 1];
      }
      ZZ[0] = YY[0];
      ZZ[NN - 1] = YY[NN - 1];

//  add smoothed XX and smoothed residuals

      for  (ii = 0; ii < NN; ii++) {
         XX[ii] = RR[ii] + ZZ[ii];
      }
   }
   delete [] YY;
   delete [] ZZ;
   delete [] RR;
}

// ------------------------------------------------------------------------
Double_t  TH1::SmoothMedian(Int_t n, Double_t *a)
{
// return the median of a vector a in monotonic order with length n
// where median is a number which divides sequence of n numbers
// into 2 halves. When n is odd, the median is kth element k = (n + 1) / 2.
// when n is even the median is a mean of the elements k = n/2 and k = n/2 + 1.

  Int_t in, imin, imax;
  Double_t  xm;

  if  (n%2 == 0)  in = n / 2;
  else            in = n / 2 + 1;

  // find array element with maximum content
  imax = TMath::LocMax(n,a);
  xm = a[imax];

  while (in < n) {
     imin = TMath::LocMin(n,a);  // find array element with minimum content
     a[imin] = xm;
     in++;
  }
  imin = TMath::LocMin(n,a);
  return a[imin];
}


// ------------------------------------------------------------------------
void  TH1::Smooth(Int_t ntimes)
{
// Smooth bin contents of this histogram.
// bin contents are replaced by their smooth values.
// Errors (if any) are not modified.
// algorithm can only be applied to 1-d histograms

   if (fDimension != 1) {
      Error("Smooth","Smooth only supported for 1-d histograms");
      return;
   }
   Int_t nbins = fXaxis.GetNbins();
   Double_t *XX = new Double_t[nbins];
   Int_t i;
   for (i=0;i<nbins;i++) {
      XX[i] = GetBinContent(i+1);
   }

   TH1::SmoothArray(nbins,XX,ntimes);

   for (i=0;i<nbins;i++) {
      SetBinContent(i+1,XX[i]);
   }
   if (gPad) gPad->Modified();
}


//_______________________________________________________________________
void TH1::Streamer(TBuffer &b)
{
//*-*-*-*-*-*-*-*-*Stream a class object*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*              =====================
   UInt_t R__s, R__c;
   if (b.IsReading()) {
      Version_t R__v = b.ReadVersion(&R__s, &R__c);
      TNamed::Streamer(b);
      TAttLine::Streamer(b);
      TAttFill::Streamer(b);
      TAttMarker::Streamer(b);
      b >> fNcells;
      fXaxis.Streamer(b);
      fYaxis.Streamer(b);
      fZaxis.Streamer(b);
      fXaxis.SetParent(this);
      fYaxis.SetParent(this);
      fZaxis.SetParent(this);
      b >> fBarOffset;
      b >> fBarWidth;
      b >> fEntries;
      b >> fTsumw;
      b >> fTsumw2;
      b >> fTsumwx;
      b >> fTsumwx2;
      if (R__v < 2) {
         Float_t maximum, minimum, norm;
         Float_t *contour=0;
         b >> maximum; fMaximum = maximum;
         b >> minimum; fMinimum = minimum;
         b >> norm;    fNormFactor = norm;
         Int_t n = b.ReadArray(contour);
         fContour.Set(n);
         for (Int_t i=0;i<n;i++) fContour.fArray[i] = contour[i];
         delete [] contour;
      } else {
         b >> fMaximum;
         b >> fMinimum;
         b >> fNormFactor;
         fContour.Streamer(b);
      }
      fSumw2.Streamer(b);
      fOption.Streamer(b);
      fFunctions->Delete();
      fFunctions->Streamer(b);
      if (!gROOT->ReadingObject()) {
         fDirectory = gDirectory;
         if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
      }
      b.CheckByteCount(R__s, R__c, TH1::IsA());
   } else {
      R__c = b.WriteVersion(TH1::IsA(), kTRUE);
      TNamed::Streamer(b);
      TAttLine::Streamer(b);
      TAttFill::Streamer(b);
      TAttMarker::Streamer(b);
      b << fNcells;
      fXaxis.Streamer(b);
      fYaxis.Streamer(b);
      fZaxis.Streamer(b);
      b << fBarOffset;
      b << fBarWidth;
      b << fEntries;
      b << fTsumw;
      b << fTsumw2;
      b << fTsumwx;
      b << fTsumwx2;
      b << fMaximum;
      b << fMinimum;
      b << fNormFactor;
      fContour.Streamer(b);
      fSumw2.Streamer(b);
      fOption.Streamer(b);
      fFunctions->Streamer(b);
      b.SetByteCount(R__c, kTRUE);
   }
}

//______________________________________________________________________________
void TH1::Print(Option_t *option)
{
//*-*-*-*-*-*-*Print some global quantities for this histogram*-*-*-*-*-*-*-*
//*-*          ===============================================
//
//  If option "range" is given, bin contents and errors are also printed
//                     for all bins in the current range (default 1-->nbins)
//  If option "all" is given, bin contents and errors are also printed
//                     for all bins including under and overflows.
//
   printf( "TH1.Print Name= %s, Entries= %d, Total sum= %gn",GetName(),Int_t(fEntries),GetSumOfWeights());
   TString opt = option;
   opt.ToLower();
   Int_t all;
   if      (opt.Contains("all"))   all = 0;
   else if (opt.Contains("range")) all = 1;
   else                            return;

   Int_t bin, binx, biny, binz;
   Int_t firstx=0,lastx=0,firsty=0,lasty=0,firstz=0,lastz=0;
   if (all == 0) {
      lastx  = fXaxis.GetNbins()+1;
      if (fDimension > 1) lasty = fYaxis.GetNbins()+1;
      if (fDimension > 2) lastz = fZaxis.GetNbins()+1;
   } else {
      firstx = fXaxis.GetFirst(); lastx  = fXaxis.GetLast();
      if (fDimension > 1) {firsty = fYaxis.GetFirst(); lasty = fYaxis.GetLast();}
      if (fDimension > 2) {firstz = fZaxis.GetFirst(); lastz = fZaxis.GetLast();}
   }
   Stat_t w,e;
   Axis_t x,y,z;
   if (fDimension == 1) {
      for (binx=firstx;binx<=lastx;binx++) {
         x = fXaxis.GetBinCenter(binx);
         w = GetBinContent(binx);
         e = GetBinError(binx);
         if(fSumw2.fN) printf(" fSumw[%d]=%g, x=%g, error=%gn",binx,w,x,e);
         else          printf(" fSumw[%d]=%g, x=%gn",binx,w,x);
      }
   }
   if (fDimension == 2) {
      for (biny=firsty;biny<=lasty;biny++) {
         y = fYaxis.GetBinCenter(biny);
         for (binx=firstx;binx<=lastx;binx++) {
            bin = GetBin(binx,biny);
            x = fXaxis.GetBinCenter(binx);
            w = GetBinContent(bin);
            e = GetBinError(bin);
            if(fSumw2.fN) printf(" fSumw[%d][%d]=%g, x=%g, y=%g, error=%gn",binx,biny,w,x,y,e);
            else          printf(" fSumw[%d][%d]=%g, x=%g, y=%gn",binx,biny,w,x,y);
         }
      }
   }
   if (fDimension == 3) {
      for (binz=firstz;binz<=lastz;binz++) {
         z = fZaxis.GetBinCenter(binz);
         for (biny=firsty;biny<=lasty;biny++) {
            y = fYaxis.GetBinCenter(biny);
            for (binx=firstx;binx<=lastx;binx++) {
               bin = GetBin(binx,biny,binz);
               x = fXaxis.GetBinCenter(binx);
               w = GetBinContent(bin);
               e = GetBinError(bin);
               if(fSumw2.fN) printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%g, error=%gn",binx,biny,binz,w,x,y,z,e);
               else          printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%gn",binx,biny,binz,w,x,y,z);
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Reset(Option_t *option)
{
//*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-*
//*-*            ===========================================
//
// if option "ICE" is specified, resets only Integral, Contents and Errors

   TString opt = option;
   opt.ToUpper();
   fSumw2.Reset();
   if (fIntegral) {delete [] fIntegral; fIntegral = 0;}

   if (opt.Contains("ICE")) return;
   fTsumw       = 0;
   fTsumw2      = 0;
   fTsumwx      = 0;
   fTsumwx2     = 0;
   fEntries     = 0;

   fFunctions->Delete();
   fContour.Set(0);
}

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

   char quote = '"';
   out<<"   "<<endl;

   SaveFillAttributes(out,GetName(),0,1001);
   SaveLineAttributes(out,GetName(),1,1,1);
   SaveMarkerAttributes(out,GetName(),1,1,1);
   fXaxis.SaveAttributes(out,GetName(),"->GetXaxis()");
   fYaxis.SaveAttributes(out,GetName(),"->GetYaxis()");
   fZaxis.SaveAttributes(out,GetName(),"->GetZaxis()");
   out<<"   "<<GetName()<<"->Draw("
      <<quote<<option<<quote<<");"<<endl;
}

//______________________________________________________________________________
void TH1::UseCurrentStyle()
{
//*-*-*-*-*-*Replace current attributes by current style*-*-*-*-*
//*-*        ===========================================

   fXaxis.ResetAttAxis("X");
   fYaxis.ResetAttAxis("Y");
   SetBarOffset(gStyle->GetBarOffset());
   SetBarWidth(gStyle->GetBarWidth());
   SetFillColor(gStyle->GetHistFillColor());
   SetFillStyle(gStyle->GetHistFillStyle());
   SetLineColor(gStyle->GetHistLineColor());
   SetLineStyle(gStyle->GetHistLineStyle());
   SetLineWidth(gStyle->GetHistLineWidth());
   SetMarkerColor(gStyle->GetMarkerColor());
   SetMarkerStyle(gStyle->GetMarkerStyle());
   SetMarkerSize(gStyle->GetMarkerSize());
}

//______________________________________________________________________________
Stat_t TH1::GetMean(Int_t axis)
{
//*-*-*-*-*-*-*-*Return mean value of this histogram along the X axis*-*-*-*-*
//*-*            ====================================================

  if (axis <1 || axis > 3) return 0;
  Stat_t stats[10];
  GetStats(stats);
  if (stats[0] == 0) return 0;
  Int_t ax[3] = {2,4,7};
  return stats[ax[axis-1]]/stats[0];
}

//______________________________________________________________________________
Stat_t TH1::GetRMS(Int_t axis)
{
//*-*-*-*-*-*-*-*Return the Root Mean Square value of this histogram*-*-*-*-*
//*-*            ===================================================

  if (axis <1 || axis > 3) return 0;
  Stat_t x, rms2, stats[10];
  GetStats(stats);
  if (stats[0] == 0) return 0;
  Int_t ax[3] = {2,4,7};
  Int_t axm = ax[axis-1];
  x    = stats[axm]/stats[0];
  rms2 = TMath::Abs(stats[axm+1]/stats[0] -x*x);
  return TMath::Sqrt(rms2);
}

//______________________________________________________________________________
void TH1::GetStats(Stat_t *stats)
{
   // fill the array stats from the contents of this histogram
   // The array stats must be correctly dimensionned in the calling program.
   // stats[0] = sumw
   // stats[1] = sumw2
   // stats[2] = sumwx
   // stats[3] = sumwx2
   //
   // 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.

   // Loop on bins (possibly including underflows/overflows)
   Int_t bin, binx;
   Stat_t w;
   Axis_t x;
   if (fTsumw == 0 || fXaxis.TestBit(TAxis::kAxisRange)) {
      for (bin=0;bin<4;bin++) stats[bin] = 0;
      for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
         x = fXaxis.GetBinCenter(binx);
         w = TMath::Abs(GetBinContent(binx));
         stats[0] += w;
         stats[1] += w*w;
         stats[2] += w*x;
         stats[3] += w*x*x;
      }
   } else {
      stats[0] = fTsumw;
      stats[1] = fTsumw2;
      stats[2] = fTsumwx;
      stats[3] = fTsumwx2;
   }
}

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

   fTsumw   = stats[0];
   fTsumw2  = stats[1];
   fTsumwx  = stats[2];
   fTsumwx2 = stats[3];
}

//______________________________________________________________________________
Stat_t TH1::GetSumOfWeights()
{
//*-*-*-*-*-*-*-*Return the sum of weights excluding under/overflows*-*-*-*-*
//*-*            ===================================================
  Int_t bin,binx,biny,binz;
  Stat_t sum =0;
  for(binz=1; binz<=fZaxis.GetNbins(); binz++) {
     for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
        for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
           bin = GetBin(binx,biny,binz);
           sum += GetBinContent(bin);
        }
     }
  }
  return sum;
}


//______________________________________________________________________________
Stat_t TH1::Integral()
{
//Return integral of bin contents. Only bins in the bins range are considered.

   return Integral(fXaxis.GetFirst(),fXaxis.GetLast());
}

//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2)
{
//Return integral of bin contents between binx1 and binx2 for a 1-D histogram

   if (binx1 < 0) binx1 = 0;
   Stat_t integral = 0;

//*-*- Loop on bins in specified range
   Int_t binx;
   for (binx=binx1;binx<=binx2;binx++) {
      integral += GetBinContent(binx);
   }
   return integral;
}

//______________________________________________________________________________
Double_t TH1::KolmogorovTest(TH1 *h2, Option_t *option)
{
//  Statistical test of compatibility in shape between
//  THIS histogram and h2, using Kolmogorov test.
//     Default: Ignore under- and overflow bins in comparison
//
//     option is a character string to specify options
//         "U" include Underflows in test  (also for 2-dim)
//         "O" include Overflows     (also valid for 2-dim)
//         "N" include comparison of normalizations
//         "D" Put out a line of "Debug" printout
//         "M" Return the Maximum Kolmogorov distance instead of prob
//
//   The returned function value is the probability of test
//       (much less than one means NOT compatible)
//
//  Code adapted by Rene Brun from original HBOOK routine HDIFF

   TString opt = option;
   opt.ToUpper();

   Double_t prb = 0;
   TH1 *h1 = this;
   if (h2 == 0) return 0;
   TAxis *axis1 = h1->GetXaxis();
   TAxis *axis2 = h2->GetXaxis();
   Int_t ncx1   = axis1->GetNbins();
   Int_t ncx2   = axis2->GetNbins();

     // Check consistency of dimensions
   if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
      Error("KolmogorovTest","Histograms must be 1-Dn");
      return 0;
   }

     // Check consistency in number of channels
   if (ncx1 != ncx2) {
      Error("KolmogorovTest","Number of channels is different, %d and %dn",ncx1,ncx2);
      return 0;
   }

     // Check consistency in channel edges
   Double_t difprec = 1e-5;
   Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
   Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
   if (diff1 > difprec || diff2 > difprec) {
      Error("KolmogorovTest","histograms with different binning");
      return 0;
   }

   Bool_t afunc1 = kFALSE;
   Bool_t afunc2 = kFALSE;
   Double_t sum1 = 0, sum2 = 0;
   Int_t bin;
   for (bin=1;bin<=ncx1;bin++) {
      sum1 += h1->GetBinContent(bin);
      sum2 += h2->GetBinContent(bin);
   }
   if (sum1 == 0) {
      Error("KolmogorovTest","Histogram1 %s integral is zeron",h1->GetName());
      return 0;
   }
   if (sum2 == 0) {
      Error("KolmogorovTest","Histogram2 %s integral is zeron",h2->GetName());
      return 0;
   }
   Double_t tsum1 = sum1;
   Double_t tsum2 = sum2;
   if (opt.Contains("U")) {
      tsum1 += h1->GetBinContent(0);
      tsum2 += h2->GetBinContent(0);
   }
   if (opt.Contains("O")) {
      tsum1 += h1->GetBinContent(ncx1+1);
      tsum2 += h2->GetBinContent(ncx1+1);
   }

     // Check if histograms are weighted.
     // If number of entries = number of channels, probably histograms were
     // not filled via Fill(), but via SetBinContent()
   Double_t ne1 = h1->GetEntries();
   Double_t ne2 = h2->GetEntries();
     // look at first histogram
   Double_t difsum1 = (ne1-tsum1)/tsum1;
   Double_t esum1 = sum1;
   if (difsum1 > difprec && Int_t(ne1) != ncx1) {
      if (opt.Contains("U") || opt.Contains("O")) {
         Warning("KolmogorovTest","U/O option with weighted events for hist:%sn",h1->GetName());
      }
      if (h1->GetSumw2N() == 0) {
         Warning("KolmogorovTest","Weighted events and no Sumw2, hist:%sn",h1->GetName());
      } else {
         esum1 = h1->GetSumOfWeights();
      }
   }
     // look at second histogram
   Double_t difsum2 = (ne2-tsum2)/tsum2;
   Double_t esum2   = sum2;
   if (difsum2 > difprec && Int_t(ne2) != ncx1) {
      if (opt.Contains("U") || opt.Contains("O")) {
         Warning("KolmogorovTest","U/O option with weighted events for hist:%sn",h2->GetName());
      }
      if (h2->GetSumw2N() == 0) {
         Warning("KolmogorovTest","Weighted events and no Sumw2, hist:%sn",h2->GetName());
      } else {
         esum1 = h2->GetSumOfWeights();
      }
   }

   Double_t s1 = 1/sum1;
   Double_t s2 = 1/sum2;

      // Find largest difference for Kolmogorov Test
   Double_t dfmax =0, rsum1 = 0, rsum2 = 0;

   Int_t first = 1;
   Int_t last  = ncx1;
   if (opt.Contains("U")) first = 0;
   if (opt.Contains("O")) last  = ncx1+1;
   for (bin=first;bin<=last;bin++) {
      rsum1 += s1*h1->GetBinContent(bin);
      rsum2 += s2*h2->GetBinContent(bin);
      dfmax = TMath::Max(dfmax,TMath::Abs(rsum1-rsum2));
   }

      // Get Kolmogorov probability
   Double_t z, prb1=0, prb2=0;
   if (afunc1)      z = dfmax*TMath::Sqrt(esum2);
   else if (afunc2) z = dfmax*TMath::Sqrt(esum1);
   else             z = dfmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));

   prb = TMath::KolmogorovProb(z);

   if (opt.Contains("N")) {
      // Combine probabilities for shape and normalization,
      prb1 = prb;
      Double_t resum1 = esum1;  if (afunc1) resum1 = 0;
      Double_t resum2 = esum2;  if (afunc2) resum2 = 0;
      Double_t d12    = esum1-esum2;
      Double_t chi2   = d12*d12/(resum1+resum2);
      prb2 = TMath::Prob(chi2,1);
      // see Eadie et al., section 11.6.2
      if (prb > 0 && prb2 > 0) prb *= prb2*(1-TMath::Log(prb*prb2));
      else                     prb = 0;
   }
      // debug printout
   if (opt.Contains("D")) {
      printf(" Kolmo Prob  h1 = %s, sum1=%gn",h1->GetName(),sum1);
      printf(" Kolmo Prob  h2 = %s, sum2=%gn",h2->GetName(),sum2);
      printf(" Kolmo Probabil = %g, Max Dist = %gn",prb,dfmax);
      if (opt.Contains("N"))
      printf(" Kolmo Probabil = %f for shape alone, =%f for normalisation alonen",prb1,prb2);
   }
      // This numerical error condition should never occur:
   if (TMath::Abs(rsum1-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h1=%sn",h1->GetName());
   if (TMath::Abs(rsum2-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h2=%sn",h2->GetName());

   if(opt.Contains("M"))  return dfmax;
   else                   return prb;
}

//______________________________________________________________________________
void TH1::SetContent(Stat_t *content)
{
//*-*-*-*-*-*-*-*Replace bin contents by the contents of array content*-*-*-*
//*-*            =====================================================
   Int_t bin;
   Stat_t bincontent;
   for (bin=0; bin<fNcells; bin++) {
      bincontent = *(content + bin);
      SetBinContent(bin, bincontent);
   }
}

//______________________________________________________________________________
Int_t TH1::GetContour(Double_t *levels)
{
//  Return contour values into array levels if pointer levels is non zero
//
//  The function returns the number of contour levels.
//  see GetContourLevel to return one contour only
//

  Int_t nlevels = fContour.fN;
  if (levels) {
     if (nlevels == 0) {
        nlevels = 20;
        SetContour(nlevels);
     } else {
        if (TestBit(kUserContour) == 0) SetContour(nlevels);
     }
     for (Int_t level=0; level<nlevels; level++) levels[level] = fContour.fArray[level];
  }
  return nlevels;
}

//______________________________________________________________________________
Double_t TH1::GetContourLevel(Int_t level)
{
// Return value of contour number level
// see GetContour to return the array of all contour levels

  if (level <0 || level >= fContour.fN) return 0;
  Double_t zlevel = fContour.fArray[level];
  return zlevel;
}


//______________________________________________________________________________
void TH1::SetContour(Int_t  nlevels, Double_t *levels)
{
//*-*-*-*-*-*-*-*Set the number and values of contour levels*-*-*-*-*-*-*-*-*
//*-*            ===========================================
//
//  By default the number of contour levels is set to 20.
//
//  if argument levels = 0 or issing, equidistant contours are computed
//

  Int_t level;
  ResetBit(kUserContour);
  if (nlevels <=0 ) {
     fContour.Set(0);
     return;
  }
  fContour.Set(nlevels);

//*-*-  Contour levels are specified
  if (levels) {
     SetBit(kUserContour);
     for (level=0; level<nlevels; level++) fContour.fArray[level] = levels[level];
  } else {
//*-*- contour levels are computed automatically as equidistant contours
     Double_t zmin = GetMinimum();
     Double_t zmax = GetMaximum();
     Double_t dz   = (zmax-zmin)/Double_t(nlevels);
     if (gPad && gPad->GetLogz()) {
        if (zmax <= 0) return;
        if (zmin <= 0) zmin = 0.001*zmax;
        zmin = TMath::Log10(zmin);
        zmax = TMath::Log10(zmax);
        dz   = (zmax-zmin)/Double_t(nlevels);
     }
     for (level=0; level<nlevels; level++) {
        fContour.fArray[level] = zmin + dz*Double_t(level);
     }
  }
}

//______________________________________________________________________________
void TH1::SetContourLevel(Int_t level, Double_t value)
{
//*-*-*-*-*-*-*-*-*-*-*Set value for one contour level*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ===============================
  if (level <0 || level >= fContour.fN) return;
  SetBit(kUserContour);
  fContour.fArray[level] = value;
}

//______________________________________________________________________________
Double_t TH1::GetMaximum()
{
//*-*-*-*-*-*-*-*-*-*-*Return maximum value of bins in the range*-*-*-*-*-*
//*-*                  =========================================
  if (fMaximum != -1111) return fMaximum;
  Int_t bin, binx, biny, binz;
  Int_t xfirst  = fXaxis.GetFirst();
  Int_t xlast   = fXaxis.GetLast();
  Int_t yfirst  = fYaxis.GetFirst();
  Int_t ylast   = fYaxis.GetLast();
  Int_t zfirst  = fZaxis.GetFirst();
  Int_t zlast   = fZaxis.GetLast();
  Double_t maximum = -FLT_MAX, value;
  for (binz=zfirst;binz<=zlast;binz++) {
     for (biny=yfirst;biny<=ylast;biny++) {
        for (binx=xfirst;binx<=xlast;binx++) {
           bin = GetBin(binx,biny,binz);
           value = GetBinContent(bin);
           if (value > maximum) maximum = value;
        }
     }
  }
  return maximum;
}

//______________________________________________________________________________
Int_t TH1::GetMaximumBin()
{
//*-*-*-*-*-*-*Return location of bin with maximum value in the range*-*
//*-*          ======================================================
  Int_t locmax, locmay, locmaz;
  return GetMaximumBin(locmax, locmay, locmaz);
}

//______________________________________________________________________________
Int_t TH1::GetMaximumBin(Int_t &locmax, Int_t &locmay, Int_t &locmaz)
{
//*-*-*-*-*-*-*Return location of bin with maximum value in the range*-*
//*-*          ======================================================
  Int_t bin, binx, biny, binz;
  Int_t locm;
  Int_t xfirst  = fXaxis.GetFirst();
  Int_t xlast   = fXaxis.GetLast();
  Int_t yfirst  = fYaxis.GetFirst();
  Int_t ylast   = fYaxis.GetLast();
  Int_t zfirst  = fZaxis.GetFirst();
  Int_t zlast   = fZaxis.GetLast();
  Double_t maximum = -FLT_MAX, value;
  locm = locmax = locmay = locmaz = 0;
  for (binz=zfirst;binz<=zlast;binz++) {
     for (biny=yfirst;biny<=ylast;biny++) {
        for (binx=xfirst;binx<=xlast;binx++) {
           bin = GetBin(binx,biny,binz);
           value = GetBinContent(bin);
           if (value > maximum) {
              maximum = value;
              locm    = bin;
              locmax  = binx;
              locmay  = biny;
              locmaz  = binz;
           }
        }
     }
  }
  return locm;
}

//______________________________________________________________________________
Double_t TH1::GetMinimum()
{
//*-*-*-*-*-*-*-*-*-*-*Return minimum value of bins in the range-*-*-*-*-*
//*-*                  =========================================
  if (fMinimum != -1111) return fMinimum;
  Int_t bin, binx, biny, binz;
  Int_t xfirst  = fXaxis.GetFirst();
  Int_t xlast   = fXaxis.GetLast();
  Int_t yfirst  = fYaxis.GetFirst();
  Int_t ylast   = fYaxis.GetLast();
  Int_t zfirst  = fZaxis.GetFirst();
  Int_t zlast   = fZaxis.GetLast();
  Double_t minimum=FLT_MAX, value;
  for (binz=zfirst;binz<=zlast;binz++) {
     for (biny=yfirst;biny<=ylast;biny++) {
        for (binx=xfirst;binx<=xlast;binx++) {
           bin = GetBin(binx,biny,binz);
           value = GetBinContent(bin);
           if (value < minimum) minimum = value;
        }
     }
  }
  return minimum;
}

//______________________________________________________________________________
Int_t TH1::GetMinimumBin()
{
//*-*-*-*-*-*-*Return location of bin with minimum value in the range*-*
//*-*          ======================================================
  Int_t locmix, locmiy, locmiz;
  return GetMinimumBin(locmix, locmiy, locmiz);
}

//______________________________________________________________________________
Int_t TH1::GetMinimumBin(Int_t &locmix, Int_t &locmiy, Int_t &locmiz)
{
//*-*-*-*-*-*-*Return location of bin with minimum value in the range*-*
//*-*          ======================================================
  Int_t bin, binx, biny, binz;
  Int_t locm;
  Int_t xfirst  = fXaxis.GetFirst();
  Int_t xlast   = fXaxis.GetLast();
  Int_t yfirst  = fYaxis.GetFirst();
  Int_t ylast   = fYaxis.GetLast();
  Int_t zfirst  = fZaxis.GetFirst();
  Int_t zlast   = fZaxis.GetLast();
  Double_t minimum = FLT_MAX, value;
  locm = locmix = locmiy = locmiz = 0;
  for (binz=zfirst;binz<=zlast;binz++) {
     for (biny=yfirst;biny<=ylast;biny++) {
        for (binx=xfirst;binx<=xlast;binx++) {
           bin = GetBin(binx,biny,binz);
           value = GetBinContent(bin);
           if (value < minimum) {
              minimum = value;
              locm    = bin;
              locmix  = binx;
              locmiy  = biny;
              locmiz  = binz;
           }
        }
     }
  }
  return locm;
}

//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Axis_t xmin, Axis_t xmax)
{
//*-*-*-*-*-*-*-*-*Redefine  x axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ===========================
// The X axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange

   if (GetDimension() != 1) {
      Error("SetBins","Operation only valid for 1-d histograms");
      return;
   }
   fXaxis.Set(nx,xmin,xmax);
   fNcells = nx+2;
   SetBinsLength(fNcells);
   if (fSumw2.fN) {
      fSumw2.Set(fNcells);
   }
}

//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Axis_t xmin, Axis_t xmax, Int_t ny, Axis_t ymin, Axis_t ymax)
{
//*-*-*-*-*-*-*-*-*Redefine  x and y axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-*              =================================
// The X and Y axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange

   if (GetDimension() != 2) {
      Error("SetBins","Operation only valid for 2-d histograms");
      return;
   }
   fXaxis.Set(nx,xmin,xmax);
   fYaxis.Set(ny,ymin,ymax);
   fNcells = (nx+2)*(ny+2);
   SetBinsLength(fNcells);
   if (fSumw2.fN) {
      fSumw2.Set(fNcells);
   }
}

//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Axis_t xmin, Axis_t xmax, Int_t ny, Axis_t ymin, Axis_t ymax, Int_t nz, Axis_t zmin, Axis_t zmax)
{
//*-*-*-*-*-*-*-*-*Redefine  x, y and z axis parameters*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ====================================
// The X, Y and Z axis parameters are modified.
// The bins content array is resized
// if errors (Sumw2) the errors array is resized
// The previous bin contents are lost
// To change only the axis limits, see TAxis::SetRange

   if (GetDimension() != 3) {
      Error("SetBins","Operation only valid for 3-d histograms");
      return;
   }
   fXaxis.Set(nx,xmin,xmax);
   fYaxis.Set(ny,ymin,ymax);
   fZaxis.Set(nz,zmin,zmax);
   fNcells = (nx+2)*(ny+2)*(nz+2);
   SetBinsLength(fNcells);
   if (fSumw2.fN) {
      fSumw2.Set(fNcells);
   }
}

//______________________________________________________________________________
void TH1::SetMaximum(Double_t maximum)
{
//*-*-*-*-*-*-*-*-*Set the maximum value for the Y axis*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ====================================
// By default the maximum value is automatically set to the maximum
// bin content plus a margin of 10 per cent.
// Use TH1::GetMaximum to find the maximum value of an histogram
// Use TH1::GetMaximumBin to find the bin with the maximum value of an histogram
//
   fMaximum = maximum;
}


//______________________________________________________________________________
void TH1::SetMinimum(Double_t minimum)
{
//*-*-*-*-*-*-*-*-*Set the minimum value for the Y axis*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ====================================
// By default the minimum value is automatically set to zero if all bin contents
// are positive or the minimum - 10 per cent otherwise.
// Use TH1::GetMinimum to find the minimum value of an histogram
// Use TH1::GetMinimumBin to find the bin with the minimum value of an histogram
//
   fMinimum = minimum;
}

//______________________________________________________________________________
void TH1::SetDirectory(TDirectory *dir)
{
   // By default when an histogram is created, it is added to the list
   // of histogram objects in the current directory in memory.
   // Remove reference to this histogram from current directory and add
   // reference to new directory dir. dir can be 0 in which case the
   // histogram does not belong to any directory.

   if (fDirectory == dir) return;
   if (fDirectory) fDirectory->GetList()->Remove(this);
   fDirectory = dir;
   if (fDirectory) fDirectory->GetList()->Add(this);
}


//______________________________________________________________________________
void TH1::SetError(Stat_t *error)
{
//*-*-*-*-*-*-*-*-*Replace bin errors by values in array error*-*-*-*-*-*-*-*-*
//*-*              ===========================================
  Int_t bin;
  Stat_t binerror;
  for (bin=0; bin<fNcells; bin++) {
     binerror = error[bin];
     SetBinError(bin, binerror);
  }
}

//______________________________________________________________________________
void TH1::SetName(const char *name)
{
// Change the name of this histogram
//

//  Histograms are named objects in a THashList.
//  We must update the hashlist if we change the name
   if (fDirectory) fDirectory->GetList()->Remove(this);
   fName = name;
   if (fDirectory) fDirectory->GetList()->Add(this);
}

//______________________________________________________________________________
void TH1::SetObject(const char *name, const char *title)
{
// Change the name and title of this histogram
//

//  Histograms are named objects in a THashList.
//  We must update the hashlist if we change the name
   if (fDirectory) fDirectory->GetList()->Remove(this);
   fName  = name;
   fTitle = title;
   if (fDirectory) fDirectory->GetList()->Add(this);
}

//______________________________________________________________________________
void TH1::SetStats(Bool_t stats)
{
//*-*-*-*-*-*-*-*-*Set statistics option on/off
//*-*              ============================
//  By default, the statistics box is drawn.
//  The paint options can be selected via gStyle->SetOptStats.
//  This function sets/resets the kNoStats bin in the histogram object.
//  It has priority over the Style option.

   ResetBit(kNoStats);
   if (!stats) SetBit(kNoStats);
}

//______________________________________________________________________________
void TH1::Sumw2()
{
//*-*-*-*-*Create structure to store sum of squares of weights*-*-*-*-*-*-*-*
//*-*      ===================================================
//*-*
//*-*  if histogram is already filled, the sum of squares of weights
//*-*  is filled with the existing bin contents
//*-*
//*-*  The error per bin will be computed as sqrt(sum of squares of weight)
//*-*  for each bin.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

  if (fSumw2.fN) {
     Warning("Sumw2","Sum of squares of weights structure already created");
     return;
  }

  fSumw2.Set(fNcells);

  for (Int_t bin=0; bin<fNcells; bin++) {
     fSumw2.fArray[bin] = GetBinContent(bin);
  }
}

//______________________________________________________________________________
TF1 *TH1::GetFunction(const char *name)
{
//*-*-*-*-*Return pointer to function with name*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*      ===================================
//
// Functions such as TH1::Fit store the fitted function in the list of
// functions of this histogram.

   return (TF1*)fFunctions->FindObject(name);
}

//______________________________________________________________________________
Stat_t TH1::GetBinError(Int_t bin)
{
//*-*-*-*-*-*-*Return value of error associated to bin number bin*-*-*-*-*
//*-*          ==================================================
//*-*
//*-* if the sum of squares of weights has been defined (via Sumw2),
//*-* this function returns the sqrt(sum of w2).
//*-* otherwise it returns the sqrt(contents) for this bin.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

  if (bin < 0) bin = 0;
  if (bin >= fNcells) bin = fNcells-1;
  if (fSumw2.fN) return TMath::Sqrt(fSumw2.fArray[bin]);
  Stat_t error2 = TMath::Abs(GetBinContent(bin));
  return TMath::Sqrt(error2);
}

//______________________________________________________________________________
Stat_t TH1::GetCellContent(Int_t binx, Int_t biny)
{
  Stat_t binval;
  if (GetDimension() > 1) {
     binval = GetBinContent(biny*(fXaxis.GetNbins()+2) + binx);
     return binval;
  }
  return GetBinContent(binx);
}

//______________________________________________________________________________
Stat_t TH1::GetCellError(Int_t binx, Int_t biny)
{
  if (GetDimension() > 1) return GetBinError(biny*(fXaxis.GetNbins()+2) + binx);
  return GetBinError(binx);
}

//______________________________________________________________________________
void TH1::SetBinError(Int_t bin, Stat_t error)
{
  if (!fSumw2.fN) Sumw2();
  if (bin <0 || bin>= fSumw2.fN) return;
  fSumw2.fArray[bin] = error*error;
}

//______________________________________________________________________________
void TH1::SetCellContent(Int_t binx, Int_t biny, Stat_t content)
{
  if (binx <0 || binx>fXaxis.GetNbins()+1) return;
  if (biny <0 || biny>fYaxis.GetNbins()+1) return;
  SetBinContent(biny*(fXaxis.GetNbins()+2) + binx,content);
}

//______________________________________________________________________________
void TH1::SetCellError(Int_t binx, Int_t biny, Stat_t error)
{
  if (binx <0 || binx>fXaxis.GetNbins()+1) return;
  if (biny <0 || biny>fYaxis.GetNbins()+1) return;
  if (!fSumw2.fN) Sumw2();
  Int_t bin = biny*(fXaxis.GetNbins()+2) + binx;
  fSumw2.fArray[bin] = error*error;
}

//______________________________________________________________________________
void TH1::SetBinContent(Int_t, Stat_t)
{
   AbstractMethod("SetBinContent");
}

ClassImp(TH1C)

//______________________________________________________________________________
//                     TH1C methods
//______________________________________________________________________________
TH1C::TH1C(): TH1(), TArrayC()
{
   fDimension = 1;
}

//______________________________________________________________________________
TH1C::TH1C(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayC(nbins+2)
{
//
//    Create a 1-Dim histogram with fix bins of type char (one byte per channel)
//    ==========================================================================
//                    (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1C::TH1C(const char *name,const char *title,Int_t nbins,Float_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayC(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type char (one byte per channel)
//    ==========================================================================
//                    (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1C::TH1C(const char *name,const char *title,Int_t nbins,Double_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayC(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type char (one byte per channel)
//    ==========================================================================
//                    (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1C::~TH1C()
{

}

//______________________________________________________________________________
TH1C::TH1C(const TH1C &h1c)
{
   ((TH1C&)h1c).Copy(*this);
}

//______________________________________________________________________________
void TH1C::AddBinContent(Int_t bin)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                ==========================

   if (fArray[bin] < 127) fArray[bin]++;
}

//______________________________________________________________________________
void TH1C::AddBinContent(Int_t bin, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by w*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                ==========================

   Int_t newval = fArray[bin] + Int_t(w);
   if (newval > -128 && newval < 128) {fArray[bin] = Char_t(newval); return;}
   if (newval < -127) fArray[bin] = -127;
   if (newval >  127) fArray[bin] =  127;
}

//______________________________________________________________________________
void TH1C::Copy(TObject &newth1)
{
   TH1::Copy(newth1);
   TArrayC::Copy((TH1C&)newth1);
}

//______________________________________________________________________________
TH1 *TH1C::DrawCopy(Option_t *option)
{

   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TH1C *newth1 = (TH1C*)Clone();
   newth1->SetDirectory(0);
   newth1->SetBit(kCanDelete);
   newth1->AppendPad(opt.Data());
   return newth1;
}

//______________________________________________________________________________
Stat_t TH1C::GetBinContent(Int_t bin)
{
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
   return Stat_t (fArray[bin]);
}

//______________________________________________________________________________
void TH1C::Reset(Option_t *option)
{
   TH1::Reset(option);
   TArrayC::Reset();
}

//______________________________________________________________________________
TH1C& TH1C::operator=(const TH1C &h1)
{
   if (this != &h1)  ((TH1C&)h1).Copy(*this);
   return *this;
}


//______________________________________________________________________________
TH1C operator*(Double_t c1, TH1C &h1)
{
   TH1C hnew = h1;
   hnew.Scale(c1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1C operator+(TH1C &h1, TH1C &h2)
{
   TH1C hnew = h1;
   hnew.Add(&h2,1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1C operator-(TH1C &h1, TH1C &h2)
{
   TH1C hnew = h1;
   hnew.Add(&h2,-1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1C operator*(TH1C &h1, TH1C &h2)
{
   TH1C hnew = h1;
   hnew.Multiply(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1C operator/(TH1C &h1, TH1C &h2)
{
   TH1C hnew = h1;
   hnew.Divide(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

ClassImp(TH1S)

//______________________________________________________________________________
//                     TH1S methods
//______________________________________________________________________________
 TH1S::TH1S(): TH1(), TArrayS()
{
   fDimension = 1;
}

//______________________________________________________________________________
 TH1S::TH1S(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayS(nbins+2)
{
//
//    Create a 1-Dim histogram with fix bins of type short
//    ====================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
 TH1S::TH1S(const char *name,const char *title,Int_t nbins,Float_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayS(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type short
//    =========================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
 TH1S::TH1S(const char *name,const char *title,Int_t nbins,Double_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayS(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type short
//    =========================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
 TH1S::~TH1S()
{

}

//______________________________________________________________________________
 TH1S::TH1S(const TH1S &h1s)
{
   ((TH1S&)h1s).Copy(*this);
}

//______________________________________________________________________________
 void TH1S::AddBinContent(Int_t bin)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by 1*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                ==========================

   if (fArray[bin] < 32767) fArray[bin]++;
}

//______________________________________________________________________________
 void TH1S::AddBinContent(Int_t bin, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*Increment bin content by w*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*                ==========================

   Int_t newval = fArray[bin] + Int_t(w);
   if (newval > -32768 && newval < 32768) {fArray[bin] = Short_t(newval); return;}
   if (newval < -32767) fArray[bin] = -32767;
   if (newval >  32767) fArray[bin] =  32767;
}

//______________________________________________________________________________
 void TH1S::Copy(TObject &newth1)
{
   TH1::Copy(newth1);
   TArrayS::Copy((TH1S&)newth1);
}

//______________________________________________________________________________
 TH1 *TH1S::DrawCopy(Option_t *option)
{
   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TH1S *newth1 = (TH1S*)Clone();
   newth1->SetDirectory(0);
   newth1->SetBit(kCanDelete);
   newth1->AppendPad(opt.Data());
   return newth1;
}

//______________________________________________________________________________
 Stat_t TH1S::GetBinContent(Int_t bin)
{
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
   return Stat_t (fArray[bin]);
}

//______________________________________________________________________________
 void TH1S::Reset(Option_t *option)
{
   TH1::Reset(option);
   TArrayS::Reset();
}

//______________________________________________________________________________
TH1S& TH1S::operator=(const TH1S &h1)
{
   if (this != &h1)  ((TH1S&)h1).Copy(*this);
   return *this;
}


//______________________________________________________________________________
TH1S operator*(Double_t c1, TH1S &h1)
{
   TH1S hnew = h1;
   hnew.Scale(c1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1S operator+(TH1S &h1, TH1S &h2)
{
   TH1S hnew = h1;
   hnew.Add(&h2,1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1S operator-(TH1S &h1, TH1S &h2)
{
   TH1S hnew = h1;
   hnew.Add(&h2,-1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1S operator*(TH1S &h1, TH1S &h2)
{
   TH1S hnew = h1;
   hnew.Multiply(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1S operator/(TH1S &h1, TH1S &h2)
{
   TH1S hnew = h1;
   hnew.Divide(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

ClassImp(TH1F)

//______________________________________________________________________________
//                     TH1F methods
//______________________________________________________________________________
TH1F::TH1F(): TH1(), TArrayF()
{
   fDimension = 1;
}

//______________________________________________________________________________
TH1F::TH1F(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayF(nbins+2)
{
//
//    Create a 1-Dim histogram with fix bins of type float
//    ====================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1F::TH1F(const char *name,const char *title,Int_t nbins,Float_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayF(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type float
//    =========================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1F::TH1F(const char *name,const char *title,Int_t nbins,Double_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayF(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type float
//    =========================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1F::TH1F(const TH1F &h)
{
   ((TH1F&)h).Copy(*this);
}

//______________________________________________________________________________
TH1F::~TH1F()
{

}

//______________________________________________________________________________
void TH1F::Copy(TObject &newth1)
{
   TH1::Copy(newth1);
   TArrayF::Copy((TH1F&)newth1);
}

//______________________________________________________________________________
TH1 *TH1F::DrawCopy(Option_t *option)
{
   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TH1F *newth1 = (TH1F*)Clone();
   newth1->SetDirectory(0);
   newth1->SetBit(kCanDelete);
   newth1->AppendPad(opt.Data());
   return newth1;
}

//______________________________________________________________________________
Stat_t TH1F::GetBinContent(Int_t bin)
{
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
   return Stat_t (fArray[bin]);
}

//______________________________________________________________________________
void TH1F::Reset(Option_t *option)
{
   TH1::Reset(option);
   TArrayF::Reset();
}

//______________________________________________________________________________
TH1F& TH1F::operator=(const TH1F &h1)
{
   if (this != &h1)  ((TH1F&)h1).Copy(*this);
   return *this;
}


//______________________________________________________________________________
TH1F operator*(Double_t c1, TH1F &h1)
{
   TH1F hnew = h1;
   hnew.Scale(c1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1F operator+(TH1F &h1, TH1F &h2)
{
   TH1F hnew = h1;
   hnew.Add(&h2,1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1F operator-(TH1F &h1, TH1F &h2)
{
   TH1F hnew = h1;
   hnew.Add(&h2,-1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1F operator*(TH1F &h1, TH1F &h2)
{
   TH1F hnew = h1;
   hnew.Multiply(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1F operator/(TH1F &h1, TH1F &h2)
{
   TH1F hnew = h1;
   hnew.Divide(&h2);
   hnew.SetDirectory(0);
   return hnew;
}


ClassImp(TH1D)

//______________________________________________________________________________
//                     TH1D methods
//______________________________________________________________________________
TH1D::TH1D(): TH1(), TArrayD()
{
   fDimension = 1;
}

//______________________________________________________________________________
TH1D::TH1D(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayD(nbins+2)
{
//
//    Create a 1-Dim histogram with fix bins of type double
//    =====================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1D::TH1D(const char *name,const char *title,Int_t nbins,Float_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayD(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type double
//    =====================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1D::TH1D(const char *name,const char *title,Int_t nbins,Double_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayD(nbins+2)
{
//
//    Create a 1-Dim histogram with variable bins of type double
//    =====================================================
//           (see TH1::TH1 for explanation of parameters)
//
   fDimension = 1;
}

//______________________________________________________________________________
TH1D::~TH1D()
{

}

//______________________________________________________________________________
TH1D::TH1D(const TH1D &h1d)
{
   ((TH1D&)h1d).Copy(*this);
}

//______________________________________________________________________________
void TH1D::Copy(TObject &newth1)
{
   TH1::Copy(newth1);
   TArrayD::Copy((TH1D&)newth1);
}

//______________________________________________________________________________
TH1 *TH1D::DrawCopy(Option_t *option)
{
   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TH1D *newth1 = (TH1D*)Clone();
   newth1->SetDirectory(0);
   newth1->SetBit(kCanDelete);
   newth1->AppendPad(opt.Data());
   return newth1;
}

//______________________________________________________________________________
Stat_t TH1D::GetBinContent(Int_t bin)
{
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
   return Stat_t (fArray[bin]);
}

//______________________________________________________________________________
void TH1D::Reset(Option_t *option)
{
   TH1::Reset(option);
   TArrayD::Reset();
}

//______________________________________________________________________________
TH1D& TH1D::operator=(const TH1D &h1)
{
   if (this != &h1)  ((TH1D&)h1).Copy(*this);
   return *this;
}

//______________________________________________________________________________
TH1D operator*(Double_t c1, TH1D &h1)
{
   TH1D hnew = h1;
   hnew.Scale(c1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1D operator+(TH1D &h1, TH1D &h2)
{
   TH1D hnew = h1;
   hnew.Add(&h2,1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1D operator-(TH1D &h1, TH1D &h2)
{
   TH1D hnew = h1;
   hnew.Add(&h2,-1);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1D operator*(TH1D &h1, TH1D &h2)
{
   TH1D hnew = h1;
   hnew.Multiply(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1D operator/(TH1D &h1, TH1D &h2)
{
   TH1D hnew = h1;
   hnew.Divide(&h2);
   hnew.SetDirectory(0);
   return hnew;
}

//______________________________________________________________________________
TH1 *R__H(Int_t hid)
{
   //return pointer to histogram with name
   //   hid if id >=0
   //   h_id if id <0

   char hname[20];
   if(hid >= 0) sprintf(hname,"h%d",hid);
   else         sprintf(hname,"h_%d",hid);
   return (TH1*)gDirectory->Get(hname);
}

//______________________________________________________________________________
TH1 *R__H(const char * hname)
{
   //return pointer to histogram with name hname

   return (TH1*)gDirectory->Get(hname);
}




ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.