//*CMZ :  2.23/08 31/10/99  18.32.25  by  Rene Brun
//*CMZ :  2.23/07 25/10/99  08.20.53  by  Rene Brun
//*CMZ :  2.23/04 13/10/99  17.58.35  by  Rene Brun
//*CMZ :  2.23/03 27/09/99  14.25.41  by  Rene Brun
//*CMZ :  2.23/02 07/09/99  14.53.38  by  Rene Brun
//*CMZ :  2.23/01 31/08/99  16.59.36  by  Rene Brun
//*CMZ :  2.23/00 30/07/99  20.04.37  by  Fons Rademakers
//*-- Author :    Rene Brun   26/12/94

//*KEEP,CopyRight,T=C.
/*************************************************************************
 * Copyright(c) 1995-1999, The ROOT System, All rights reserved.         *
 * Authors: Rene Brun and Fons Rademakers.                               *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/AA_LICENSE.                      *
 * For the list of contributors see $ROOTSYS/AA_CREDITS.                 *
 *************************************************************************/
//*KEND.

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

//*KEEP,TROOT.
#include "TROOT.h"
//*KEEP,TH1.
#include "TH1.h"
//*KEEP,TH2.
#include "TH2.h"
//*KEEP,TF2.
#include "TF2.h"
//*KEEP,TF3.
#include "TF3.h"
//*KEEP,TVirtualPad.
#include "TVirtualPad.h"
//*KEEP,Foption.
#include "Foption.h"
//*KEEP,TMath.
#include "TMath.h"
//*KEEP,TRandom,T=C++.
#include "TRandom.h"
//*KEEP,TVirtualFitter,T=C++.
#include "TVirtualFitter.h"
//*KEEP,TProfile.
#include "TProfile.h"
//*KEEP,TStyle.
#include "TStyle.h"
//*KEND.

//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*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 class TProfile
//*-*
//*-*- All histogram classes are derived from the base class TH1
//*-*   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.
//*-*
//*-*   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
//*-*
//
/*

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

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 Float_t xmin, xmax, ymin, ymax, zmin, zmax, binwidx, binwidy, binwidz;
const Int_t kCloseDirectory = BIT(7);
const Int_t kNoStats     = BIT(9);
const Int_t kUserContour = BIT(10);
const Int_t kAxisRange   = BIT(11);

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(this);
   fIntegral  = 0;
   fPainter   = 0;
}

//______________________________________________________________________________
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(kCloseDirectory)) fDirectory->GetList()->Remove(this);
   }
   delete fPainter;
   fDirectory = 0;
   fFunctions = 0;
}

//______________________________________________________________________________
TH1::TH1(const Text_t *name,const Text_t *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 Text_t *name,const Text_t *title,Int_t nbins,Axis_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.);

   UseCurrentStyle();

   if (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(this);
}

//______________________________________________________________________________
void TH1::Add(TH1 *h1, Float_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
   Float_t ac1 = TMath::Abs(c1);
   fEntries += ac1*h1->GetEntries();
   fTsumw   += ac1*h1->fTsumw;
   fTsumw2  += ac1*h1->fTsumw2;
   fTsumwx  += ac1*h1->fTsumwx;
   fTsumwx2 += ac1*h1->fTsumwx2;

//*-*- 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, Float_t c1, Float_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
   Float_t ac1 = TMath::Abs(c1);
   Float_t ac2 = TMath::Abs(c2);
   fEntries = ac1*h1->GetEntries() + ac2*h2->GetEntries();
   fTsumw   = ac1*h1->fTsumw       + ac2*h2->fTsumw;
   fTsumw2  = ac1*h1->fTsumw2      + ac2*h2->fTsumw2;
   fTsumwx  = ac1*h1->fTsumwx      + ac2*h2->fTsumwx;
   fTsumwx2 = ac1*h1->fTsumwx2     + ac2*h2->fTsumwx2;

//*-*- 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");
}

//______________________________________________________________________________
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;
   for (binz=1;binz<=nbinsz;binz++) {
      for (biny=1;biny<=nbinsy;biny++) {
         for (binx=1;binx<=nbinsx;binx++) {
            ibin++;
            bin  = GetBin(binx, biny, binz);
            fIntegral[ibin] = fIntegral[ibin-1] + GetBinContent(bin);
         }
      }
   }

//*-*- 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*-*-*-*-*-*-*-*-*-*-*-*
//*-*          =======================================

   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);
   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(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.
//

   if (!h1) {
      Error("Divide","Attempt to divide 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   = fTsumw2 = fTsumwx = fTsumwx2 = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t c0,c1,w,z,x;
   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);
            z = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            fEntries++;
            fTsumw   += z;
            fTsumw2  += z*z;
            fTsumwx  += z*x;
            fTsumwx2 += z*x*x;
            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);
            }
         }
      }
   }
}


//______________________________________________________________________________
void TH1::Divide(TH1 *h1, TH1 *h2, Float_t c1, Float_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.
//

   TString opt = option;
   opt.ToLower();
   Bool_t binomial = kFALSE;
   if (opt.Contains("b")) binomial = kTRUE;
   if (!h1 || !h2) {
      Error("Divide","Attempt to divide a non-existing 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   = fTsumw2 = fTsumwx = fTsumwx2 = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t b1,b2,w,z,x,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);
            z = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            fEntries++;
            fTsumw   += z;
            fTsumw2  += z*z;
            fTsumwx  += z*x;
            fTsumwx2 += z*x*x;
            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);
               }
            }
         }
      }
   }
}

//______________________________________________________________________________
void TH1::Draw(Option_t *option)
{
//*-*-*-*-*-*-*-*-*-*-*Draw this histogram with options*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ================================
//*-*
//*-*  This histogram is added in the list of objects to be drawn in the
//*-*  current pad.
//*-*  It is important to note that only a reference to this histogram is added.
//*-*  This means that if the histogram contents change between pad redraws,
//*-*  the current histogram contents will be shown.
//*-*  Use function DrawCopy to force a copy of the histogram at the time
//*-*  of the Draw operation in the Pad list.
//*-*
//*-*  The following options are supported on all types:
//*-*    "SAME"   : Superimpose on previous picture in the same pad
//*-*    "CYL"    : Use Cylindrical coordinates
//*-*    "POL"    : Use Polar coordinates
//*-*    "SPH"    : Use Spherical coordinates
//*-*    "PSR"    : Use PseudoRapidity/Phi coordinates
//*-*    "LEGO"   : Draw a lego plot with hidden line removal
//*-*    "LEGO1"  : Draw a lego plot with hidden surface removal
//*-*    "LEGO2"  : Draw a lego plot using colors to show the cell contents
//*-*    "SURF"   : Draw a surface plot with hidden line removal
//*-*    "SURF1"  : Draw a surface plot with hidden surface removal
//*-*    "SURF2"  : Draw a surface plot using colors to show the cell contents
//*-*    "SURF3"  : same as SURF with in addition a contour view drawn on the top
//*-*    "SURF4"  : Draw a surface using Gouraud shading
//*-*
//*-*  The following options are supported for 1-D types:
//*-*    "AH"     : Draw histogram, but not the axis labels and tick marks
//*-*    "B"      : Bar chart option
//*-*    "C"      : Draw a smooth Curve througth the histogram bins
//*-*    "E"      : Draw error bars
//*-*    "E0"     : Draw error bars including bins with o contents
//*-*    "E1"     : Draw error bars with perpendicular lines at the edges
//*-*    "E2"     : Draw error bars with rectangles
//*-*    "E3"     : Draw a fill area througth the end points of the vertical error bars
//*-*    "E4"     : Draw a smoothed filled area through the end points of the error bars
//*-*    "L"      : Draw a line througth the bin contents
//*-*    "P"      : Draw current marker at each bin
//*-*    "*H"     : Draw histogram with a * at each bin
//*-*
//*-*
//*-*  The following options are supported for 2-D types:
//*-*    "ARR"    : arrow mode. Shows gradient between adjacent cells
//*-*    "BOX"    : a box is drawn for each cell with surface proportional to contents
//*-*    "COL"    : a box is drawn for each cell with a color scale varying with contents
//*-*    "COLZ"   : same as "COL". In addition the color mapping is also drawn
//*-*    "CONT"   : Draw a contour plot (same as CONT0)
//*-*    "CONT0"  : Draw a contour plot using colors to distinguish contours
//*-*    "CONT1"  : Draw a contour plot using line styles to distinguish contours
//*-*    "CONT2"  : Draw a contour plot using the same line style for all contours
//*-*    "CONT3"  : Draw a contour plot using fill area colors
//*-*    "FB"     : With LEGO or SURFACE, suppress the Front-Box
//*-*    "BB"     : With LEGO or SURFACE, suppress the Back-Box
//*-*    "SCAT"   : Draw a scatter-plot (default)
//*-*
//*-*  Note that most options above can be concatenated, example:
//*-*  h-Draw("e1same");
//*-*  Options are case insensitive
//*-*
//*-*  When using the options "BOX", "COL" or "COLZ", the color palette used
//*-*  is the one defined in the current style (see TStyle::SetPalette)
//*-*
//*-*  When using the "CONT" or "SURF" or "LEGO" options, the number
//*-*  of contour levels can be changed via TH1::SetContour.
//*-*  (default is 20 equidistant levels)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   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(option);
}

//______________________________________________________________________________
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;
   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();

   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 (add) AddBinContent(bin,fu);
            else     SetBinContent(bin,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.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   //authorize calling this function for 2-d histograms. assume weight=1
   if (fDimension == 2) {
      return Fill12(x, w);
   }

   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];
   }
}

//______________________________________________________________________________
Int_t TH1::Fill12(Axis_t x,Axis_t y)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y by 1*-*-*-*-*-*-*-*-*-*
//*-*                  ==================================
//*-*
//*-* if x or/and y is less than the low-edge of the corresponding axis first bin,
//*-*   the Underflow cell is incremented.
//*-* if x or/and y is greater than the upper edge of corresponding axis last bin,
//*-*   the Overflow cell is incremented.
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by 1in the cell corresponding to x,y.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Int_t binx, biny, bin;
   fEntries++;
   binx = fXaxis.FindBin(x);
   biny = fYaxis.FindBin(y);
   bin  = biny*(fXaxis.GetNbins()+2) + binx;
   AddBinContent(bin);
   if (fSumw2.fN) ++fSumw2.fArray[bin];
   if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
   if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   return bin;
}

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

   //authorize calling this function for 3-d histograms. assume weight=1
   if (fDimension == 3) {
      return Fill(x, y, w, 1);
   }

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

//______________________________________________________________________________
Int_t TH1::Fill(Axis_t x, Axis_t y, Axis_t z, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y,z by a weight 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^2 in the cell corresponding to x,y,z.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   Int_t binx, biny, binz, bin;
   fEntries++;
   binx = fXaxis.FindBin(x);
   biny = fYaxis.FindBin(y);
   binz = fZaxis.FindBin(z);
   bin  =  binx + (fXaxis.GetNbins()+2)*(biny + (fYaxis.GetNbins()+2)*binz);
   AddBinContent(bin,w);
   if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
   if (binx == 0 || binx > fXaxis.GetNbins()) return -1;
   if (biny == 0 || biny > fYaxis.GetNbins()) return -1;
   if (binz == 0 || binz > fZaxis.GetNbins()) return -1;
   return bin;
}

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

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

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

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

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

//*-*--------------Start main loop ntimes
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   Bool_t profile = kFALSE;
   if (IsA() == TProfile::Class()) profile = kTRUE;
   for (loop=0;loop<ntimes;loop++) {
      r1 = gRandom->Rndm(loop);
      ibin = TMath::BinarySearch(nbins,&integral[0],r1);
      binz = ibin/nxy;
      biny = (ibin - nxy*binz)/nbinsx;
      binx = 1 + ibin - nbinsx*(biny + nbinsy*binz);
      if (nbinsz) binz++;
      if (nbinsy) biny++;
      x    = GetBinCenter(binx);
      y    = GetBinCenter(biny);
      z    = GetBinCenter(binz);
      if (fDimension == 1) {
         if (profile) Fill(x,y,1.);
         else         Fill(x, 1.);
      }
      else if (fDimension == 2) Fill(x,y, 1.);
      else if (fDimension == 3) Fill(x,y,z, 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,y,z;
   if (fDimension == 1) {
      for (loop=0;loop<ntimes;loop++) {
         x = h->GetRandom();
         Fill(x);
      }
   } else if (fDimension == 2) {
      for (loop=0;loop<ntimes;loop++) {
         h->GetRandom2(x,y);
         Fill(x,y,1.);
      }
   } else {
      for (loop=0;loop<ntimes;loop++) {
         h->GetRandom3(x,y,z);
         Fill(x,y,z,1.);
      }
   }
}


//______________________________________________________________________________
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 Text_t *fname ,Option_t *option ,Option_t *goption, Float_t xxmin, Float_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
//*-*             = "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.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   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?
   Float_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.FindBin(xmin); if (xfirst < 1) xfirst = 1;
      xlast  = fXaxis.FindBin(xmax); if (xlast > fXaxis.GetLast()) xlast = fXaxis.GetLast();
      yfirst = fYaxis.FindBin(ymin); if (yfirst < 1) yfirst = 1;
      ylast  = fYaxis.FindBin(ymax); if (ylast > fYaxis.GetLast()) ylast = fYaxis.GetLast();
      zfirst = fZaxis.FindBin(zmin); if (zfirst < 1) zfirst = 1;
      zlast  = fZaxis.FindBin(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 = 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);
   }

//*-*- Perform minimization
   arglist[0] = hFitter->GetMaxIterations();
   arglist[1] = 1; //epsilon
   hFitter->ExecuteCommand("MIGRAD",arglist,2);
   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) fFunctions->Delete();
      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();
}

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

   //only valid for 2-D histogram
   if (fDimension != 2) {
      Error("FitSlicesX", "Function only valid for 2-D histograms");
      return;
   }

   Int_t nbins  = fYaxis.GetNbins();
   if (binmin < 1) binmin = 1;
   if (binmax > nbins) binmax = nbins;
   if (binmax < binmin) {binmin = 1; binmax = nbins;}

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

   //Create one histogram for each function parameter
   Int_t ipar;
   char name[80], title[80];
   TH1D *hlist[25];
   for (ipar=0;ipar<npar;ipar++) {
      sprintf(name,"%s_%d",GetName(),ipar);
      sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
      hlist[ipar] = new TH1D(name,title, nbins, fYaxis.GetXmin(), fYaxis.GetXmax());
   }
   sprintf(name,"%s_chi2",GetName());
   TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fYaxis.GetXmin(), fYaxis.GetXmax());

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

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

*/
//

   //only valid for 2-D histogram
   if (fDimension != 2) {
      Error("FitSlicesY", "Function only valid for 2-D histograms");
      return;
   }

   Int_t nbins  = fXaxis.GetNbins();
   if (binmin < 1) binmin = 1;
   if (binmax > nbins) binmax = nbins;
   if (binmax < binmin) {binmin = 1; binmax = nbins;}

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

   //Create one histogram for each function parameter
   Int_t ipar;
   char name[80], title[80];
   TH1D *hlist[25];
   for (ipar=0;ipar<npar;ipar++) {
      sprintf(name,"%s_%d",GetName(),ipar);
      sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
      hlist[ipar] = new TH1D(name,title, nbins, fXaxis.GetXmin(), fXaxis.GetXmax());
   }
   sprintf(name,"%s_chi2",GetName());
   TH1D *hchi2 = new TH1D(name,"chisquare", nbins, fXaxis.GetXmin(), fXaxis.GetXmax());

   //Loop on all bins in X, generate a projection along Y
   Int_t bin;
   Int_t nentries;
   for (bin=binmin;bin<=binmax;bin++) {
      TH1D *hpy = ProjectionY("_temp",bin,bin,"e");
      if (hpy == 0) continue;
      nentries = Int_t(hpy->GetEntries());
      if (nentries == 0 || nentries < cut) {delete hpy; continue;}
      f1->SetParameters(parsave);
      hpy->Fit(fname,option);
      Int_t npfits = f1->GetNumberFitPoints();
      if (npfits > npar && npfits >= cut) {
         for (ipar=0;ipar<npar;ipar++) {
            hlist[ipar]->Fill(GetBinCenter(bin),f1->GetParameter(ipar));
            hlist[ipar]->SetBinError(bin,f1->GetParError(ipar));
         }
         hchi2->Fill(GetBinCenter(bin),f1->GetChisquare()/(npfits-npar));
      }
      delete hpy;
   }
   delete [] parsave;
}

//______________________________________________________________________________
void TH1::FitSlicesZ(TF1 *f1, Int_t binminx, Int_t binmaxx, Int_t binminy, Int_t binmaxy, Int_t cut, Option_t *option)
{
// Project slices along Z in case of a 3-D histogram, then fit each slice
// with function f1 and make a 2-d histogram for each fit parameter
// Only cells in the bin range [binminx,binmaxx] and [binminy,binmaxy] are considered.
// if f1=0, a gaussian is assumed
// Before invoking this function, one can set a subrange to be fitted along Z
// via f1->SetRange(zmin,zmax)
// The argument option (default="QNR") can be used to change the fit options.
//     "Q" means Quiet mode
//     "N" means do not show the result of the fit
//     "R" means fit the function in the specified function range
//
//
//  Example: Assume a 3-d histogram h3
//   Root > h3->FitSlicesZ(); produces 4 TH2D histograms
//          with h3_0 containing parameter 0(Constant) for a Gaus fit
//                    of each cell in X,Y projected along Z
//          with h3_1 containing parameter 1(Mean) for a gaus fit
//          with h3_2 containing parameter 2(RMS)  for a gaus fit
//          with h3_chi2 containing the chisquare/number of degrees of freedom for a gaus fit
//
//   Root > h3->Fit(0,15,22,0,0,10);
//          same as above, but only for bins 15 to 22 along X
//          and only for cells in X,Y for which the corresponding projection
//          along Z has more than cut bins filled.

   //only valid for 3-D histogram
   if (fDimension != 3) {
      Error("FitSlicesZ", "Function only valid for 3-D histograms");
      return;
   }

   Int_t nbinsx  = fXaxis.GetNbins();
   Int_t nbinsy  = fYaxis.GetNbins();
   Int_t nbinsz  = fZaxis.GetNbins();
   if (binminx < 1) binminx = 1;
   if (binmaxx > nbinsx) binmaxx = nbinsx;
   if (binmaxx < binminx) {binminx = 1; binmaxx = nbinsx;}
   if (binminy < 1) binminy = 1;
   if (binmaxy > nbinsy) binmaxy = nbinsy;
   if (binmaxy < binminy) {binminy = 1; binmaxy = nbinsy;}

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

   //Create one 2-d histogram for each function parameter
   Int_t ipar;
   char name[80], title[80];
   TH2D *hlist[25];
   for (ipar=0;ipar<npar;ipar++) {
      sprintf(name,"%s_%d",GetName(),ipar);
      sprintf(title,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
      hlist[ipar] = new TH2D(name,title, nbinsx, fXaxis.GetXmin(), fXaxis.GetXmax()
                                       , nbinsy, fYaxis.GetXmin(), fYaxis.GetXmax());
   }
   sprintf(name,"%s_chi2",GetName());
   TH2D *hchi2 = new TH2D(name,"chisquare", nbinsx, fXaxis.GetXmin(), fXaxis.GetXmax()
                                          , nbinsy, fYaxis.GetXmin(), fYaxis.GetXmax());

   //Loop on all cells in X,Y generate a projection along Z
   TH1D *hpz = new TH1D("R_temp","_temp",nbinsz, fZaxis.GetXmin(), fZaxis.GetXmax());
   Int_t bin,binx,biny,binz;
   for (biny=binminy;biny<=binmaxy;biny++) {
      Float_t y = fYaxis.GetBinCenter(biny);
      for (binx=binminx;binx<=binmaxx;binx++) {
         Float_t x = fXaxis.GetBinCenter(binx);
         hpz->Reset();
         Int_t nfill = 0;
         for (binz=1;binz<=nbinsz;binz++) {
            bin = GetBin(binx,biny,binz);
            Float_t w = GetBinContent(bin);
            if (w == 0) continue;
            hpz->Fill(fZaxis.GetBinCenter(binz),w);
            hpz->SetBinError(binz,GetBinError(bin));
            nfill++;
         }
         if (nfill < cut) continue;
         f1->SetParameters(parsave);
         hpz->Fit(fname,option);
         Int_t npfits = f1->GetNumberFitPoints();
         if (npfits > npar && npfits >= cut) {
            for (ipar=0;ipar<npar;ipar++) {
               hlist[ipar]->Fill(x,y,f1->GetParameter(ipar));
               hlist[ipar]->SetCellError(binx,biny,f1->GetParError(ipar));
            }
            hchi2->Fill(x,y,f1->GetChisquare()/(npfits-npar));
         }
      }
   }
   delete [] parsave;
   delete hpz;
}

//______________________________________________________________________________
Text_t *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.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,"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;
   Float_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("binx=%d, x=%g, cu=%g, fu=%g, eu=%g, fsum=%g, f=%g, u1=%g, u2=%g, u3=%gn",binx,x[0],cu,fu,eu,fsum,f,u[0],u[1],u[2]);

         }
      }
   }
   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;
   Float_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
   Float_t xmin = curHist->GetXaxis()->GetXmin();
   Float_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()/Float_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.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    static Double_t zero = 0.;
    static Double_t one = 1.;
    static Int_t idim = 20;

    static Double_t  b[400]	/* was [20][20] */;
    static Int_t i, k, l, ifail;
    static Double_t power;
    static 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]  = Float_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.)
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

    static Double_t xbar, ybar, x2bar;
    static Int_t i, n;
    static Double_t xybar;
    static Float_t fn, xk, yk;
    static 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    = Float_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;
    static Int_t nmjp1, i, j, l;
    static Int_t im1, jp1, nm1, nmi;
    static Double_t s1, s21, s22;
    static 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.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
   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.

   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;
   }
   Float_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]);
}

//______________________________________________________________________________
void TH1::GetRandom2(Axis_t &x, Axis_t &y)
{
// return 2 random numbers along axis x and y distributed according
// the cellcontents of a 2-dim histogram

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

//______________________________________________________________________________
void TH1::GetRandom3(Axis_t &x, Axis_t &y, Axis_t &z)
{
// return 3 random numbers along axis x , y and z distributed according
// the cellcontents of a 3-dim histogram

   if (fDimension != 3) {
      Error("GetRandom","Function only valid for 3-d histograms");
      return;
   }
   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   Int_t nxy    = nbinsx*nbinsy;
   Int_t nbins  = nxy*nbinsz;
   Double_t integral;
   if (fIntegral) {
      if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral();
   } else {
      integral = ComputeIntegral();
      if (integral == 0 || fIntegral == 0) return;
   }
   Float_t r1 = gRandom->Rndm();
   Int_t ibin = TMath::BinarySearch(nbins,&fIntegral[0],r1);
   Int_t binz = ibin/nxy;
   Int_t biny = (ibin - nxy*binz)/nbinsx;
   Int_t binx = ibin - nbinsx*(biny + nbinsy*binz);
   x = fXaxis.GetBinLowEdge(binx+1)
      +fXaxis.GetBinWidth(binx+1)*(fIntegral[ibin+1]-r1)/(fIntegral[ibin+1] - fIntegral[ibin]);
   y = fYaxis.GetBinLowEdge(biny+1) + fYaxis.GetBinWidth(biny+1)*gRandom->Rndm();
   z = fZaxis.GetBinLowEdge(binz+1) + fZaxis.GetBinWidth(binz+1)*gRandom->Rndm();
}

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

//______________________________________________________________________________
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   = fTsumw2 = fTsumwx = fTsumwx2 = 0;

//*-*- Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t c0,c1,w,z,x;
   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);
            z = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            fEntries++;
            fTsumw   += z;
            fTsumw2  += z*z;
            fTsumwx  += z*x;
            fTsumwx2 += z*x*x;
            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);
            }
         }
      }
   }
}


//______________________________________________________________________________
void TH1::Multiply(TH1 *h1, TH1 *h2, Float_t c1, Float_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,z,x,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);
            z = TMath::Abs(w);
            x = fXaxis.GetBinCenter(binx);
            fEntries++;
            fTsumw   += z;
            fTsumw2  += z*z;
            fTsumwx  += z*x;
            fTsumwx2 += z*x*x;
            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);
            }
         }
      }
   }
}

//______________________________________________________________________________
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) fPainter->Paint(option);
}

//______________________________________________________________________________
TProfile *TH1::ProfileX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a profile histogram along X*-*-*-*-*-*
//*-*      ========================================================
//
//   The projection is made from the channels along the Y axis
//   ranging from firstybin to lastybin included.
//

// Check if this is a 2-d histogram
  if (GetDimension() != 2) {
     Error("ProfileX","Not a 2-D histogram");
     return 0;
  }
  TString opt = option;
  opt.ToLower();
  Int_t nx = fXaxis.GetNbins();
  Int_t ny = fYaxis.GetNbins();
  if (firstybin < 0) firstybin = 0;
  if (lastybin > ny) lastybin = ny;

// Create the profile histogram
  char *pname = (char*)name;
  if (strcmp(name,"_pfx") == 0) {
     Int_t nch = strlen(GetName()) + 5;
     pname = new char[nch];
     sprintf(pname,"%s%s",GetName(),name);
  }
  TProfile *h1 = new TProfile(pname,GetTitle(),nx,fXaxis.GetXmin(),fXaxis.GetXmax(),option);
  if (pname != name)  delete [] pname;

// Fill the profile histogram
  Double_t cont;
  Double_t entries = 0;
  for (Int_t binx =0;binx<=nx;binx++) {
     for (Int_t biny=firstybin;biny<=lastybin;biny++) {
        cont =  GetCellContent(binx,biny);
        if (cont) {
           h1->Fill(fXaxis.GetBinCenter(binx),fYaxis.GetBinCenter(biny), cont);
           entries += cont;
        }
     }
  }
  h1->SetEntries(entries);
  return h1;
}

//______________________________________________________________________________
TProfile *TH1::ProfileY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a profile histogram along Y*-*-*-*-*-*
//*-*      ========================================================
//
//   The projection is made from the channels along the X axis
//   ranging from firstxbin to lastxbin included.
//

// Check if this is a 2-d histogram
  if (GetDimension() != 2) {
     Error("ProfileY","Not a 2-D histogram");
     return 0;
  }
  TString opt = option;
  opt.ToLower();
  Int_t nx = fXaxis.GetNbins();
  Int_t ny = fYaxis.GetNbins();
  if (firstxbin < 0) firstxbin = 0;
  if (lastxbin > nx) lastxbin = nx;

// Create the projection histogram
  char *pname = (char*)name;
  if (strcmp(name,"_pfy") == 0) {
     Int_t nch = strlen(GetName()) + 5;
     pname = new char[nch];
     sprintf(pname,"%s%s",GetName(),name);
  }
  TProfile *h1 = new TProfile(pname,GetTitle(),nx,fYaxis.GetXmin(),fYaxis.GetXmax(),option);
  if (pname != name)  delete [] pname;

// Fill the profile histogram
  Double_t cont;
  Double_t entries = 0;
  for (Int_t biny =0;biny<=ny;biny++) {
     for (Int_t binx=firstxbin;binx<=lastxbin;binx++) {
        cont =  GetCellContent(binx,biny);
        if (cont) {
           h1->Fill(fYaxis.GetBinCenter(biny),fXaxis.GetBinCenter(binx), cont);
           entries += cont;
        }
     }
  }
  h1->SetEntries(entries);
  return h1;
}

//______________________________________________________________________________
TH1D *TH1::ProjectionX(const char *name, Int_t firstybin, Int_t lastybin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a 1-D histogram along X*-*-*-*-*-*-*
//*-*      ====================================================
//
//   The projection is always of the type TH1D.
//   The projection is made from the channels along the Y axis
//   ranging from firstybin to lastybin included.
//
//   if option "E" is specified, the errors are computed.
//

// Check if this is a 2-d histogram
  if (GetDimension() != 2) {
     Error("ProjectionX","Not a 2-D histogram");
     return 0;
  }
  TString opt = option;
  opt.ToLower();
  Int_t nx = fXaxis.GetNbins();
  Int_t ny = fYaxis.GetNbins();
  if (firstybin < 0) firstybin = 0;
  if (lastybin > ny) lastybin = ny;

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

// Fill the projected histogram
  Double_t cont,err,err2;
  Double_t entries = 0;
  for (Int_t binx =0;binx<=nx+1;binx++) {
     err2 = 0;
     for (Int_t biny=firstybin;biny<=lastybin;biny++) {
        cont  = GetCellContent(binx,biny);
        err   = GetCellError(binx,biny);
        err2 += err*err;
        if (cont) {
           h1->Fill(fXaxis.GetBinCenter(binx), cont);
           entries += cont;
        }
     }
     if (computeErrors) h1->SetBinError(binx,TMath::Sqrt(err2));
  }
  h1->SetEntries(entries);
  return h1;
}

//______________________________________________________________________________
TH1D *TH1::ProjectionY(const char *name, Int_t firstxbin, Int_t lastxbin, Option_t *option)
{
//*-*-*-*-*Project a 2-D histogram into a 1-D histogram along Y*-*-*-*-*-*-*
//*-*      ====================================================
//
//   The projection is always of the type TH1D.
//   The projection is made from the channels along the X axis
//   ranging from firstxbin to lastxbin included.
//
//   if option "E" is specified, the errors are computed.
//

// Check if this is a 2-d histogram
  if (GetDimension() != 2) {
     Error("ProjectionY","Not a 2-D histogram");
     return 0;
  }
  TString opt = option;
  opt.ToLower();
  Int_t nx = fXaxis.GetNbins();
  Int_t ny = fYaxis.GetNbins();
  if (firstxbin < 0) firstxbin = 0;
  if (lastxbin > nx) lastxbin = nx;

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

// Fill the projected histogram
  Double_t cont,err,err2;
  Double_t entries = 0;
  for (Int_t biny =0;biny<=ny+1;biny++) {
     err2 = 0;
     for (Int_t binx=firstxbin;binx<=lastxbin;binx++) {
        cont  = GetCellContent(binx,biny);
        err   = GetCellError(binx,biny);
        err2 += err*err;
        if (cont) {
           h1->Fill(fYaxis.GetBinCenter(biny), cont);
           entries += cont;
        }
     }
     if (computeErrors) h1->SetBinError(biny,TMath::Sqrt(err2));
  }
  h1->SetEntries(entries);
  return h1;
}

//______________________________________________________________________________
TH1D *TH1::ProjectionZ(const char *name, Int_t ixmin, Int_t ixmax, Int_t iymin, Int_t iymax, Option_t *option)
{
//*-*-*-*-*Project a 3-D histogram into a 1-D histogram along Z*-*-*-*-*-*-*
//*-*      ====================================================
//
//   The projection is always of the type TH1D.
//   The projection is made from the cells along the X axis
//   ranging from ixmin to ixmax and iymin to iymax included.
//
//   if option "E" is specified, the errors are computed.
//
//   code from Paola Collins & Hans Dijkstra

  if (GetDimension() != 3){
     Error("ProjectionZ", "Not a 3-D histogram");
     return 0;
  }

  TString opt = option;
  opt.ToLower();
  Int_t nx = GetNbinsX();
  Int_t ny = GetNbinsY();
  Int_t nz = GetNbinsZ();
  if (ixmin < 0)  ixmin = 0;
  if (ixmax > nx) ixmax = nx;
  if (iymin < 0)  iymin = 0;
  if (iymax > nx) iymax = ny;

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

// Fill the projected histogram
  Float_t cont,e,e1;
  Double_t entries  = 0;
  Double_t newerror = 0;
  for (Int_t ixbin=ixmin;ixbin<=ixmax;ixbin++){
     for (Int_t iybin=iymin;iybin<=iymax;iybin++){
        for (Int_t binz=0;binz<=(nz+1);binz++){
           Int_t bin = GetBin(ixbin,iybin,binz);
           cont = GetBinContent(bin);
           if (computeErrors) {
              e        = GetBinError(bin);
              e1       = h1->GetBinError(binz);
              newerror = TMath::Sqrt(e*e + e1*e1);
           }
           if (cont) {
              h1->Fill(fZaxis.GetBinCenter(binz), cont);
              entries += cont;
           }
           if (computeErrors) h1->SetBinError(binz,newerror);
        }
     }
  }
  h1->SetEntries(entries);
  return h1;
}

//______________________________________________________________________________
TH1 *TH1::Project3D(Option_t *option)
{
   // Project a 3-d histogram into 1 or 2-d histograms depending on the
   // option parameter
   // option may contain a combination of the characters x,y,z,e
   // option = "x" return the x projection into a TH1D histogram
   // option = "y" return the y projection into a TH1D histogram
   // option = "z" return the z projection into a TH1D histogram
   // option = "xy" return the x versus y projection into a TH2D histogram
   // option = "yx" return the y versus x projection into a TH2D histogram
   // option = "xz" return the x versus z projection into a TH2D histogram
   // option = "zx" return the z versus x projection into a TH2D histogram
   // option = "yz" return the y versus z projection into a TH2D histogram
   // option = "zy" return the z versus y projection into a TH2D histogram
   //
   // If option contains the string "e", errors are computed
   //
   // The projection is made for the selected bins only

  if (GetDimension() != 3) {
     Error("Project3D", "Not a 3-D histogram");
     return 0;
  }

  TString opt = option; opt.ToLower();
  Int_t ixmin = fXaxis.GetFirst();
  Int_t ixmax = fXaxis.GetLast();
  Int_t iymin = fYaxis.GetFirst();
  Int_t iymax = fYaxis.GetLast();
  Int_t izmin = fZaxis.GetFirst();
  Int_t izmax = fZaxis.GetLast();
  Int_t nx = ixmax-ixmin+1;
  Int_t ny = iymax-iymin+1;
  Int_t nz = izmax-izmin+1;
  Int_t pcase = 0;
  if (opt.Contains("x"))  pcase = 1;
  if (opt.Contains("y"))  pcase = 2;
  if (opt.Contains("z"))  pcase = 3;
  if (opt.Contains("xy")) pcase = 4;
  if (opt.Contains("yx")) pcase = 5;
  if (opt.Contains("xz")) pcase = 6;
  if (opt.Contains("zx")) pcase = 7;
  if (opt.Contains("yz")) pcase = 8;
  if (opt.Contains("zy")) pcase = 9;

// Create the projection histogram
  TH1 *h = 0;
  Int_t nch = strlen(GetName()) +opt.Length() +2;
  char *name = new char[nch];
  sprintf(name,"%s_%s",GetName(),option);
  nch = strlen(GetTitle()) +opt.Length() +2;
  char *title = new char[nch];
  sprintf(title,"%s_%s",GetTitle(),option);
  switch (pcase) {
     case 1:
        h = new TH1D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
        break;

     case 2:
        h = new TH1D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
        break;

     case 3:
        h = new TH1D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
        break;

     case 4:
        h = new TH2D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax)
                               ,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
        break;

     case 5:
        h = new TH2D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax)
                               ,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
        break;

     case 6:
        h = new TH2D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax)
                               ,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax));
        break;

     case 7:
        h = new TH2D(name,title,nx,fXaxis.GetBinLowEdge(ixmin),fXaxis.GetBinUpEdge(ixmax)
                               ,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
        break;

     case 8:
        h = new TH2D(name,title,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax)
                               ,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax));
        break;

     case 9:
        h = new TH2D(name,title,ny,fYaxis.GetBinLowEdge(iymin),fYaxis.GetBinUpEdge(iymax)
                               ,nz,fZaxis.GetBinLowEdge(izmin),fZaxis.GetBinUpEdge(izmax));
        break;

     }
  delete [] name;
  delete [] title;
  if (h == 0) return h;

  Bool_t computeErrors = kFALSE;
  if (opt.Contains("e")) {h->Sumw2(); computeErrors = kTRUE;}

// Fill the projected histogram
  Float_t cont,e,e1;
  Double_t entries  = 0;
  Double_t newerror = 0;
  for (Int_t ixbin=ixmin;ixbin<=ixmax;ixbin++){
     Int_t ix = ixbin-ixmin+1;
     for (Int_t iybin=iymin;iybin<=iymax;iybin++){
        Int_t iy = iybin-iymin+1;
        for (Int_t izbin=izmin;izbin<=izmax;izbin++){
           Int_t iz = izbin-izmin+1;
           Int_t bin = GetBin(ixbin,iybin,izbin);
           cont = GetBinContent(bin);
           switch (pcase) {
           case 1:
              if (cont) h->Fill(fXaxis.GetBinCenter(ixbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetBinError(ix);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetBinError(ix,newerror);
              }
              break;

           case 2:
              if (cont) h->Fill(fYaxis.GetBinCenter(iybin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetBinError(iy);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetBinError(iy,newerror);
              }
              break;

           case 3:
              if (cont) h->Fill(fZaxis.GetBinCenter(izbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetBinError(iz);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetBinError(iz,newerror);
              }
              break;

           case 4:
              if (cont) h->Fill(fYaxis.GetBinCenter(iybin),fXaxis.GetBinCenter(ixbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(iy,ix);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(iy,ix,newerror);
              }
              break;

           case 5:
              if (cont) h->Fill(fXaxis.GetBinCenter(ixbin),fYaxis.GetBinCenter(iybin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(ix,iy);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(ix,iy,newerror);
              }
              break;

           case 6:
              if (cont) h->Fill(fZaxis.GetBinCenter(izbin),fXaxis.GetBinCenter(ixbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(iz,ix);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(iz,ix,newerror);
              }
              break;

           case 7:
              if (cont) h->Fill(fXaxis.GetBinCenter(ixbin),fZaxis.GetBinCenter(izbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(ix,iz);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(ix,iz,newerror);
              }
              break;

           case 8:
              if (cont) h->Fill(fZaxis.GetBinCenter(izbin),fYaxis.GetBinCenter(iybin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(iz,iy);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(iz,iy,newerror);
              }
              break;

           case 9:
              if (cont) h->Fill(fYaxis.GetBinCenter(iybin),fZaxis.GetBinCenter(izbin), cont);
              if (computeErrors) {
                 e        = GetBinError(bin);
                 e1       = h->GetCellError(iy,iz);
                 newerror = TMath::Sqrt(e*e + e1*e1);
                 h->SetCellError(iy,iz,newerror);
              }
              break;
           }
           if (cont) {
              entries += cont;
           }
        }
     }
  }
  h->SetEntries(entries);
  return h;
}

//______________________________________________________________________________
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
//   NOTE: This function is currently implemented only for 1-D histograms

   Int_t nbins   = fXaxis.GetNbins();
   Float_t xmin  = fXaxis.GetXmin();
   Float_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
   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::Scale(Float_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;
   Float_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*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*              =====================
   if (b.IsReading()) {
      b.ReadVersion();  //Version_t v = b.ReadVersion();
      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->Delete();
      fFunctions->Streamer(b);
      if (!gROOT->ReadingBasket()) {
         fDirectory = gDirectory;
         if (!gDirectory->GetList()->FindObject(this)) gDirectory->Append(this);
      }
   } else {
      b.WriteVersion(TH1::IsA());
      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);
   }
}

//______________________________________________________________________________
void TH1::Print(Option_t *option)
{
//*-*-*-*-*-*-*Print some global quantities for this histogram*-*-*-*-*-*-*-*
//*-*          ===============================================
//
//  If option "all" is given, bin contents and errors are also printed.
//
   printf( "TH1.Print Name= %s, Entries= %d, Total sum= %gn",GetName(),Int_t(fEntries),GetSumOfWeights());
   if (strcasecmp(option,"all")) return;

   Int_t bin, binx, biny, binz;
   Stat_t w,e;
   Float_t x,y,z;
   if (fDimension == 1) {
      for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();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=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
         y = fYaxis.GetBinCenter(biny);
         for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
            bin = GetBin(binx,biny);
            x = fXaxis.GetBinCenter(binx);
            w = 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=fZaxis.GetFirst();binz<=fZaxis.GetLast();binz++) {
         z = fZaxis.GetBinCenter(binz);
         for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
            y = fYaxis.GetBinCenter(biny);
            for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
               bin = GetBin(binx,biny,binz);
               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 *)
{
//*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-*
//*-*            ===========================================

   fTsumw       = 0;
   fTsumw2      = 0;
   fTsumwx      = 0;
   fTsumwx2     = 0;
   fEntries     = 0;

   fSumw2.Reset();
   fFunctions->Delete();
   fContour.Set(0);
   if (fIntegral) {delete [] fIntegral; fIntegral = 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 = 2*axis;
  return stats[ax]/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 = 2*axis;
  x    = stats[ax]/stats[0];
  rms2 = TMath::Abs(stats[ax+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
   // stats[4] = sumwy   if 2-d or 3-d
   // stats[5] = sumwy2
   // stats[6] = sumwz   if 3-d
   // stats[7] = sumwz2
   //
   // 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, biny, binz;
   Stat_t w;
   Float_t x,y,z;
   if (fDimension == 1) {
      if (fTsumw == 0 || fXaxis.TestBit(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;
      }
   }
   if (fDimension == 2) {
      for (bin=0;bin<6;bin++) stats[bin] = 0;
      for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
         y = fYaxis.GetBinCenter(biny);
         for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
            bin = GetBin(binx,biny);
            x = fXaxis.GetBinCenter(binx);
            w = TMath::Abs(GetBinContent(bin));
            stats[0] += w;
            stats[1] += w*w;
            stats[2] += w*x;
            stats[3] += w*x*x;
            stats[4] += w*y;
            stats[5] += w*y*y;
         }
      }
   }
   if (fDimension == 3) {
      for (bin=0;bin<8;bin++) stats[bin] = 0;
      for (binz=fZaxis.GetFirst();binz<=fZaxis.GetLast();binz++) {
         z = fZaxis.GetBinCenter(binz);
         for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
            y = fYaxis.GetBinCenter(biny);
            for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
               bin = GetBin(binx,biny,binz);
               x = fXaxis.GetBinCenter(binx);
               w = TMath::Abs(GetBinContent(bin));
               stats[0] += w;
               stats[1] += w*w;
               stats[2] += w*x;
               stats[3] += w*x*x;
               stats[4] += w*y;
               stats[5] += w*y*y;
               stats[6] += w*z;
               stats[7] += w*z*z;
            }
         }
      }
   }
}

//______________________________________________________________________________
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.

  Stat_t stats[10];
  GetStats(stats);
  return stats[0];
}

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

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (binx1 < 0) binx1 = 0;
   if (binx2 > nbinsx+1) binx2 = nbinsx+1;
   if (binx2 < binx1)    binx2 = nbinsx;
   Stat_t integral = 0;

//*-*- Loop on bins in specified range
   Int_t bin, binx, biny, binz;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=binx1;binx<=binx2;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            integral += GetBinContent(bin);
         }
      }
   }
   return integral;
}

//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2)
{
//Return integral of bin contents in range [binx1,binx2],[biny1,biny2]
// for a 2-D histogram

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (binx1 < 0) binx1 = 0;
   if (binx2 > nbinsx+1) binx2 = nbinsx+1;
   if (binx2 < binx1)    binx2 = nbinsx;
   if (biny1 < 0) biny1 = 0;
   if (biny2 > nbinsy+1) biny2 = nbinsy+1;
   if (biny2 < biny1)    biny2 = nbinsy;
   Stat_t integral = 0;

//*-*- Loop on bins in specified range
   Int_t bin, binx, biny, binz;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=biny1;biny<=biny2;biny++) {
         for (binx=binx1;binx<=binx2;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            integral += GetBinContent(bin);
         }
      }
   }
   return integral;
}

//______________________________________________________________________________
Stat_t TH1::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2)
{
//Return integral of bin contents in range [binx1,binx2],[biny1,biny2],[binz1,binz2]
// for a 3-D histogram

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   if (fDimension < 2) nbinsy = -1;
   if (fDimension < 3) nbinsz = -1;
   if (binx1 < 0) binx1 = 0;
   if (binx2 > nbinsx+1) binx2 = nbinsx+1;
   if (biny1 < 0) biny1 = 0;
   if (biny2 > nbinsy+1) biny2 = nbinsy+1;
   if (binz1 < 0) binz1 = 0;
   if (binz2 > nbinsz+1) binz2 = nbinsz+1;
   Stat_t integral = 0;

//*-*- Loop on bins in specified range
   Int_t bin, binx, biny, binz;
   for (binz=binz1;binz<=binz2;binz++) {
      for (biny=biny1;biny<=biny2;biny++) {
         for (binx=binx1;binx<=binx2;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            integral += GetBinContent(bin);
         }
      }
   }
   return integral;
}

//______________________________________________________________________________
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(Float_t *levels)
{
//*-*-*-*-*-*-*-*Return contour values into array levels*-*-*-*-*-*-*-*-*-*
//*-*            =======================================
//*-*
//*-*  The number of contour levels can be returned by getContourLevel
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  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;
}

//______________________________________________________________________________
Float_t TH1::GetContourLevel(Int_t level)
{
//*-*-*-*-*-*-*-*Return the number of contour levels*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*            ===================================
  if (level <0 || level >= fContour.fN) return 0;
  Float_t zlevel = fContour.fArray[level];
  return zlevel;
}


//______________________________________________________________________________
void TH1::SetContour(Int_t  nlevels, Float_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
     Float_t zmin = GetMinimum();
     Float_t zmax = GetMaximum();
     Float_t dz   = (zmax-zmin)/Float_t(nlevels);
     if (gPad->GetLogz()) {
        if (zmax <= 0) return;
        if (zmin <= 0) zmin = 0.001*zmax;
        zmin = TMath::Log10(zmin);
        zmax = TMath::Log10(zmax);
        dz   = (zmax-zmin)/Float_t(nlevels);
     }
     for (level=0; level<nlevels; level++) {
        fContour.fArray[level] = zmin + dz*Float_t(level);
     }
  }
}

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

//______________________________________________________________________________
Float_t TH1::GetMaximum()
{
//*-*-*-*-*-*-*-*-*-*-*Return maximum value of all bins*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ================================
  Int_t binx, biny;
  Float_t maximum, value;
  if (fMaximum != -1111) return fMaximum;
  if (GetDimension() < 2) {
     maximum = GetBinContent(1);
     for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
        value = GetBinContent(binx);
        if (value > maximum) maximum = value;
     }
  } else {
     maximum = GetCellContent(1,1);
     for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
        for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
           value = GetCellContent(binx, biny);
           if (value > maximum) maximum = value;
        }
     }
  }
  return maximum;
}

//______________________________________________________________________________
Float_t TH1::GetMinimum()
{
//*-*-*-*-*-*-*-*-*-*-*Return minimum value of all bins*-*-*-*-*-*-*-*-*-*-*-*
//*-*                  ================================
  Int_t binx, biny;
  Float_t minimum, value;
  if (fMinimum != -1111) return fMinimum;
  if (GetDimension() < 2) {
     minimum = GetBinContent(1);
     for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
        value = GetBinContent(binx);
        if (value < minimum) minimum = value;
     }
  } else {
     minimum = GetCellContent(1,1);
     for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
        for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
           value = GetCellContent(binx, biny);
           if (value < minimum) minimum = value;
        }
     }
  }
  return minimum;
}

//______________________________________________________________________________
void TH1::SetBins(Int_t nx, Float_t xmin, Float_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, Float_t xmin, Float_t xmax, Int_t ny, Float_t ymin, Float_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, Float_t xmin, Float_t xmax, Int_t ny, Float_t ymin, Float_t ymax, Int_t nz, Float_t zmin, Float_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(Float_t maximum)
{
//*-*-*-*-*-*-*-*-*Set the maximum value for the Y axis*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ====================================
// By default the maximum value is automatically set to the maximum
// bin content plyus a margin of 10 per cent.
//
   fMaximum = maximum;
}


//______________________________________________________________________________
void TH1::SetMinimum(Float_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.
//
   fMinimum = minimum;
}

//______________________________________________________________________________
void TH1::SetDirectory(TDirectory *dir)
{
   // 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 Text_t *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::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
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

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

  fSumw2.Set(fNcells);

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

//______________________________________________________________________________
TF1 *TH1::GetFunction(const Text_t *name)
{
//*-*-*-*-*Return pointer to function with name*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*      ===================================

   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(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayC()
{
   fDimension = dim;
}

//______________________________________________________________________________
 TH1C::TH1C(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayC()
{
   fDimension = dim;
}

//______________________________________________________________________________
 TH1C::TH1C(const Text_t *name,const Text_t *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 Text_t *name,const Text_t *title,Int_t nbins,Axis_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(option);
   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*(Float_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(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayS()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1S::TH1S(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayS()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1S::TH1S(const Text_t *name,const Text_t *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 Text_t *name,const Text_t *title,Int_t nbins,Axis_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(option);
   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*(Float_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(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayF()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1F::TH1F(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayF()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1F::TH1F(const Text_t *name,const Text_t *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 Text_t *name,const Text_t *title,Int_t nbins,Axis_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(option);
   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*(Float_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(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t xlow,Axis_t xup)
     : TH1(name,title,nbins,xlow,xup), TArrayD()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1D::TH1D(Int_t dim,const Text_t *name,const Text_t *title,Int_t nbins,Axis_t *xbins)
     : TH1(name,title,nbins,xbins), TArrayD()
{
   fDimension = dim;
}

//______________________________________________________________________________
TH1D::TH1D(const Text_t *name,const Text_t *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 Text_t *name,const Text_t *title,Int_t nbins,Axis_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(option);
   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*(Float_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;
}


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.