TH1.cxx

Go to the documentation of this file.

// @(#)root/hist:$Id$
// Author: Rene Brun   26/12/94
 
/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
 
#include <array>
#include <cctype>
#include <climits>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <sstream>
 
#include "TROOT.h"
#include "TBuffer.h"
#include "TEnv.h"
#include "TClass.h"
#include "TMath.h"
#include "THashList.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TF2.h"
#include "TF3.h"
#include "TPluginManager.h"
#include "TVirtualPad.h"
#include "TRandom.h"
#include "TVirtualFitter.h"
#include "THLimitsFinder.h"
#include "TProfile.h"
#include "TStyle.h"
#include "TVectorF.h"
#include "TVectorD.h"
#include "TBrowser.h"
#include "TError.h"
#include "TVirtualHistPainter.h"
#include "TVirtualFFT.h"
#include "TVirtualPaveStats.h"
 
#include "HFitInterface.h"
#include "Fit/DataRange.h"
#include "Fit/BinData.h"
#include "Math/GoFTest.h"
#include "Math/MinimizerOptions.h"
#include "Math/QuantFuncMathCore.h"
 
#include "TH1Merger.h"
 
/** \addtogroup Histograms
@{
\class TH1C
\brief 1-D histogram with a byte per channel (see TH1 documentation)
\class TH1S
\brief 1-D histogram with a short per channel (see TH1 documentation)
\class TH1I
\brief 1-D histogram with an int per channel (see TH1 documentation)}
\class TH1F
\brief 1-D histogram with a float per channel (see TH1 documentation)}
\class TH1D
\brief 1-D histogram with a double per channel (see TH1 documentation)}
@}
*/
 
/** \class TH1
    \ingroup Histograms
TH1 is the base class of all histogram classes in %ROOT.
 
It provides the common interface for operations such as binning, filling, drawing, which
will be detailed below.
 
-# [Creating histograms](\ref creating-histograms)
   -  [Labelling axes](\ref labelling-axis)
-# [Binning](\ref binning)
   - [Fix or variable bin size](\ref fix-var)
   - [Convention for numbering bins](\ref convention)
   - [Alphanumeric Bin Labels](\ref alpha)
   - [Histograms with automatic bins](\ref auto-bin)
   - [Rebinning](\ref rebinning)
-# [Filling histograms](\ref filling-histograms)
   - [Associated errors](\ref associated-errors)
   - [Associated functions](\ref associated-functions)
   - [Projections of histograms](\ref prof-hist)
   - [Random Numbers and histograms](\ref random-numbers)
   - [Making a copy of a histogram](\ref making-a-copy)
   - [Normalizing histograms](\ref normalizing)
-# [Drawing histograms](\ref drawing-histograms)
   - [Setting Drawing histogram contour levels (2-D hists only)](\ref cont-level)
   - [Setting histogram graphics attributes](\ref graph-att)
   - [Customising how axes are drawn](\ref axis-drawing)
-# [Fitting histograms](\ref fitting-histograms)
-# [Saving/reading histograms to/from a ROOT file](\ref saving-histograms)
-# [Operations on histograms](\ref operations-on-histograms)
-# [Miscellaneous operations](\ref misc)
 
ROOT supports the following histogram types:
 
  - 1-D histograms:
      - TH1C : histograms with one byte per channel.   Maximum bin content = 127
      - TH1S : histograms with one short per channel.  Maximum bin content = 32767
      - TH1I : histograms with one int per channel.    Maximum bin content = INT_MAX (\ref intmax "*")
      - 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 = 127
      - TH2S : histograms with one short per channel.  Maximum bin content = 32767
      - TH2I : histograms with one int per channel.    Maximum bin content = INT_MAX (\ref intmax "*")
      - 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 = 127
      - TH3S : histograms with one short per channel.  Maximum bin content = 32767
      - TH3I : histograms with one int per channel.    Maximum bin content = INT_MAX (\ref intmax "*")
      - TH3F : histograms with one float per channel.  Maximum precision 7 digits
      - TH3D : histograms with one double per channel. Maximum precision 14 digits
  - Profile histograms: See classes  TProfile, TProfile2D and TProfile3D.
      Profile histograms are used to display the mean value of Y and its standard deviation
      for each bin in X. Profile histograms are in many cases an elegant
      replacement of two-dimensional histograms : the inter-relation of two
      measured quantities X and Y can always be visualized by a two-dimensional
      histogram or scatter-plot; If Y is an unknown (but single-valued)
      approximate function of X, this function is displayed by a profile
      histogram with much better precision than by a scatter-plot.
 
<sup>
\anchor intmax (*) INT_MAX = 2147483647 is the [maximum value for a variable of type int.](https://docs.microsoft.com/en-us/cpp/c-language/cpp-integer-limits)
</sup>
 
The inheritance hierarchy looks as follows:
 
\image html classTH1__inherit__graph_org.svg width=100%
 
\anchor creating-histograms
## Creating histograms
 
Histograms are created by invoking one of the constructors, e.g.
~~~ {.cpp}
       TH1F *h1 = new TH1F("h1", "h1 title", 100, 0, 4.4);
       TH2F *h2 = new TH2F("h2", "h2 title", 40, 0, 4, 30, -3, 3);
~~~
Histograms may also be created by:
 
  -  calling the Clone() function, see below
  -  making a projection from a 2-D or 3-D histogram, see below
  -  reading a histogram from a file
 
 When a histogram is created, a reference to it is automatically added
 to the list of in-memory objects for the current file or directory.
 Then the pointer to this histogram in the current directory can be found
 by its name, doing:
~~~ {.cpp}
       TH1F *h1 = (TH1F*)gDirectory->FindObject(name);
~~~
 
 This default behaviour can be changed by:
~~~ {.cpp}
       h->SetDirectory(nullptr);          // for the current histogram h
       TH1::AddDirectory(kFALSE);   // sets a global switch disabling the referencing
~~~
 When the histogram is deleted, the reference to it is removed from
 the list of objects in memory.
 When a file is closed, all histograms in memory associated with this file
 are automatically deleted.
 
\anchor labelling-axis
### Labelling axes
 
 Axis titles can be specified in the title argument of the constructor.
 They must be separated by ";":
~~~ {.cpp}
        TH1F* h=new TH1F("h", "Histogram title;X Axis;Y Axis", 100, 0, 1);
~~~
 The histogram title and the axis titles can be any TLatex string, and
 are persisted if a histogram is written to a file.
 
 Any title can be omitted:
~~~ {.cpp}
        TH1F* h=new TH1F("h", "Histogram title;;Y Axis", 100, 0, 1);
        TH1F* h=new TH1F("h", ";;Y Axis", 100, 0, 1);
~~~
 The method SetTitle() has the same syntax:
~~~ {.cpp}
        h->SetTitle("Histogram title;Another X title Axis");
~~~
Alternatively, the title of each axis can be set directly:
~~~ {.cpp}
       h->GetXaxis()->SetTitle("X axis title");
       h->GetYaxis()->SetTitle("Y axis title");
~~~
For bin labels see \ref binning.
 
\anchor binning
## Binning
 
\anchor fix-var
### Fix or variable bin size
 
 All histogram types support either fix or variable bin sizes.
 2-D histograms may have fix size bins along X and variable size bins
 along Y or vice-versa. The functions to fill, manipulate, draw or access
 histograms are identical in both cases.
 
 Each histogram always contains 3 axis objects of type TAxis: fXaxis, fYaxis and fZaxis.
 To access the axis parameters, use:
~~~ {.cpp}
        TAxis *xaxis = h->GetXaxis(); etc.
        Double_t binCenter = xaxis->GetBinCenter(bin), etc.
~~~
 See class TAxis for a description of all the access functions.
 The axis range is always stored internally in double precision.
 
\anchor convention
### Convention for numbering bins
 
 For all histogram types: nbins, xlow, xup
~~~ {.cpp}
        bin = 0;       underflow bin
        bin = 1;       first bin with low-edge xlow INCLUDED
        bin = nbins;   last bin with upper-edge xup EXCLUDED
        bin = nbins+1; overflow bin
~~~
 In case of 2-D or 3-D histograms, a "global bin" number is defined.
 For example, assuming a 3-D histogram with (binx, biny, binz), the function
~~~ {.cpp}
        Int_t gbin = h->GetBin(binx, biny, binz);
~~~
 returns a global/linearized gbin number. This global gbin is useful
 to access the bin content/error information independently of the dimension.
 Note that to access the information other than bin content and errors
 one should use the TAxis object directly with e.g.:
~~~ {.cpp}
         Double_t xcenter = h3->GetZaxis()->GetBinCenter(27);
~~~
 returns the center along z of bin number 27 (not the global bin)
 in the 3-D histogram h3.
 
\anchor alpha
### Alphanumeric Bin Labels
 
 By default, a histogram axis is drawn with its numeric bin labels.
 One can specify alphanumeric labels instead with:
 
  - call TAxis::SetBinLabel(bin, label);
    This can always be done before or after filling.
    When the histogram is drawn, bin labels will be automatically drawn.
    See examples labels1.C and labels2.C
  - call to a Fill function with one of the arguments being a string, e.g.
~~~ {.cpp}
           hist1->Fill(somename, weight);
           hist2->Fill(x, somename, weight);
           hist2->Fill(somename, y, weight);
           hist2->Fill(somenamex, somenamey, weight);
~~~
    See examples hlabels1.C and hlabels2.C
  - via TTree::Draw. see for example cernstaff.C
~~~ {.cpp}
           tree.Draw("Nation::Division");
~~~
    where "Nation" and "Division" are two branches of a Tree.
 
When using the options 2 or 3 above, the labels are automatically
 added to the list (THashList) of labels for a given axis.
 By default, an axis is drawn with the order of bins corresponding
 to the filling sequence. It is possible to reorder the axis
 
  - alphabetically
  - by increasing or decreasing values
 
 The reordering can be triggered via the TAxis context menu by selecting
 the menu item "LabelsOption" or by calling directly
 TH1::LabelsOption(option, axis) where
 
  - axis may be "X", "Y" or "Z"
  - option may be:
    - "a" sort by alphabetic order
    - ">" sort by decreasing values
    - "<" sort by increasing values
    - "h" draw labels horizontal
    - "v" draw labels vertical
    - "u" draw labels up (end of label right adjusted)
    - "d" draw labels down (start of label left adjusted)
 
 When using the option 2 above, new labels are added by doubling the current
 number of bins in case one label does not exist yet.
 When the Filling is terminated, it is possible to trim the number
 of bins to match the number of active labels by calling
~~~ {.cpp}
           TH1::LabelsDeflate(axis) with axis = "X", "Y" or "Z"
~~~
 This operation is automatic when using TTree::Draw.
 Once bin labels have been created, they become persistent if the histogram
 is written to a file or when generating the C++ code via SavePrimitive.
 
\anchor auto-bin
### Histograms with automatic bins
 
 When a histogram is created with an axis lower limit greater or equal
 to its upper limit, the SetBuffer is automatically called with an
 argument fBufferSize equal to fgBufferSize (default value=1000).
 fgBufferSize may be reset via the static function TH1::SetDefaultBufferSize.
 The axis limits will be automatically computed when the buffer will
 be full or when the function BufferEmpty is called.
 
\anchor rebinning
### Rebinning
 
 At any time, a histogram can be rebinned via TH1::Rebin. This function
 returns a new histogram with the rebinned contents.
 If bin errors were stored, they are recomputed during the rebinning.
 
 
\anchor filling-histograms
## Filling histograms
 
 A histogram is typically filled with statements like:
~~~ {.cpp}
       h1->Fill(x);
       h1->Fill(x, w); //fill with weight
       h2->Fill(x, y)
       h2->Fill(x, y, w)
       h3->Fill(x, y, z)
       h3->Fill(x, y, z, w)
~~~
 or via one of the Fill functions accepting names described above.
 The Fill functions compute the bin number corresponding to the given
 x, y or z argument and increment this bin by the given weight.
 The Fill functions return the bin number for 1-D histograms or global
 bin number for 2-D and 3-D histograms.
 If TH1::Sumw2 has been called before filling, the sum of squares of
 weights is also stored.
 One can also increment directly a bin number via TH1::AddBinContent
 or replace the existing content via TH1::SetBinContent.
 To access the bin content of a given bin, do:
~~~ {.cpp}
       Double_t binContent = h->GetBinContent(bin);
~~~
 
 By default, the bin number is computed using the current axis ranges.
 If the automatic binning option has been set via
~~~ {.cpp}
       h->SetCanExtend(TH1::kAllAxes);
~~~
 then, the Fill Function will automatically extend the axis range to
 accomodate the new value specified in the Fill argument. The method
 used is to double the bin size until the new value fits in the range,
 merging bins two by two. This automatic binning options is extensively
 used by the TTree::Draw function when histogramming Tree variables
 with an unknown range.
 This automatic binning option is supported for 1-D, 2-D and 3-D histograms.
 
 During filling, some statistics parameters are incremented to compute
 the mean value and Root Mean Square with the maximum precision.
 
 In case of histograms of type TH1C, TH1S, TH2C, TH2S, TH3C, TH3S
 a check is made that the bin contents do not exceed the maximum positive
 capacity (127 or 32767). Histograms of all types may have positive
 or/and negative bin contents.
 
\anchor associated-errors
### Associated errors
 By default, for each bin, the sum of weights is computed at fill time.
 One can also call TH1::Sumw2 to force the storage and computation
 of the sum of the square of weights per bin.
 If Sumw2 has been called, the error per bin is computed as the
 sqrt(sum of squares of weights), otherwise the error is set equal
 to the sqrt(bin content).
 To return the error for a given bin number, do:
~~~ {.cpp}
        Double_t error = h->GetBinError(bin);
~~~
 
\anchor associated-functions
### Associated functions
 One or more object (typically a TF1*) can be added to the list
 of functions (fFunctions) associated to each histogram.
 When TH1::Fit is invoked, the fitted function is added to this list.
 Given a histogram h, one can retrieve an associated function
 with:
~~~ {.cpp}
        TF1 *myfunc = h->GetFunction("myfunc");
~~~
 
 
\anchor operations-on-histograms
## Operations on histograms
 
 Many types of operations are supported on histograms or between histograms
 
  -  Addition of a histogram to the current histogram.
  -  Additions of two histograms with coefficients and storage into the current
     histogram.
  -  Multiplications and Divisions are supported in the same way as additions.
  -  The Add, Divide and Multiply functions also exist to add, divide or multiply
     a histogram by a function.
 
 If a histogram has associated error bars (TH1::Sumw2 has been called),
 the resulting error bars are also computed assuming independent histograms.
 In case of divisions, Binomial errors are also supported.
 One can mark a histogram to be an "average" histogram by setting its bit kIsAverage via
   myhist.SetBit(TH1::kIsAverage);
 When adding (see TH1::Add) average histograms, the histograms are averaged and not summed.
 
 
\anchor prof-hist
### Projections of histograms
 
 One can:
 
  -  make a 1-D projection of a 2-D histogram or Profile
     see functions TH2::ProjectionX,Y, TH2::ProfileX,Y, TProfile::ProjectionX
  -  make a 1-D, 2-D or profile out of a 3-D histogram
     see functions TH3::ProjectionZ, TH3::Project3D.
 
 One can fit these projections via:
~~~ {.cpp}
      TH2::FitSlicesX,Y, TH3::FitSlicesZ.
~~~
 
\anchor random-numbers
### Random Numbers and histograms
 
 TH1::FillRandom can be used to randomly fill a histogram using
 the contents of an existing TF1 function or another
 TH1 histogram (for all dimensions).
 For example, the following two statements create and fill a histogram
 10000 times with a default gaussian distribution of mean 0 and sigma 1:
~~~ {.cpp}
       TH1F h1("h1", "histo from a gaussian", 100, -3, 3);
       h1.FillRandom("gaus", 10000);
~~~
 TH1::GetRandom can be used to return a random number distributed
 according to the contents of a histogram.
 
\anchor making-a-copy
### Making a copy of a histogram
 Like for any other ROOT object derived from TObject, one can use
 the Clone() function. This makes an identical copy of the original
 histogram including all associated errors and functions, e.g.:
~~~ {.cpp}
       TH1F *hnew = (TH1F*)h->Clone("hnew");
~~~
 
\anchor normalizing
### Normalizing histograms
 
 One can scale a histogram such that the bins integral is equal to
 the normalization parameter via TH1::Scale(Double_t norm), where norm
 is the desired normalization divided by the integral of the histogram.
 
 
\anchor drawing-histograms
## Drawing histograms
 
 Histograms are drawn via the THistPainter class. Each histogram has
 a pointer to its own painter (to be usable in a multithreaded program).
 Many drawing options are supported.
 See THistPainter::Paint() for more details.
 
 The same histogram can be drawn with different options in different pads.
 When a histogram drawn in a pad is deleted, the histogram is
 automatically removed from the pad or pads where it was drawn.
 If a histogram is drawn in a pad, then filled again, the new status
 of the histogram will be automatically shown in the pad next time
 the pad is updated. One does not need to redraw the histogram.
 To draw the current version of a histogram in a pad, one can use
~~~ {.cpp}
        h->DrawCopy();
~~~
 This makes a clone (see Clone below) of the histogram. Once the clone
 is drawn, the original histogram may be modified or deleted without
 affecting the aspect of the clone.
 
 One can use TH1::SetMaximum() and TH1::SetMinimum() to force a particular
 value for the maximum or the minimum scale on the plot. (For 1-D
 histograms this means the y-axis, while for 2-D histograms these
 functions affect the z-axis).
 
 TH1::UseCurrentStyle() can be used to change all histogram graphics
 attributes to correspond to the current selected style.
 This function must be called for each histogram.
 In case one reads and draws many histograms from a file, one can force
 the histograms to inherit automatically the current graphics style
 by calling before gROOT->ForceStyle().
 
\anchor cont-level
### Setting Drawing histogram contour levels (2-D hists only)
 
 By default contours are automatically generated at equidistant
 intervals. A default value of 20 levels is used. This can be modified
 via TH1::SetContour() or TH1::SetContourLevel().
 the contours level info is used by the drawing options "cont", "surf",
 and "lego".
 
\anchor graph-att
### Setting histogram graphics attributes
 
 The histogram classes inherit from the attribute classes:
 TAttLine, TAttFill, and TAttMarker.
 See the member functions of these classes for the list of options.
 
\anchor axis-drawing
### Customizing how axes are drawn
 
 Use the functions of TAxis, such as
~~~ {.cpp}
 histogram.GetXaxis()->SetTicks("+");
 histogram.GetYaxis()->SetRangeUser(1., 5.);
~~~
 
\anchor fitting-histograms
## Fitting histograms
 
 Histograms (1-D, 2-D, 3-D and Profiles) can be fitted with a user
 specified function or a pre-defined function via TH1::Fit.
 See TH1::Fit(TF1*, Option_t *, Option_t *, Double_t, Double_t) for the fitting documentation and the possible [fitting options](\ref HFitOpt)
 
 The FitPanel can also be used for fitting an histogram. See the [FitPanel documentation](https://root.cern/manual/fitting/#using-the-fit-panel).
 
\anchor saving-histograms
## Saving/reading histograms to/from a ROOT file
 
 The following statements create a ROOT file and store a histogram
 on the file. Because TH1 derives from TNamed, the key identifier on
 the file is the histogram name:
~~~ {.cpp}
        TFile f("histos.root", "new");
        TH1F h1("hgaus", "histo from a gaussian", 100, -3, 3);
        h1.FillRandom("gaus", 10000);
        h1->Write();
~~~
 To read this histogram in another Root session, do:
~~~ {.cpp}
        TFile f("histos.root");
        TH1F *h = (TH1F*)f.Get("hgaus");
~~~
 One can save all histograms in memory to the file by:
~~~ {.cpp}
        file->Write();
~~~
 
 
\anchor misc
## Miscellaneous operations
 
~~~ {.cpp}
        TH1::KolmogorovTest(): statistical test of compatibility in shape
                             between two histograms
        TH1::Smooth() smooths the bin contents of a 1-d histogram
        TH1::Integral() returns the integral of bin contents in a given bin range
        TH1::GetMean(int axis) returns the mean value along axis
        TH1::GetStdDev(int axis)  returns the sigma distribution along axis
        TH1::GetEntries() returns the number of entries
        TH1::Reset() resets the bin contents and errors of a histogram
~~~
 IMPORTANT NOTE: The returned values for GetMean and GetStdDev depend on how the
 histogram statistics are calculated. By default, if no range has been set, the
 returned values are the (unbinned) ones calculated at fill time. If a range has been
 set, however, the values are calculated using the bins in range; THIS IS TRUE EVEN
 IF THE RANGE INCLUDES ALL BINS--use TAxis::SetRange(0, 0) to unset the range.
 To ensure that the returned values are always those of the binned data stored in the
 histogram, call TH1::ResetStats. See TH1::GetStats.
*/
 
TF1 *gF1=nullptr;  //left for back compatibility (use TVirtualFitter::GetUserFunc instead)
 
Int_t  TH1::fgBufferSize   = 1000;
Bool_t TH1::fgAddDirectory = kTRUE;
Bool_t TH1::fgDefaultSumw2 = kFALSE;
Bool_t TH1::fgStatOverflows= kFALSE;
 
extern void H1InitGaus();
extern void H1InitExpo();
extern void H1InitPolynom();
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);
 
namespace {
 
/// Enumeration specifying inconsistencies between two histograms,
/// in increasing severity.
enum EInconsistencyBits {
   kFullyConsistent = 0,
   kDifferentLabels = BIT(0),
   kDifferentBinLimits = BIT(1),
   kDifferentAxisLimits = BIT(2),
   kDifferentNumberOfBins = BIT(3),
   kDifferentDimensions = BIT(4)
};
 
} // namespace
 
ClassImp(TH1);
 
////////////////////////////////////////////////////////////////////////////////
/// Histogram default constructor.
 
TH1::TH1(): TNamed(), TAttLine(), TAttFill(), TAttMarker()
{
   fDirectory     = nullptr;
   fFunctions     = new TList;
   fNcells        = 0;
   fIntegral      = nullptr;
   fPainter       = nullptr;
   fEntries       = 0;
   fNormFactor    = 0;
   fTsumw         = fTsumw2=fTsumwx=fTsumwx2=0;
   fMaximum       = -1111;
   fMinimum       = -1111;
   fBufferSize    = 0;
   fBuffer        = nullptr;
   fBinStatErrOpt = kNormal;
   fStatOverflows = EStatOverflows::kNeutral;
   fXaxis.SetName("xaxis");
   fYaxis.SetName("yaxis");
   fZaxis.SetName("zaxis");
   fXaxis.SetParent(this);
   fYaxis.SetParent(this);
   fZaxis.SetParent(this);
   UseCurrentStyle();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Histogram default destructor.
 
TH1::~TH1()
{
   if (ROOT::Detail::HasBeenDeleted(this)) {
      return;
   }
   delete[] fIntegral;
   fIntegral = nullptr;
   delete[] fBuffer;
   fBuffer = nullptr;
   if (fFunctions) {
      R__WRITE_LOCKGUARD(ROOT::gCoreMutex);
 
      fFunctions->SetBit(kInvalidObject);
      TObject* obj = nullptr;
      //special logic to support the case where the same object is
      //added multiple times in fFunctions.
      //This case happens when the same object is added with different
      //drawing modes
      //In the loop below we must be careful with objects (eg TCutG) that may
      // have been added to the list of functions of several histograms
      //and may have been already deleted.
      while ((obj  = fFunctions->First())) {
         while(fFunctions->Remove(obj)) { }
         if (ROOT::Detail::HasBeenDeleted(obj)) {
            break;
         }
         delete obj;
         obj = nullptr;
      }
      delete fFunctions;
      fFunctions = nullptr;
   }
   if (fDirectory) {
      fDirectory->Remove(this);
      fDirectory = nullptr;
   }
   delete fPainter;
   fPainter = nullptr;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor for fix bin size histograms.
/// Creates the main histogram structure.
///
/// \param[in] name name of histogram (avoid blanks)
/// \param[in] title histogram title.
///            If title is of the form `stringt;stringx;stringy;stringz`,
///            the histogram title is set to `stringt`,
///            the x axis title to `stringx`, the y axis title to `stringy`, etc.
/// \param[in] nbins number of bins
/// \param[in] xlow low edge of first bin
/// \param[in] xup upper edge of last bin (not included in last bin)
 
 
TH1::TH1(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
   Build();
   if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
   fXaxis.Set(nbins,xlow,xup);
   fNcells = fXaxis.GetNbins()+2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor for variable bin size histograms using an input array of type float.
/// Creates the main histogram structure.
///
/// \param[in] name name of histogram (avoid blanks)
/// \param[in] title histogram title.
///            If title is of the form `stringt;stringx;stringy;stringz`
///            the histogram title is set to `stringt`,
///            the x axis title to `stringx`, the y axis title to `stringy`, etc.
/// \param[in] nbins number of bins
/// \param[in] xbins array of low-edges for each bin.
///            This is an array of type float and size nbins+1
 
TH1::TH1(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
   Build();
   if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
   if (xbins) fXaxis.Set(nbins,xbins);
   else       fXaxis.Set(nbins,0,1);
   fNcells = fXaxis.GetNbins()+2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor for variable bin size histograms using an input array of type double.
///
/// \param[in] name name of histogram (avoid blanks)
/// \param[in] title histogram title.
///        If title is of the form `stringt;stringx;stringy;stringz`
///        the histogram title is set to `stringt`,
///        the x axis title to `stringx`, the y axis title to `stringy`, etc.
/// \param[in] nbins number of bins
/// \param[in] xbins array of low-edges for each bin.
///            This is an array of type double and size nbins+1
 
TH1::TH1(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
    :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
{
   Build();
   if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
   if (xbins) fXaxis.Set(nbins,xbins);
   else       fXaxis.Set(nbins,0,1);
   fNcells = fXaxis.GetNbins()+2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Static function: cannot be inlined on Windows/NT.
 
Bool_t TH1::AddDirectoryStatus()
{
   return fgAddDirectory;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Browse the Histogram object.
 
void TH1::Browse(TBrowser *b)
{
   Draw(b ? b->GetDrawOption() : "");
   gPad->Update();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Creates histogram basic data structure.
 
void TH1::Build()
{
   fDirectory     = nullptr;
   fPainter       = nullptr;
   fIntegral      = nullptr;
   fEntries       = 0;
   fNormFactor    = 0;
   fTsumw         = fTsumw2=fTsumwx=fTsumwx2=0;
   fMaximum       = -1111;
   fMinimum       = -1111;
   fBufferSize    = 0;
   fBuffer        = nullptr;
   fBinStatErrOpt = kNormal;
   fStatOverflows = EStatOverflows::kNeutral;
   fXaxis.SetName("xaxis");
   fYaxis.SetName("yaxis");
   fZaxis.SetName("zaxis");
   fYaxis.Set(1,0.,1.);
   fZaxis.Set(1,0.,1.);
   fXaxis.SetParent(this);
   fYaxis.SetParent(this);
   fZaxis.SetParent(this);
 
   SetTitle(fTitle.Data());
 
   fFunctions = new TList;
 
   UseCurrentStyle();
 
   if (TH1::AddDirectoryStatus()) {
      fDirectory = gDirectory;
      if (fDirectory) {
         fFunctions->UseRWLock();
         fDirectory->Append(this,kTRUE);
      }
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Performs the operation: `this = this + c1*f1`
/// if errors are defined (see TH1::Sumw2), errors are also recalculated.
///
/// By default, the function is computed at the centre of the bin.
/// if option "I" is specified (1-d histogram only), the integral of the
/// function in each bin is used instead of the value of the function at
/// the centre of the bin.
///
/// Only bins inside the function range are recomputed.
///
/// IMPORTANT NOTE: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Add
///
/// The function return kFALSE if the Add operation failed
 
Bool_t TH1::Add(TF1 *f1, Double_t c1, Option_t *option)
{
   if (!f1) {
      Error("Add","Attempt to add a non-existing function");
      return kFALSE;
   }
 
   TString opt = option;
   opt.ToLower();
   Bool_t integral = kFALSE;
   if (opt.Contains("i") && fDimension == 1) integral = kTRUE;
 
   Int_t ncellsx = GetNbinsX() + 2; // cells = normal bins + underflow bin + overflow bin
   Int_t ncellsy = GetNbinsY() + 2;
   Int_t ncellsz = GetNbinsZ() + 2;
   if (fDimension < 2) ncellsy = 1;
   if (fDimension < 3) ncellsz = 1;
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   //   - Add statistics
   Double_t s1[10];
   for (Int_t i = 0; i < 10; ++i) s1[i] = 0;
   PutStats(s1);
   SetMinimum();
   SetMaximum();
 
   //   - Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu=0;
   Double_t xx[3];
   Double_t *params = nullptr;
   f1->InitArgs(xx,params);
   for (binz = 0; binz < ncellsz; ++binz) {
      xx[2] = fZaxis.GetBinCenter(binz);
      for (biny = 0; biny < ncellsy; ++biny) {
         xx[1] = fYaxis.GetBinCenter(biny);
         for (binx = 0; binx < ncellsx; ++binx) {
            xx[0] = fXaxis.GetBinCenter(binx);
            if (!f1->IsInside(xx)) continue;
            TF1::RejectPoint(kFALSE);
            bin = binx + ncellsx * (biny + ncellsy * binz);
            if (integral) {
               cu = c1*f1->Integral(fXaxis.GetBinLowEdge(binx), fXaxis.GetBinUpEdge(binx), 0.) / fXaxis.GetBinWidth(binx);
            } else {
               cu  = c1*f1->EvalPar(xx);
            }
            if (TF1::RejectedPoint()) continue;
            AddBinContent(bin,cu);
         }
      }
   }
 
   return kTRUE;
}
 
int TH1::LoggedInconsistency(const char *name, const TH1 *h1, const TH1 *h2, bool useMerge) const
{
   const auto inconsistency = CheckConsistency(h1, h2);
 
   if (inconsistency & kDifferentDimensions) {
      if (useMerge)
         Info(name, "Histograms have different dimensions - trying to use TH1::Merge");
      else {
         Error(name, "Histograms have different dimensions");
      }
   } else if (inconsistency & kDifferentNumberOfBins) {
      if (useMerge)
         Info(name, "Histograms have different number of bins - trying to use TH1::Merge");
      else {
         Error(name, "Histograms have different number of bins");
      }
   } else if (inconsistency & kDifferentAxisLimits) {
      if (useMerge)
         Info(name, "Histograms have different axis limits - trying to use TH1::Merge");
      else
         Warning(name, "Histograms have different axis limits");
   } else if (inconsistency & kDifferentBinLimits) {
      if (useMerge)
         Info(name, "Histograms have different bin limits - trying to use TH1::Merge");
      else
         Warning(name, "Histograms have different bin limits");
   } else if (inconsistency & kDifferentLabels) {
      // in case of different labels -
      if (useMerge)
         Info(name, "Histograms have different labels - trying to use TH1::Merge");
      else
         Info(name, "Histograms have different labels");
   }
 
   return inconsistency;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
///
/// Note also that adding histogram with labels is not supported, histogram will be
/// added merging them by bin number independently of the labels.
/// For adding histogram with labels one should use TH1::Merge
///
/// SPECIAL CASE (Average/Efficiency histograms)
/// For histograms representing averages or efficiencies, one should compute the average
/// of the two histograms and not the sum. One can mark a histogram to be an average
/// histogram by setting its bit kIsAverage with
/// myhist.SetBit(TH1::kIsAverage);
/// Note that the two histograms must have their kIsAverage bit set
///
/// IMPORTANT NOTE1: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Add
///
/// IMPORTANT NOTE2: if h1 has a normalisation factor, the normalisation factor
/// is used , ie  this = this + c1*factor*h1
/// Use the other TH1::Add function if you do not want this feature
///
/// IMPORTANT NOTE3: You should be careful about the statistics of the
/// returned histogram, whose statistics may be binned or unbinned,
/// depending on whether c1 is negative, whether TAxis::kAxisRange is true,
/// and whether TH1::ResetStats has been called on either this or h1.
/// See TH1::GetStats.
///
/// The function return kFALSE if the Add operation failed
 
Bool_t TH1::Add(const TH1 *h1, Double_t c1)
{
   if (!h1) {
      Error("Add","Attempt to add a non-existing histogram");
      return kFALSE;
   }
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   bool useMerge = false;
   const bool considerMerge = (c1 == 1. &&  !this->TestBit(kIsAverage) && !h1->TestBit(kIsAverage) );
   const auto inconsistency = LoggedInconsistency("Add", this, h1, considerMerge);
   // If there is a bad inconsistency and we can't even consider merging, just give up
   if(inconsistency >= kDifferentNumberOfBins && !considerMerge) {
       return false;
   }
   // If there is an inconsistency, we try to use merging
   if(inconsistency > kFullyConsistent) {
      useMerge = considerMerge;
   }
 
   if (useMerge) {
      TList l;
      l.Add(const_cast<TH1*>(h1));
      auto iret = Merge(&l);
      return (iret >= 0);
   }
 
   //    Create Sumw2 if h1 has Sumw2 set
   if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
 
   //   - Add statistics
   Double_t entries = TMath::Abs( GetEntries() + c1 * h1->GetEntries() );
 
   // statistics can be preserved only in case of positive coefficients
   // otherwise with negative c1 (histogram subtraction) one risks to get negative variances
   Bool_t resetStats = (c1 < 0);
   Double_t s1[kNstat] = {0};
   Double_t s2[kNstat] = {0};
   if (!resetStats) {
      // need to initialize to zero s1 and s2 since
      // GetStats fills only used elements depending on dimension and type
      GetStats(s1);
      h1->GetStats(s2);
   }
 
   SetMinimum();
   SetMaximum();
 
   //   - Loop on bins (including underflows/overflows)
   Double_t factor = 1;
   if (h1->GetNormFactor() != 0) factor = h1->GetNormFactor()/h1->GetSumOfWeights();;
   Double_t c1sq = c1 * c1;
   Double_t factsq = factor * factor;
 
   for (Int_t bin = 0; bin < fNcells; ++bin) {
      //special case where histograms have the kIsAverage bit set
      if (this->TestBit(kIsAverage) && h1->TestBit(kIsAverage)) {
         Double_t y1 = h1->RetrieveBinContent(bin);
         Double_t y2 = this->RetrieveBinContent(bin);
         Double_t e1sq = h1->GetBinErrorSqUnchecked(bin);
         Double_t e2sq = this->GetBinErrorSqUnchecked(bin);
         Double_t w1 = 1., w2 = 1.;
 
         // consider all special cases  when bin errors are zero
         // see http://root-forum.cern.ch/viewtopic.php?f=3&t=13299
         if (e1sq) w1 = 1. / e1sq;
         else if (h1->fSumw2.fN) {
            w1 = 1.E200; // use an arbitrary huge value
            if (y1 == 0) {
               // use an estimated error from the global histogram scale
               double sf = (s2[0] != 0) ? s2[1]/s2[0] : 1;
               w1 = 1./(sf*sf);
            }
         }
         if (e2sq) w2 = 1. / e2sq;
         else if (fSumw2.fN) {
            w2 = 1.E200; // use an arbitrary huge value
            if (y2 == 0) {
               // use an estimated error from the global histogram scale
               double sf = (s1[0] != 0) ? s1[1]/s1[0] : 1;
               w2 = 1./(sf*sf);
            }
         }
 
         double y =  (w1*y1 + w2*y2)/(w1 + w2);
         UpdateBinContent(bin, y);
         if (fSumw2.fN) {
            double err2 =  1./(w1 + w2);
            if (err2 < 1.E-200) err2 = 0;  // to remove arbitrary value when e1=0 AND e2=0
            fSumw2.fArray[bin] = err2;
         }
      } else { // normal case of addition between histograms
         AddBinContent(bin, c1 * factor * h1->RetrieveBinContent(bin));
         if (fSumw2.fN) fSumw2.fArray[bin] += c1sq * factsq * h1->GetBinErrorSqUnchecked(bin);
      }
   }
 
   // update statistics (do here to avoid changes by SetBinContent)
   if (resetStats)  {
      // statistics need to be reset in case coefficient are negative
      ResetStats();
   }
   else {
      for (Int_t i=0;i<kNstat;i++) {
         if (i == 1) s1[i] += c1*c1*s2[i];
         else        s1[i] += c1*s2[i];
      }
      PutStats(s1);
      SetEntries(entries);
   }
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
///
/// Note also that adding histogram with labels is not supported, histogram will be
/// added merging them by bin number independently of the labels.
/// For adding histogram ith labels one should use TH1::Merge
///
/// SPECIAL CASE (Average/Efficiency histograms)
/// For histograms representing averages or efficiencies, one should compute the average
/// of the two histograms and not the sum. One can mark a histogram to be an average
/// histogram by setting its bit kIsAverage with
/// myhist.SetBit(TH1::kIsAverage);
/// Note that the two histograms must have their kIsAverage bit set
///
/// IMPORTANT NOTE: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Add
///
/// IMPORTANT NOTE2: You should be careful about the statistics of the
/// returned histogram, whose statistics may be binned or unbinned,
/// depending on whether c1 is negative, whether TAxis::kAxisRange is true,
/// and whether TH1::ResetStats has been called on either this or h1.
/// See TH1::GetStats.
///
/// ANOTHER SPECIAL CASE : h1 = h2 and c2 < 0
/// do a scaling   this = c1 * h1 / (bin Volume)
///
/// The function returns kFALSE if the Add operation failed
 
Bool_t TH1::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
{
 
   if (!h1 || !h2) {
      Error("Add","Attempt to add a non-existing histogram");
      return kFALSE;
   }
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   Bool_t normWidth = kFALSE;
   if (h1 == h2 && c2 < 0) {c2 = 0; normWidth = kTRUE;}
 
   if (h1 != h2) {
      bool useMerge = false;
      const bool considerMerge = (c1 == 1. && c2 == 1. &&  !this->TestBit(kIsAverage) && !h1->TestBit(kIsAverage) );
 
      // We can combine inconsistencies like this, since they are ordered and a
      // higher inconsistency is worse
      auto const inconsistency = std::max(LoggedInconsistency("Add", this, h1, considerMerge),
                                          LoggedInconsistency("Add", h1, h2, considerMerge));
 
      // If there is a bad inconsistency and we can't even consider merging, just give up
      if(inconsistency >= kDifferentNumberOfBins && !considerMerge) {
          return false;
      }
      // If there is an inconsistency, we try to use merging
      if(inconsistency > kFullyConsistent) {
         useMerge = considerMerge;
      }
 
      if (useMerge) {
         TList l;
         // why TList takes non-const pointers ????
         l.Add(const_cast<TH1*>(h1));
         l.Add(const_cast<TH1*>(h2));
         Reset("ICE");
         auto iret = Merge(&l);
         return (iret >= 0);
      }
   }
 
   //    Create Sumw2 if h1 or h2 have Sumw2 set
   if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
 
   //   - Add statistics
   Double_t nEntries = TMath::Abs( c1*h1->GetEntries() + c2*h2->GetEntries() );
 
   // TODO remove
   // statistics can be preserved only in case of positive coefficients
   // otherwise with negative c1 (histogram subtraction) one risks to get negative variances
   // also in case of scaling with the width we cannot preserve the statistics
   Double_t s1[kNstat] = {0};
   Double_t s2[kNstat] = {0};
   Double_t s3[kNstat];
 
 
   Bool_t resetStats = (c1*c2 < 0) || normWidth;
   if (!resetStats) {
      // need to initialize to zero s1 and s2 since
      // GetStats fills only used elements depending on dimension and type
      h1->GetStats(s1);
      h2->GetStats(s2);
      for (Int_t i=0;i<kNstat;i++) {
         if (i == 1) s3[i] = c1*c1*s1[i] + c2*c2*s2[i];
         //else        s3[i] = TMath::Abs(c1)*s1[i] + TMath::Abs(c2)*s2[i];
         else        s3[i] = c1*s1[i] + c2*s2[i];
      }
   }
 
   SetMinimum();
   SetMaximum();
 
   if (normWidth) { // DEPRECATED CASE: belongs to fitting / drawing modules
 
      Int_t nbinsx = GetNbinsX() + 2; // normal bins + underflow, overflow
      Int_t nbinsy = GetNbinsY() + 2;
      Int_t nbinsz = GetNbinsZ() + 2;
 
      if (fDimension < 2) nbinsy = 1;
      if (fDimension < 3) nbinsz = 1;
 
      Int_t bin, binx, biny, binz;
      for (binz = 0; binz < nbinsz; ++binz) {
         Double_t wz = h1->GetZaxis()->GetBinWidth(binz);
         for (biny = 0; biny < nbinsy; ++biny) {
            Double_t wy = h1->GetYaxis()->GetBinWidth(biny);
            for (binx = 0; binx < nbinsx; ++binx) {
               Double_t wx = h1->GetXaxis()->GetBinWidth(binx);
               bin = GetBin(binx, biny, binz);
               Double_t w = wx*wy*wz;
               UpdateBinContent(bin, c1 * h1->RetrieveBinContent(bin) / w);
               if (fSumw2.fN) {
                  Double_t e1 = h1->GetBinError(bin)/w;
                  fSumw2.fArray[bin] = c1*c1*e1*e1;
               }
            }
         }
      }
   } else if (h1->TestBit(kIsAverage) && h2->TestBit(kIsAverage)) {
      for (Int_t i = 0; i < fNcells; ++i) { // loop on cells (bins including underflow / overflow)
         // special case where histograms have the kIsAverage bit set
         Double_t y1 = h1->RetrieveBinContent(i);
         Double_t y2 = h2->RetrieveBinContent(i);
         Double_t e1sq = h1->GetBinErrorSqUnchecked(i);
         Double_t e2sq = h2->GetBinErrorSqUnchecked(i);
         Double_t w1 = 1., w2 = 1.;
 
         // consider all special cases  when bin errors are zero
         // see http://root-forum.cern.ch/viewtopic.php?f=3&t=13299
         if (e1sq) w1 = 1./ e1sq;
         else if (h1->fSumw2.fN) {
            w1 = 1.E200; // use an arbitrary huge value
            if (y1 == 0 ) { // use an estimated error from the global histogram scale
               double sf = (s1[0] != 0) ? s1[1]/s1[0] : 1;
               w1 = 1./(sf*sf);
            }
         }
         if (e2sq) w2 = 1./ e2sq;
         else if (h2->fSumw2.fN) {
            w2 = 1.E200; // use an arbitrary huge value
            if (y2 == 0) { // use an estimated error from the global histogram scale
               double sf = (s2[0] != 0) ? s2[1]/s2[0] : 1;
               w2 = 1./(sf*sf);
            }
         }
 
         double y =  (w1*y1 + w2*y2)/(w1 + w2);
         UpdateBinContent(i, y);
         if (fSumw2.fN) {
            double err2 =  1./(w1 + w2);
            if (err2 < 1.E-200) err2 = 0;  // to remove arbitrary value when e1=0 AND e2=0
            fSumw2.fArray[i] = err2;
         }
      }
   } else { // case of simple histogram addition
      Double_t c1sq = c1 * c1;
      Double_t c2sq = c2 * c2;
      for (Int_t i = 0; i < fNcells; ++i) { // Loop on cells (bins including underflows/overflows)
         UpdateBinContent(i, c1 * h1->RetrieveBinContent(i) + c2 * h2->RetrieveBinContent(i));
         if (fSumw2.fN) {
            fSumw2.fArray[i] = c1sq * h1->GetBinErrorSqUnchecked(i) + c2sq * h2->GetBinErrorSqUnchecked(i);
         }
      }
   }
 
   if (resetStats)  {
      // statistics need to be reset in case coefficient are negative
      ResetStats();
   }
   else {
      // update statistics (do here to avoid changes by SetBinContent)  FIXME remove???
      PutStats(s3);
      SetEntries(nEntries);
   }
 
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Increment bin content by 1.
 
void TH1::AddBinContent(Int_t)
{
   AbstractMethod("AddBinContent");
}
 
////////////////////////////////////////////////////////////////////////////////
/// Increment bin content by a weight w.
 
void TH1::AddBinContent(Int_t, Double_t)
{
   AbstractMethod("AddBinContent");
}
 
////////////////////////////////////////////////////////////////////////////////
/// Sets the flag controlling the automatic add of histograms in memory
///
/// By default (fAddDirectory = kTRUE), histograms are automatically added
/// to the list of objects in memory.
/// Note that one histogram can be removed from its support directory
/// by calling h->SetDirectory(nullptr) or h->SetDirectory(dir) to add it
/// to the list of objects in the directory dir.
///
/// NOTE that this is a static function. To call it, use;
/// TH1::AddDirectory
 
void TH1::AddDirectory(Bool_t add)
{
   fgAddDirectory = add;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Auxiliary function to get the power of 2 next (larger) or previous (smaller)
/// a given x
///
///    next = kTRUE  : next larger
///    next = kFALSE : previous smaller
///
/// Used by the autobin power of 2 algorithm
 
inline Double_t TH1::AutoP2GetPower2(Double_t x, Bool_t next)
{
   Int_t nn;
   Double_t f2 = std::frexp(x, &nn);
   return ((next && x > 0.) || (!next && x <= 0.)) ? std::ldexp(std::copysign(1., f2), nn)
                                                   : std::ldexp(std::copysign(1., f2), --nn);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Auxiliary function to get the next power of 2 integer value larger then n
///
/// Used by the autobin power of 2 algorithm
 
inline Int_t TH1::AutoP2GetBins(Int_t n)
{
   Int_t nn;
   Double_t f2 = std::frexp(n, &nn);
   if (TMath::Abs(f2 - .5) > 0.001)
      return (Int_t)std::ldexp(1., nn);
   return n;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Buffer-based estimate of the histogram range using the power of 2 algorithm.
///
/// Used by the autobin power of 2 algorithm.
///
/// Works on arguments (min and max from fBuffer) and internal inputs: fXmin,
/// fXmax, NBinsX (from fXaxis), ...
/// Result save internally in fXaxis.
///
/// Overloaded by TH2 and TH3.
///
/// Return -1 if internal inputs are inconsistent, 0 otherwise.
 
Int_t TH1::AutoP2FindLimits(Double_t xmi, Double_t xma)
{
   // We need meaningful raw limits
   if (xmi >= xma)
      return -1;
 
   THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this, xmi, xma);
   Double_t xhmi = fXaxis.GetXmin();
   Double_t xhma = fXaxis.GetXmax();
 
   // Now adjust
   if (TMath::Abs(xhma) > TMath::Abs(xhmi)) {
      // Start from the upper limit
      xhma = TH1::AutoP2GetPower2(xhma);
      xhmi = xhma - TH1::AutoP2GetPower2(xhma - xhmi);
   } else {
      // Start from the lower limit
      xhmi = TH1::AutoP2GetPower2(xhmi, kFALSE);
      xhma = xhmi + TH1::AutoP2GetPower2(xhma - xhmi);
   }
 
   // Round the bins to the next power of 2; take into account the possible inflation
   // of the range
   Double_t rr = (xhma - xhmi) / (xma - xmi);
   Int_t nb = TH1::AutoP2GetBins((Int_t)(rr * GetNbinsX()));
 
   // Adjust using the same bin width and offsets
   Double_t bw = (xhma - xhmi) / nb;
   // Bins to left free on each side
   Double_t autoside = gEnv->GetValue("Hist.Binning.Auto.Side", 0.05);
   Int_t nbside = (Int_t)(nb * autoside);
 
   // Side up
   Int_t nbup = (xhma - xma) / bw;
   if (nbup % 2 != 0)
      nbup++; // Must be even
   if (nbup != nbside) {
      // Accounts also for both case: larger or smaller
      xhma -= bw * (nbup - nbside);
      nb -= (nbup - nbside);
   }
 
   // Side low
   Int_t nblw = (xmi - xhmi) / bw;
   if (nblw % 2 != 0)
      nblw++; // Must be even
   if (nblw != nbside) {
      // Accounts also for both case: larger or smaller
      xhmi += bw * (nblw - nbside);
      nb -= (nblw - nbside);
   }
 
   // Set everything and project
   SetBins(nb, xhmi, xhma);
 
   // Done
   return 0;
}
 
/// Fill histogram with all entries in the buffer.
///
///  - action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
///  - action =  0 histogram is reset and filled from the buffer. When the histogram is filled from the
///                buffer the value fBuffer[0] is set to a negative number (= - number of entries)
///                When calling with action == 0 the histogram is NOT refilled when fBuffer[0] is < 0
///                While when calling with action = -1 the histogram is reset and ALWAYS refilled independently if
///                the histogram was filled before. This is needed when drawing the histogram
///  - action =  1 histogram is filled and buffer is deleted
///                The buffer is automatically deleted when filling the histogram and the entries is
///                larger than the buffer size
 
Int_t TH1::BufferEmpty(Int_t action)
{
   // do we need to compute the bin size?
   if (!fBuffer) return 0;
   Int_t nbentries = (Int_t)fBuffer[0];
 
   // nbentries correspond to the number of entries of histogram
 
   if (nbentries == 0) {
      // if action is 1 we delete the buffer
      // this will avoid infinite recursion
      if (action > 0) {
         delete [] fBuffer;
         fBuffer = nullptr;
         fBufferSize = 0;
      }
      return 0;
   }
   if (nbentries < 0 && action == 0) return 0;    // case histogram has been already filled from the buffer
 
   Double_t *buffer = fBuffer;
   if (nbentries < 0) {
      nbentries  = -nbentries;
      //  a reset might call BufferEmpty() giving an infinite recursion
      // Protect it by setting fBuffer = nullptr
      fBuffer = nullptr;
       //do not reset the list of functions
      Reset("ICES");
      fBuffer = buffer;
   }
   if (CanExtendAllAxes() || (fXaxis.GetXmax() <= fXaxis.GetXmin())) {
      //find min, max of entries in buffer
      Double_t xmin = TMath::Infinity();
      Double_t xmax = -TMath::Infinity();
      for (Int_t i=0;i<nbentries;i++) {
         Double_t x = fBuffer[2*i+2];
         // skip infinity or NaN values
         if (!std::isfinite(x)) continue;
         if (x < xmin) xmin = x;
         if (x > xmax) xmax = x;
      }
      if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
         Int_t rc = -1;
         if (TestBit(TH1::kAutoBinPTwo)) {
            if ((rc = AutoP2FindLimits(xmin, xmax)) < 0)
               Warning("BufferEmpty",
                       "inconsistency found by power-of-2 autobin algorithm: fallback to standard method");
         }
         if (rc < 0)
            THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this, xmin, xmax);
      } else {
         fBuffer = nullptr;
         Int_t keep = fBufferSize; fBufferSize = 0;
         if (xmin <  fXaxis.GetXmin()) ExtendAxis(xmin, &fXaxis);
         if (xmax >= fXaxis.GetXmax()) ExtendAxis(xmax, &fXaxis);
         fBuffer = buffer;
         fBufferSize = keep;
      }
   }
 
   // call DoFillN which will not put entries in the buffer as FillN does
   // set fBuffer to zero to avoid re-emptying the buffer from functions called
   // by DoFillN (e.g Sumw2)
   buffer = fBuffer; fBuffer = nullptr;
   DoFillN(nbentries,&buffer[2],&buffer[1],2);
   fBuffer = buffer;
 
   // if action == 1 - delete the buffer
   if (action > 0) {
      delete [] fBuffer;
      fBuffer = nullptr;
      fBufferSize = 0;
   } else {
      // if number of entries is consistent with buffer - set it negative to avoid
      // refilling the histogram every time BufferEmpty(0) is called
      // In case it is not consistent, by setting fBuffer[0]=0 is like resetting the buffer
      // (it will not be used anymore the next time BufferEmpty is called)
      if (nbentries == (Int_t)fEntries)
         fBuffer[0] = -nbentries;
      else
         fBuffer[0] = 0;
   }
   return nbentries;
}
 
////////////////////////////////////////////////////////////////////////////////
/// accumulate arguments in buffer. When buffer is full, empty the buffer
///
///  - `fBuffer[0]` = number of entries in buffer
///  - `fBuffer[1]` = w of first entry
///  - `fBuffer[2]` = x of first entry
 
Int_t TH1::BufferFill(Double_t x, Double_t w)
{
   if (!fBuffer) return -2;
   Int_t nbentries = (Int_t)fBuffer[0];
 
 
   if (nbentries < 0) {
      // reset nbentries to a positive value so next time BufferEmpty()  is called
      // the histogram will be refilled
      nbentries  = -nbentries;
      fBuffer[0] =  nbentries;
      if (fEntries > 0) {
         // set fBuffer to zero to avoid calling BufferEmpty in Reset
         Double_t *buffer = fBuffer; fBuffer=nullptr;
         Reset("ICES");  // do not reset list of functions
         fBuffer = buffer;
      }
   }
   if (2*nbentries+2 >= fBufferSize) {
      BufferEmpty(1);
      if (!fBuffer)
         // to avoid infinite recursion Fill->BufferFill->Fill
         return Fill(x,w);
      // this cannot happen
      R__ASSERT(0);
   }
   fBuffer[2*nbentries+1] = w;
   fBuffer[2*nbentries+2] = x;
   fBuffer[0] += 1;
   return -2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check bin limits.
 
bool TH1::CheckBinLimits(const TAxis* a1, const TAxis * a2)
{
   const TArrayD * h1Array = a1->GetXbins();
   const TArrayD * h2Array = a2->GetXbins();
   Int_t fN = h1Array->fN;
   if ( fN != 0 ) {
      if ( h2Array->fN != fN ) {
         return false;
      }
      else {
         for ( int i = 0; i < fN; ++i ) {
            // for i==fN (nbin+1) a->GetBinWidth() returns last bin width
            // we do not need to exclude that case
            double binWidth = a1->GetBinWidth(i);
            if ( ! TMath::AreEqualAbs( h1Array->GetAt(i), h2Array->GetAt(i), binWidth*1E-10 ) ) {
               return false;
            }
         }
      }
   }
 
   return true;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check that axis have same labels.
 
bool TH1::CheckBinLabels(const TAxis* a1, const TAxis * a2)
{
   THashList *l1 = a1->GetLabels();
   THashList *l2 = a2->GetLabels();
 
   if (!l1 && !l2 )
      return true;
   if (!l1 ||  !l2 ) {
      return false;
   }
   // check now labels sizes  are the same
   if (l1->GetSize() != l2->GetSize() ) {
      return false;
   }
   for (int i = 1; i <= a1->GetNbins(); ++i) {
      TString label1 = a1->GetBinLabel(i);
      TString label2 = a2->GetBinLabel(i);
      if (label1 != label2) {
         return false;
      }
   }
 
   return true;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check that the axis limits of the histograms are the same.
/// If a first and last bin is passed the axis is compared between the given range
 
bool TH1::CheckAxisLimits(const TAxis *a1, const TAxis *a2 )
{
   double firstBin = a1->GetBinWidth(1);
   double lastBin = a1->GetBinWidth( a1->GetNbins() );
   if ( ! TMath::AreEqualAbs(a1->GetXmin(), a2->GetXmin(), firstBin* 1.E-10) ||
        ! TMath::AreEqualAbs(a1->GetXmax(), a2->GetXmax(), lastBin*1.E-10) ) {
      return false;
   }
   return true;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check that the axis are the same
 
bool TH1::CheckEqualAxes(const TAxis *a1, const TAxis *a2 )
{
   if (a1->GetNbins() != a2->GetNbins() ) {
      ::Info("CheckEqualAxes","Axes have different number of bins : nbin1 = %d nbin2 = %d",a1->GetNbins(),a2->GetNbins() );
      return false;
   }
   if(!CheckAxisLimits(a1,a2)) {
      ::Info("CheckEqualAxes","Axes have different limits");
      return false;
   }
   if(!CheckBinLimits(a1,a2)) {
      ::Info("CheckEqualAxes","Axes have different bin limits");
      return false;
   }
 
   // check labels
   if(!CheckBinLabels(a1,a2)) {
      ::Info("CheckEqualAxes","Axes have different labels");
      return false;
   }
 
   return true;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check that two sub axis are the same.
/// The limits are defined by first bin and last bin
/// N.B. no check is done in this case for variable bins
 
bool TH1::CheckConsistentSubAxes(const TAxis *a1, Int_t firstBin1, Int_t lastBin1, const TAxis * a2, Int_t firstBin2, Int_t lastBin2 )
{
   // By default is assumed that no bins are given for the second axis
   Int_t nbins1   = lastBin1-firstBin1 + 1;
   Double_t xmin1 = a1->GetBinLowEdge(firstBin1);
   Double_t xmax1 = a1->GetBinUpEdge(lastBin1);
 
   Int_t nbins2 = a2->GetNbins();
   Double_t xmin2 = a2->GetXmin();
   Double_t xmax2 = a2->GetXmax();
 
   if (firstBin2 <  lastBin2) {
      // in this case assume no bins are given for the second axis
      nbins2   = lastBin1-firstBin1 + 1;
      xmin2 = a1->GetBinLowEdge(firstBin1);
      xmax2 = a1->GetBinUpEdge(lastBin1);
   }
 
   if (nbins1 != nbins2 ) {
      ::Info("CheckConsistentSubAxes","Axes have different number of bins");
      return false;
   }
 
   Double_t firstBin = a1->GetBinWidth(firstBin1);
   Double_t lastBin = a1->GetBinWidth(lastBin1);
   if ( ! TMath::AreEqualAbs(xmin1,xmin2,1.E-10 * firstBin) ||
        ! TMath::AreEqualAbs(xmax1,xmax2,1.E-10 * lastBin) ) {
      ::Info("CheckConsistentSubAxes","Axes have different limits");
      return false;
   }
 
   return true;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check histogram compatibility.
 
int TH1::CheckConsistency(const TH1* h1, const TH1* h2)
{
   if (h1 == h2) return kFullyConsistent;
 
   if (h1->GetDimension() != h2->GetDimension() ) {
      return kDifferentDimensions;
   }
   Int_t dim = h1->GetDimension();
 
   // returns kTRUE if number of bins and bin limits are identical
   Int_t nbinsx = h1->GetNbinsX();
   Int_t nbinsy = h1->GetNbinsY();
   Int_t nbinsz = h1->GetNbinsZ();
 
   // Check whether the histograms have the same number of bins.
   if (nbinsx != h2->GetNbinsX() ||
       (dim > 1 && nbinsy != h2->GetNbinsY())  ||
       (dim > 2 && nbinsz != h2->GetNbinsZ()) ) {
      return kDifferentNumberOfBins;
   }
 
   bool ret = true;
 
   // check axis limits
   ret &= CheckAxisLimits(h1->GetXaxis(), h2->GetXaxis());
   if (dim > 1) ret &= CheckAxisLimits(h1->GetYaxis(), h2->GetYaxis());
   if (dim > 2) ret &= CheckAxisLimits(h1->GetZaxis(), h2->GetZaxis());
   if (!ret) return kDifferentAxisLimits;
 
   // check bin limits
   ret &= CheckBinLimits(h1->GetXaxis(), h2->GetXaxis());
   if (dim > 1) ret &= CheckBinLimits(h1->GetYaxis(), h2->GetYaxis());
   if (dim > 2) ret &= CheckBinLimits(h1->GetZaxis(), h2->GetZaxis());
   if (!ret) return kDifferentBinLimits;
 
   // check labels if histograms are both not empty
   if ( !h1->IsEmpty() && !h2->IsEmpty() ) {
      ret &= CheckBinLabels(h1->GetXaxis(), h2->GetXaxis());
      if (dim > 1) ret &= CheckBinLabels(h1->GetYaxis(), h2->GetYaxis());
      if (dim > 2) ret &= CheckBinLabels(h1->GetZaxis(), h2->GetZaxis());
      if (!ret) return kDifferentLabels;
   }
 
   return kFullyConsistent;
}
 
////////////////////////////////////////////////////////////////////////////////
/// \f$ \chi^{2} \f$ test for comparing weighted and unweighted histograms
///
/// Function: Returns p-value. Other return values are specified by the 3rd parameter
///
/// \param[in] h2 the second histogram
/// \param[in] option
///   - "UU" = experiment experiment comparison (unweighted-unweighted)
///   - "UW" = experiment MC comparison (unweighted-weighted). Note that
///      the first histogram should be unweighted
///   - "WW" = MC MC comparison (weighted-weighted)
///   - "NORM" = to be used when one or both of the histograms is scaled
///              but the histogram originally was unweighted
///   - by default underflows and overflows are not included:
///      * "OF" = overflows included
///      * "UF" = underflows included
///   - "P" = print chi2, ndf, p_value, igood
///   - "CHI2" = returns chi2 instead of p-value
///   - "CHI2/NDF" = returns \f$ \chi^{2} \f$/ndf
/// \param[in] res not empty - computes normalized residuals and returns them in this array
///
/// The current implementation is based on the papers \f$ \chi^{2} \f$ test for comparison
/// of weighted and unweighted histograms" in Proceedings of PHYSTAT05 and
/// "Comparison weighted and unweighted histograms", arXiv:physics/0605123
/// by N.Gagunashvili. This function has been implemented by Daniel Haertl in August 2006.
///
/// #### Introduction:
///
/// A frequently used technique in data analysis is the comparison of
/// histograms. First suggested by Pearson [1] the \f$ \chi^{2} \f$  test of
/// homogeneity is used widely for comparing usual (unweighted) histograms.
/// This paper describes the implementation modified \f$ \chi^{2} \f$ tests
/// for comparison of weighted and unweighted  histograms and two weighted
/// histograms [2] as well as usual Pearson's \f$ \chi^{2} \f$ test for
/// comparison two usual (unweighted) histograms.
///
/// #### Overview:
///
/// Comparison of two histograms expect hypotheses that two histograms
/// represent identical distributions. To make a decision p-value should
/// be calculated. The hypotheses of identity is rejected if the p-value is
/// lower then some significance level. Traditionally significance levels
/// 0.1, 0.05 and 0.01 are used. The comparison procedure should include an
/// analysis of the residuals which is often helpful in identifying the
/// bins of histograms responsible for a significant overall \f$ \chi^{2} \f$ value.
/// Residuals are the difference between bin contents and expected bin
/// contents. Most convenient for analysis are the normalized residuals. If
/// hypotheses of identity are valid then normalized residuals are
/// approximately independent and identically distributed random variables
/// having N(0,1) distribution. Analysis of residuals expect test of above
/// mentioned properties of residuals. Notice that indirectly the analysis
/// of residuals increase the power of \f$ \chi^{2} \f$ test.
///
/// #### Methods of comparison:
///
/// \f$ \chi^{2} \f$ test for comparison two (unweighted) histograms:
/// Let us consider two  histograms with the  same binning and the  number
/// of bins equal to r. Let us denote the number of events in the ith bin
/// in the first histogram as ni and as mi in the second one. The total
/// number of events in the first histogram is equal to:
/// \f[
///  N = \sum_{i=1}^{r} n_{i}
/// \f]
/// and
/// \f[
///  M = \sum_{i=1}^{r} m_{i}
/// \f]
/// in the second histogram. The hypothesis of identity (homogeneity) [3]
/// is that the two histograms represent random values with identical
/// distributions. It is equivalent that there exist r constants p1,...,pr,
/// such that
/// \f[
///\sum_{i=1}^{r} p_{i}=1
/// \f]
/// and the probability of belonging to the ith bin for some measured value
/// in both experiments is equal to pi. The number of events in the ith
/// bin is a random variable with a distribution approximated by a Poisson
/// probability distribution
/// \f[
///\frac{e^{-Np_{i}}(Np_{i})^{n_{i}}}{n_{i}!}
/// \f]
///for the first histogram and with distribution
/// \f[
///\frac{e^{-Mp_{i}}(Mp_{i})^{m_{i}}}{m_{i}!}
/// \f]
/// for the second histogram. If the hypothesis of homogeneity is valid,
/// then the  maximum likelihood estimator of pi, i=1,...,r, is
/// \f[
///\hat{p}_{i}= \frac{n_{i}+m_{i}}{N+M}
/// \f]
/// and then
/// \f[
///  X^{2} = \sum_{i=1}^{r}\frac{(n_{i}-N\hat{p}_{i})^{2}}{N\hat{p}_{i}} + \sum_{i=1}^{r}\frac{(m_{i}-M\hat{p}_{i})^{2}}{M\hat{p}_{i}} =\frac{1}{MN} \sum_{i=1}^{r}\frac{(Mn_{i}-Nm_{i})^{2}}{n_{i}+m_{i}}
/// \f]
/// has approximately a \f$ \chi^{2}_{(r-1)} \f$ distribution [3].
/// The comparison procedure can include an analysis of the residuals which
/// is often helpful in identifying the bins of histograms responsible for
/// a significant overall \f$ \chi^{2} \f$ value. Most convenient for
/// analysis are the adjusted (normalized) residuals [4]
/// \f[
///  r_{i} = \frac{n_{i}-N\hat{p}_{i}}{\sqrt{N\hat{p}_{i}}\sqrt{(1-N/(N+M))(1-(n_{i}+m_{i})/(N+M))}}
/// \f]
/// If hypotheses of  homogeneity are valid then residuals ri are
/// approximately independent and identically distributed random variables
/// having N(0,1) distribution. The application of the \f$ \chi^{2} \f$ test has
/// restrictions related to the value of the expected frequencies Npi,
/// Mpi, i=1,...,r. A conservative rule formulated in [5] is that all the
/// expectations must be 1 or greater for both histograms. In practical
/// cases when expected frequencies are not known the estimated expected
/// frequencies \f$ M\hat{p}_{i}, N\hat{p}_{i}, i=1,...,r \f$ can be used.
///
/// #### Unweighted and weighted histograms comparison:
///
/// A simple modification of the ideas described above can be used for the
/// comparison of the usual (unweighted) and weighted histograms. Let us
/// denote the number of events in the ith bin in the unweighted
/// histogram as ni and the common weight of events in the ith bin of the
/// weighted histogram as wi. The total number of events in the
/// unweighted histogram is equal to
///\f[
///  N = \sum_{i=1}^{r} n_{i}
///\f]
/// and the total weight of events in the weighted histogram is equal to
///\f[
///  W = \sum_{i=1}^{r} w_{i}
///\f]
/// Let us formulate the hypothesis of identity of an unweighted histogram
/// to a weighted histogram so that there exist r constants p1,...,pr, such
/// that
///\f[
///  \sum_{i=1}^{r} p_{i} = 1
///\f]
/// for the unweighted histogram. The weight wi is a random variable with a
/// distribution approximated by the normal probability distribution
/// \f$ N(Wp_{i},\sigma_{i}^{2}) \f$ where \f$ \sigma_{i}^{2} \f$ is the variance of the weight wi.
/// If we replace the variance \f$ \sigma_{i}^{2} \f$
/// with estimate \f$ s_{i}^{2} \f$ (sum of squares of weights of
/// events in the ith bin) and the hypothesis of identity is valid, then the
/// maximum likelihood estimator of  pi,i=1,...,r, is
///\f[
///  \hat{p}_{i} = \frac{Ww_{i}-Ns_{i}^{2}+\sqrt{(Ww_{i}-Ns_{i}^{2})^{2}+4W^{2}s_{i}^{2}n_{i}}}{2W^{2}}
///\f]
/// We may then use the test statistic
///\f[
///  X^{2} = \sum_{i=1}^{r} \frac{(n_{i}-N\hat{p}_{i})^{2}}{N\hat{p}_{i}} + \sum_{i=1}^{r} \frac{(w_{i}-W\hat{p}_{i})^{2}}{s_{i}^{2}}
///\f]
/// and it has approximately a \f$ \sigma^{2}_{(r-1)} \f$ distribution [2]. This test, as well
/// as the original one [3], has a restriction on the expected frequencies. The
/// expected frequencies recommended for the weighted histogram is more than 25.
/// The value of the minimal expected frequency can be decreased down to 10 for
/// the case when the weights of the events are close to constant. In the case
/// of a weighted histogram if the number of events is unknown, then we can
/// apply this recommendation for the equivalent number of events as
///\f[
///  n_{i}^{equiv} = \frac{ w_{i}^{2} }{ s_{i}^{2} }
///\f]
/// The minimal expected frequency for an unweighted histogram must be 1. Notice
/// that any usual (unweighted) histogram can be considered as a weighted
/// histogram with events that have constant weights equal to 1.
/// The variance \f$ z_{i}^{2} \f$ of the difference between the weight wi
/// and the estimated expectation value of the weight is approximately equal to:
///\f[
///  z_{i}^{2} = Var(w_{i}-W\hat{p}_{i}) = N\hat{p}_{i}(1-N\hat{p}_{i})\left(\frac{Ws_{i}^{2}}{\sqrt{(Ns_{i}^{2}-w_{i}W)^{2}+4W^{2}s_{i}^{2}n_{i}}}\right)^{2}+\frac{s_{i}^{2}}{4}\left(1+\frac{Ns_{i}^{2}-w_{i}W}{\sqrt{(Ns_{i}^{2}-w_{i}W)^{2}+4W^{2}s_{i}^{2}n_{i}}}\right)^{2}
///\f]
/// The  residuals
///\f[
///  r_{i} = \frac{w_{i}-W\hat{p}_{i}}{z_{i}}
///\f]
/// have approximately a normal distribution with mean equal to 0 and standard
/// deviation  equal to 1.
///
/// #### Two weighted histograms comparison:
///
/// Let us denote the common  weight of events of the ith bin in the first
/// histogram as w1i and as w2i in the second one. The total weight of events
/// in the first histogram is equal to
///\f[
///  W_{1} = \sum_{i=1}^{r} w_{1i}
///\f]
/// and
///\f[
///  W_{2} = \sum_{i=1}^{r} w_{2i}
///\f]
/// in the second histogram. Let us formulate the hypothesis of identity of
/// weighted histograms so that there exist r constants p1,...,pr, such that
///\f[
///  \sum_{i=1}^{r} p_{i} = 1
///\f]
/// and also expectation value of weight w1i equal to W1pi and expectation value
/// of weight w2i equal to W2pi. Weights in both the histograms are random
/// variables with distributions which can be approximated by a normal
/// probability distribution \f$ N(W_{1}p_{i},\sigma_{1i}^{2}) \f$ for the first histogram
/// and by a distribution \f$ N(W_{2}p_{i},\sigma_{2i}^{2}) \f$ for the second.
/// Here \f$ \sigma_{1i}^{2} \f$ and \f$ \sigma_{2i}^{2} \f$ are the variances
/// of w1i and w2i with estimators \f$ s_{1i}^{2} \f$ and \f$ s_{2i}^{2} \f$ respectively.
/// If the hypothesis of identity is valid, then the maximum likelihood and
/// Least Square Method estimator of pi,i=1,...,r, is
///\f[
///  \hat{p}_{i} = \frac{w_{1i}W_{1}/s_{1i}^{2}+w_{2i}W_{2} /s_{2i}^{2}}{W_{1}^{2}/s_{1i}^{2}+W_{2}^{2}/s_{2i}^{2}}
///\f]
/// We may then use the test statistic
///\f[
/// X^{2} = \sum_{i=1}^{r} \frac{(w_{1i}-W_{1}\hat{p}_{i})^{2}}{s_{1i}^{2}} + \sum_{i=1}^{r} \frac{(w_{2i}-W_{2}\hat{p}_{i})^{2}}{s_{2i}^{2}} = \sum_{i=1}^{r} \frac{(W_{1}w_{2i}-W_{2}w_{1i})^{2}}{W_{1}^{2}s_{2i}^{2}+W_{2}^{2}s_{1i}^{2}}
///\f]
/// and it has approximately a \f$ \chi^{2}_{(r-1)} \f$ distribution [2].
/// The normalized or studentised residuals [6]
///\f[
///  r_{i} = \frac{w_{1i}-W_{1}\hat{p}_{i}}{s_{1i}\sqrt{1 - \frac{1}{(1+W_{2}^{2}s_{1i}^{2}/W_{1}^{2}s_{2i}^{2})}}}
///\f]
/// have approximately a normal distribution with mean equal to 0 and standard
/// deviation 1. A recommended minimal expected frequency is equal to 10 for
/// the proposed test.
///
/// #### Numerical examples:
///
/// The method described herein is now illustrated with an example.
/// We take a distribution
///\f[
/// \phi(x) = \frac{2}{(x-10)^{2}+1} + \frac{1}{(x-14)^{2}+1}       (1)
///\f]
/// defined on the interval [4,16]. Events distributed according to the formula
/// (1) are simulated to create the unweighted histogram. Uniformly distributed
/// events are simulated for the weighted histogram with weights calculated by
/// formula (1). Each histogram has the same number of bins: 20. Fig.1 shows
/// the result of comparison of the unweighted histogram with 200 events
/// (minimal expected frequency equal to one) and the weighted histogram with
/// 500 events (minimal expected frequency equal to 25)
/// Begin_Macro
/// ../../../tutorials/math/chi2test.C
/// End_Macro
/// Fig 1. An example of comparison of the unweighted histogram with 200 events
/// and the weighted histogram with 500 events:
///   1. unweighted histogram;
///   2. weighted histogram;
///   3. normalized residuals plot;
///   4. normal Q-Q plot of residuals.
///
/// The value of the test statistic \f$ \chi^{2} \f$ is equal to
/// 21.09 with p-value equal to 0.33, therefore the hypothesis of identity of
/// the two histograms can be accepted for 0.05 significant level. The behavior
/// of the normalized residuals plot (see Fig. 1c) and the normal Q-Q plot
/// (see Fig. 1d) of residuals are regular and we cannot identify the outliers
/// or bins with a big influence on \f$ \chi^{2} \f$.
///
/// The second example presents the same two histograms but 17 events was added
/// to content of bin number 15 in unweighted histogram. Fig.2 shows the result
/// of comparison of the unweighted histogram with 217 events (minimal expected
/// frequency equal to one) and the weighted histogram with 500 events (minimal
/// expected frequency equal to 25)
/// Begin_Macro
/// ../../../tutorials/math/chi2test.C(17)
/// End_Macro
/// Fig 2. An example of comparison of the unweighted histogram with 217 events
/// and the weighted histogram with 500 events:
///   1. unweighted histogram;
///   2. weighted histogram;
///   3. normalized residuals plot;
///   4. normal Q-Q plot of residuals.
///
/// The value of the test statistic \f$ \chi^{2} \f$ is equal to
/// 32.33 with p-value equal to 0.029, therefore the hypothesis of identity of
/// the two histograms is rejected for 0.05 significant level. The behavior of
/// the normalized residuals plot (see Fig. 2c) and the normal Q-Q plot (see
/// Fig. 2d) of residuals are not regular and we can identify the outlier or
/// bin with a big influence on \f$ \chi^{2} \f$.
///
/// #### References:
///
///  - [1] Pearson, K., 1904. On the Theory of Contingency and Its Relation to
///    Association and Normal Correlation. Drapers' Co. Memoirs, Biometric
///    Series No. 1, London.
///  - [2] Gagunashvili, N., 2006. \f$ \sigma^{2} \f$ test for comparison
///    of weighted and unweighted histograms. Statistical Problems in Particle
///    Physics, Astrophysics and Cosmology, Proceedings of PHYSTAT05,
///    Oxford, UK, 12-15 September 2005, Imperial College Press, London, 43-44.
///    Gagunashvili,N., Comparison of weighted and unweighted histograms,
///    arXiv:physics/0605123, 2006.
///  - [3] Cramer, H., 1946. Mathematical methods of statistics.
///    Princeton University Press, Princeton.
///  - [4] Haberman, S.J., 1973. The analysis of residuals in cross-classified tables.
///    Biometrics 29, 205-220.
///  - [5] Lewontin, R.C. and Felsenstein, J., 1965. The robustness of homogeneity
///    test in 2xN tables. Biometrics 21, 19-33.
///  - [6] Seber, G.A.F., Lee, A.J., 2003, Linear Regression Analysis.
///    John Wiley & Sons Inc., New York.
 
Double_t TH1::Chi2Test(const TH1* h2, Option_t *option, Double_t *res) const
{
   Double_t chi2 = 0;
   Int_t ndf = 0, igood = 0;
 
   TString opt = option;
   opt.ToUpper();
 
   Double_t prob = Chi2TestX(h2,chi2,ndf,igood,option,res);
 
   if(opt.Contains("P")) {
      printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
   }
   if(opt.Contains("CHI2/NDF")) {
      if (ndf == 0) return 0;
      return chi2/ndf;
   }
   if(opt.Contains("CHI2")) {
      return chi2;
   }
 
   return prob;
}
 
////////////////////////////////////////////////////////////////////////////////
/// The computation routine of the Chisquare test. For the method description,
/// see Chi2Test() function.
///
/// \return p-value
/// \param[in] h2 the second histogram
/// \param[in] option
///  - "UU" = experiment experiment comparison (unweighted-unweighted)
///  - "UW" = experiment MC comparison (unweighted-weighted). Note that the first
///        histogram should be unweighted
///  - "WW" = MC MC comparison (weighted-weighted)
///  - "NORM" = if one or both histograms is scaled
///  - "OF" = overflows included
///  - "UF" = underflows included
///      by default underflows and overflows are not included
/// \param[out] igood test output
///    - igood=0 - no problems
///    - For unweighted unweighted  comparison
///      - igood=1'There is a bin in the 1st histogram with less than 1 event'
///      - igood=2'There is a bin in the 2nd histogram with less than 1 event'
///      - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
///    - For  unweighted weighted  comparison
///      - igood=1'There is a bin in the 1st histogram with less then 1 event'
///      - igood=2'There is a bin in the 2nd histogram with less then 10 effective number of events'
///      - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
///    - For  weighted weighted  comparison
///      - igood=1'There is a bin in the 1st  histogram with less then 10 effective
///        number of events'
///      - igood=2'There is a bin in the 2nd  histogram with less then 10 effective
///        number of events'
///      - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
/// \param[out] chi2 chisquare of the test
/// \param[out] ndf number of degrees of freedom (important, when both histograms have the same empty bins)
/// \param[out] res normalized residuals for further analysis
 
Double_t TH1::Chi2TestX(const TH1* h2,  Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option,  Double_t *res) const
{
 
   Int_t i_start, i_end;
   Int_t j_start, j_end;
   Int_t k_start, k_end;
 
   Double_t sum1 = 0.0, sumw1 = 0.0;
   Double_t sum2 = 0.0, sumw2 = 0.0;
 
   chi2 = 0.0;
   ndf = 0;
 
   TString opt = option;
   opt.ToUpper();
 
   if (fBuffer) const_cast<TH1*>(this)->BufferEmpty();
 
   const TAxis *xaxis1 = GetXaxis();
   const TAxis *xaxis2 = h2->GetXaxis();
   const TAxis *yaxis1 = GetYaxis();
   const TAxis *yaxis2 = h2->GetYaxis();
   const TAxis *zaxis1 = GetZaxis();
   const TAxis *zaxis2 = h2->GetZaxis();
 
   Int_t nbinx1 = xaxis1->GetNbins();
   Int_t nbinx2 = xaxis2->GetNbins();
   Int_t nbiny1 = yaxis1->GetNbins();
   Int_t nbiny2 = yaxis2->GetNbins();
   Int_t nbinz1 = zaxis1->GetNbins();
   Int_t nbinz2 = zaxis2->GetNbins();
 
   //check dimensions
   if (this->GetDimension() != h2->GetDimension() ){
      Error("Chi2TestX","Histograms have different dimensions.");
      return 0.0;
   }
 
   //check number of channels
   if (nbinx1 != nbinx2) {
      Error("Chi2TestX","different number of x channels");
   }
   if (nbiny1 != nbiny2) {
      Error("Chi2TestX","different number of y channels");
   }
   if (nbinz1 != nbinz2) {
      Error("Chi2TestX","different number of z channels");
   }
 
   //check for ranges
   i_start = j_start = k_start = 1;
   i_end = nbinx1;
   j_end = nbiny1;
   k_end = nbinz1;
 
   if (xaxis1->TestBit(TAxis::kAxisRange)) {
      i_start = xaxis1->GetFirst();
      i_end   = xaxis1->GetLast();
   }
   if (yaxis1->TestBit(TAxis::kAxisRange)) {
      j_start = yaxis1->GetFirst();
      j_end   = yaxis1->GetLast();
   }
   if (zaxis1->TestBit(TAxis::kAxisRange)) {
      k_start = zaxis1->GetFirst();
      k_end   = zaxis1->GetLast();
   }
 
 
   if (opt.Contains("OF")) {
      if (GetDimension() == 3) k_end = ++nbinz1;
      if (GetDimension() >= 2) j_end = ++nbiny1;
      if (GetDimension() >= 1) i_end = ++nbinx1;
   }
 
   if (opt.Contains("UF")) {
      if (GetDimension() == 3) k_start = 0;
      if (GetDimension() >= 2) j_start = 0;
      if (GetDimension() >= 1) i_start = 0;
   }
 
   ndf = (i_end - i_start + 1) * (j_end - j_start + 1) * (k_end - k_start + 1) - 1;
 
   Bool_t comparisonUU = opt.Contains("UU");
   Bool_t comparisonUW = opt.Contains("UW");
   Bool_t comparisonWW = opt.Contains("WW");
   Bool_t scaledHistogram  = opt.Contains("NORM");
 
   if (scaledHistogram && !comparisonUU) {
      Info("Chi2TestX", "NORM option should be used together with UU option. It is ignored");
   }
 
   // look at histo global bin content and effective entries
   Stat_t s[kNstat];
   GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
   Double_t sumBinContent1 = s[0];
   Double_t effEntries1 = (s[1] ? s[0] * s[0] / s[1] : 0.0);
 
   h2->GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
   Double_t sumBinContent2 = s[0];
   Double_t effEntries2 = (s[1] ? s[0] * s[0] / s[1] : 0.0);
 
   if (!comparisonUU && !comparisonUW && !comparisonWW ) {
      // deduce automatically from type of histogram
      if (TMath::Abs(sumBinContent1 - effEntries1) < 1) {
         if ( TMath::Abs(sumBinContent2 - effEntries2) < 1) comparisonUU = true;
         else comparisonUW = true;
      }
      else comparisonWW = true;
   }
   // check unweighted histogram
   if (comparisonUW) {
      if (TMath::Abs(sumBinContent1 - effEntries1) >= 1) {
         Warning("Chi2TestX","First histogram is not unweighted and option UW has been requested");
      }
   }
   if ( (!scaledHistogram && comparisonUU)   ) {
      if ( ( TMath::Abs(sumBinContent1 - effEntries1) >= 1) || (TMath::Abs(sumBinContent2 - effEntries2) >= 1) ) {
         Warning("Chi2TestX","Both histograms are not unweighted and option UU has been requested");
      }
   }
 
 
   //get number of events in histogram
   if (comparisonUU && scaledHistogram) {
      for (Int_t i = i_start; i <= i_end; ++i) {
         for (Int_t j = j_start; j <= j_end; ++j) {
            for (Int_t k = k_start; k <= k_end; ++k) {
 
               Int_t bin = GetBin(i, j, k);
 
               Double_t cnt1 = RetrieveBinContent(bin);
               Double_t cnt2 = h2->RetrieveBinContent(bin);
               Double_t e1sq = GetBinErrorSqUnchecked(bin);
               Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
 
               if (e1sq > 0.0) cnt1 = TMath::Floor(cnt1 * cnt1 / e1sq + 0.5); // avoid rounding errors
               else cnt1 = 0.0;
 
               if (e2sq > 0.0) cnt2 = TMath::Floor(cnt2 * cnt2 / e2sq + 0.5); // avoid rounding errors
               else cnt2 = 0.0;
 
               // sum contents
               sum1 += cnt1;
               sum2 += cnt2;
               sumw1 += e1sq;
               sumw2 += e2sq;
            }
         }
      }
      if (sumw1 <= 0.0 || sumw2 <= 0.0) {
         Error("Chi2TestX", "Cannot use option NORM when one histogram has all zero errors");
         return 0.0;
      }
 
   } else {
      for (Int_t i = i_start; i <= i_end; ++i) {
         for (Int_t j = j_start; j <= j_end; ++j) {
            for (Int_t k = k_start; k <= k_end; ++k) {
 
               Int_t bin = GetBin(i, j, k);
 
               sum1 += RetrieveBinContent(bin);
               sum2 += h2->RetrieveBinContent(bin);
 
               if ( comparisonWW ) sumw1 += GetBinErrorSqUnchecked(bin);
               if ( comparisonUW || comparisonWW ) sumw2 += h2->GetBinErrorSqUnchecked(bin);
            }
         }
      }
   }
   //checks that the histograms are not empty
   if (sum1 == 0.0 || sum2 == 0.0) {
      Error("Chi2TestX","one histogram is empty");
      return 0.0;
   }
 
   if ( comparisonWW  && ( sumw1 <= 0.0 && sumw2 <= 0.0 ) ){
      Error("Chi2TestX","Hist1 and Hist2 have both all zero errors\n");
      return 0.0;
   }
 
   //THE TEST
   Int_t m = 0, n = 0;
 
   //Experiment - experiment comparison
   if (comparisonUU) {
      Double_t sum = sum1 + sum2;
      for (Int_t i = i_start; i <= i_end; ++i) {
         for (Int_t j = j_start; j <= j_end; ++j) {
            for (Int_t k = k_start; k <= k_end; ++k) {
 
               Int_t bin = GetBin(i, j, k);
 
               Double_t cnt1 = RetrieveBinContent(bin);
               Double_t cnt2 = h2->RetrieveBinContent(bin);
 
               if (scaledHistogram) {
                  // scale bin value to effective bin entries
                  Double_t e1sq = GetBinErrorSqUnchecked(bin);
                  Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
 
                  if (e1sq > 0) cnt1 = TMath::Floor(cnt1 * cnt1 / e1sq + 0.5); // avoid rounding errors
                  else cnt1 = 0;
 
                  if (e2sq > 0) cnt2 = TMath::Floor(cnt2 * cnt2 / e2sq + 0.5); // avoid rounding errors
                  else cnt2 = 0;
               }
 
               if (Int_t(cnt1) == 0 && Int_t(cnt2) == 0) --ndf;  // no data means one degree of freedom less
               else {
 
                  Double_t cntsum = cnt1 + cnt2;
                  Double_t nexp1 = cntsum * sum1 / sum;
                  //Double_t nexp2 = binsum*sum2/sum;
 
                  if (res) res[i - i_start] = (cnt1 - nexp1) / TMath::Sqrt(nexp1);
 
                  if (cnt1 < 1) ++m;
                  if (cnt2 < 1) ++n;
 
                  //Habermann correction for residuals
                  Double_t correc = (1. - sum1 / sum) * (1. - cntsum / sum);
                  if (res) res[i - i_start] /= TMath::Sqrt(correc);
 
                  Double_t delta = sum2 * cnt1 - sum1 * cnt2;
                  chi2 += delta * delta / cntsum;
               }
            }
         }
      }
      chi2 /= sum1 * sum2;
 
      // flag error only when of the two histogram is zero
      if (m) {
         igood += 1;
         Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
      }
      if (n) {
         igood += 2;
         Info("Chi2TestX","There is a bin in h2 with less than 1 event.\n");
      }
 
      Double_t prob = TMath::Prob(chi2,ndf);
      return prob;
 
   }
 
   // unweighted - weighted  comparison
   // case of error = 0 and content not zero is treated without problems by excluding second chi2 sum
   // and can be considered as a data-theory comparison
   if ( comparisonUW ) {
      for (Int_t i = i_start; i <= i_end; ++i) {
         for (Int_t j = j_start; j <= j_end; ++j) {
            for (Int_t k = k_start; k <= k_end; ++k) {
 
               Int_t bin = GetBin(i, j, k);
 
               Double_t cnt1 = RetrieveBinContent(bin);
               Double_t cnt2 = h2->RetrieveBinContent(bin);
               Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
 
               // case both histogram have zero bin contents
               if (cnt1 * cnt1 == 0 && cnt2 * cnt2 == 0) {
                  --ndf;  //no data means one degree of freedom less
                  continue;
               }
 
               // case weighted histogram has zero bin content and error
               if (cnt2 * cnt2 == 0 && e2sq == 0) {
                  if (sumw2 > 0) {
                     // use as approximated  error as 1 scaled by a scaling ratio
                     // estimated from the total sum weight and sum weight squared
                     e2sq = sumw2 / sum2;
                  }
                  else {
                     // return error because infinite discrepancy here:
                     // bin1 != 0 and bin2 =0 in a histogram with all errors zero
                     Error("Chi2TestX","Hist2 has in bin (%d,%d,%d) zero content and zero errors\n", i, j, k);
                     chi2 = 0; return 0;
                  }
               }
 
               if (cnt1 < 1) m++;
               if (e2sq > 0 && cnt2 * cnt2 / e2sq < 10) n++;
 
               Double_t var1 = sum2 * cnt2 - sum1 * e2sq;
               Double_t var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
 
               // if cnt1 is zero and cnt2 = 1 and sum1 = sum2 var1 = 0 && var2 == 0
               // approximate by incrementing cnt1
               // LM (this need to be fixed for numerical errors)
               while (var1 * var1 + cnt1 == 0 || var1 + var2 == 0) {
                  sum1++;
                  cnt1++;
                  var1 = sum2 * cnt2 - sum1 * e2sq;
                  var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
               }
               var2 = TMath::Sqrt(var2);
 
               while (var1 + var2 == 0) {
                  sum1++;
                  cnt1++;
                  var1 = sum2 * cnt2 - sum1 * e2sq;
                  var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
                  while (var1 * var1 + cnt1 == 0 || var1 + var2 == 0) {
                     sum1++;
                     cnt1++;
                     var1 = sum2 * cnt2 - sum1 * e2sq;
                     var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
                  }
                  var2 = TMath::Sqrt(var2);
               }
 
               Double_t probb = (var1 + var2) / (2. * sum2 * sum2);
 
               Double_t nexp1 = probb * sum1;
               Double_t nexp2 = probb * sum2;
 
               Double_t delta1 = cnt1 - nexp1;
               Double_t delta2 = cnt2 - nexp2;
 
               chi2 += delta1 * delta1 / nexp1;
 
               if (e2sq > 0) {
                  chi2 += delta2 * delta2 / e2sq;
               }
 
               if (res) {
                  if (e2sq > 0) {
                     Double_t temp1 = sum2 * e2sq / var2;
                     Double_t temp2 = 1.0 + (sum1 * e2sq - sum2 * cnt2) / var2;
                     temp2 = temp1 * temp1 * sum1 * probb * (1.0 - probb) + temp2 * temp2 * e2sq / 4.0;
                     // invert sign here
                     res[i - i_start] = - delta2 / TMath::Sqrt(temp2);
                  }
                  else
                     res[i - i_start] = delta1 / TMath::Sqrt(nexp1);
               }
            }
         }
      }
 
      if (m) {
         igood += 1;
         Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
      }
      if (n) {
         igood += 2;
         Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
      }
 
      Double_t prob = TMath::Prob(chi2, ndf);
 
      return prob;
   }
 
   // weighted - weighted  comparison
   if (comparisonWW) {
      for (Int_t i = i_start; i <= i_end; ++i) {
         for (Int_t j = j_start; j <= j_end; ++j) {
            for (Int_t k = k_start; k <= k_end; ++k) {
 
               Int_t bin = GetBin(i, j, k);
               Double_t cnt1 = RetrieveBinContent(bin);
               Double_t cnt2 = h2->RetrieveBinContent(bin);
               Double_t e1sq = GetBinErrorSqUnchecked(bin);
               Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
 
               // case both histogram have zero bin contents
               // (use square of content to avoid numerical errors)
                if (cnt1 * cnt1 == 0 && cnt2 * cnt2 == 0) {
                   --ndf;  //no data means one degree of freedom less
                   continue;
                }
 
                if (e1sq == 0 && e2sq == 0) {
                   // cannot treat case of booth histogram have zero zero errors
                  Error("Chi2TestX","h1 and h2 both have bin %d,%d,%d with all zero errors\n", i,j,k);
                  chi2 = 0; return 0;
               }
 
               Double_t sigma = sum1 * sum1 * e2sq + sum2 * sum2 * e1sq;
               Double_t delta = sum2 * cnt1 - sum1 * cnt2;
               chi2 += delta * delta / sigma;
 
               if (res) {
                  Double_t temp = cnt1 * sum1 * e2sq + cnt2 * sum2 * e1sq;
                  Double_t probb = temp / sigma;
                  Double_t z = 0;
                  if (e1sq > e2sq) {
                     Double_t d1 = cnt1 - sum1 * probb;
                     Double_t s1 = e1sq * ( 1. - e2sq * sum1 * sum1 / sigma );
                     z = d1 / TMath::Sqrt(s1);
                  }
                  else {
                     Double_t d2 = cnt2 - sum2 * probb;
                     Double_t s2 = e2sq * ( 1. - e1sq * sum2 * sum2 / sigma );
                     z = -d2 / TMath::Sqrt(s2);
                  }
                  res[i - i_start] = z;
               }
 
               if (e1sq > 0 && cnt1 * cnt1 / e1sq < 10) m++;
               if (e2sq > 0 && cnt2 * cnt2 / e2sq < 10) n++;
            }
         }
      }
      if (m) {
         igood += 1;
         Info("Chi2TestX","There is a bin in h1 with less than 10 effective events.\n");
      }
      if (n) {
         igood += 2;
         Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
      }
      Double_t prob = TMath::Prob(chi2, ndf);
      return prob;
   }
   return 0;
}
////////////////////////////////////////////////////////////////////////////////
/// Compute and return the chisquare of this histogram with respect to a function
/// The chisquare is computed by weighting each histogram point by the bin error
/// By default the full range of the histogram is used.
/// Use option "R" for restricting the chisquare calculation to the given range of the function
/// Use option "L" for using the chisquare based on the poisson likelihood (Baker-Cousins Chisquare)
/// Use option "P" for using the Pearson chisquare based on the expected bin errors
 
Double_t TH1::Chisquare(TF1 * func, Option_t *option) const
{
   if (!func) {
      Error("Chisquare","Function pointer is Null - return -1");
      return -1;
   }
 
   TString opt(option); opt.ToUpper();
   bool useRange = opt.Contains("R");
   ROOT::Fit::EChisquareType type = ROOT::Fit::EChisquareType::kNeyman;  // default chi2 with observed error
   if (opt.Contains("L")) type = ROOT::Fit::EChisquareType::kPLikeRatio;
   else if (opt.Contains("P")) type = ROOT::Fit::EChisquareType::kPearson;
 
   return ROOT::Fit::Chisquare(*this, *func, useRange, type);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Remove all the content from the underflow and overflow bins, without changing the number of entries
/// After calling this method, every undeflow and overflow bins will have content 0.0
/// The Sumw2 is also cleared, since there is no more content in the bins
 
void TH1::ClearUnderflowAndOverflow()
{
   for (Int_t bin = 0; bin < fNcells; ++bin)
      if (IsBinUnderflow(bin) || IsBinOverflow(bin)) {
         UpdateBinContent(bin, 0.0);
         if (fSumw2.fN) fSumw2.fArray[bin] = 0.0;
      }
}
 
////////////////////////////////////////////////////////////////////////////////
///  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
///  If the routine is called with the onlyPositive flag set an error will
///  be produced in case of negative bin content and a NaN value returned
 
Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
{
   if (fBuffer) BufferEmpty();
 
   // 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 nbins  = nbinsx * nbinsy * nbinsz;
 
   fIntegral = new Double_t[nbins + 2];
   Int_t ibin = 0; fIntegral[ibin] = 0;
 
   for (Int_t binz=1; binz <= nbinsz; ++binz) {
      for (Int_t biny=1; biny <= nbinsy; ++biny) {
         for (Int_t binx=1; binx <= nbinsx; ++binx) {
            ++ibin;
            Double_t y = RetrieveBinContent(GetBin(binx, biny, binz));
            if (onlyPositive && y < 0) {
                 Error("ComputeIntegral","Bin content is negative - return a NaN value");
                 fIntegral[nbins] = TMath::QuietNaN();
                 break;
             }
            fIntegral[ibin] = fIntegral[ibin - 1] + y;
         }
      }
   }
 
   //   - Normalize integral to 1
   if (fIntegral[nbins] == 0 ) {
      Error("ComputeIntegral", "Integral = zero"); return 0;
   }
   for (Int_t bin=1; bin <= nbins; ++bin)  fIntegral[bin] /= fIntegral[nbins];
   fIntegral[nbins+1] = fEntries;
   return fIntegral[nbins];
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return a pointer to the array of bins integral.
/// if the pointer fIntegral is null, TH1::ComputeIntegral is called
/// The array dimension is the number of bins in the histograms
/// including underflow and overflow (fNCells)
/// the last value integral[fNCells] is set to the number of entries of
/// the histogram
 
Double_t *TH1::GetIntegral()
{
   if (!fIntegral) ComputeIntegral();
   return fIntegral;
}
 
////////////////////////////////////////////////////////////////////////////////
///  Return a pointer to a histogram containing the cumulative content.
///  The cumulative can be computed both in the forward (default) or backward
///  direction; the name of the new histogram is constructed from
///  the name of this histogram with the suffix "suffix" appended provided
///  by the user. If not provided a default suffix="_cumulative" is used.
///
/// The cumulative distribution is formed by filling each bin of the
/// resulting histogram with the sum of that bin and all previous
/// (forward == kTRUE) or following (forward = kFALSE) bins.
///
/// Note: while cumulative distributions make sense in one dimension, you
/// may not be getting what you expect in more than 1D because the concept
/// of a cumulative distribution is much trickier to define; make sure you
/// understand the order of summation before you use this method with
/// histograms of dimension >= 2.
///
/// Note 2: By default the cumulative is computed from bin 1 to Nbins
/// If an axis range is set, values between the minimum and maximum of the range
/// are set.
/// Setting an axis range can also be used for including underflow and overflow in
/// the cumulative (e.g. by setting h->GetXaxis()->SetRange(0, h->GetNbinsX()+1); )
///
 
TH1 *TH1::GetCumulative(Bool_t forward, const char* suffix) const
{
   const Int_t firstX = fXaxis.GetFirst();
   const Int_t lastX  = fXaxis.GetLast();
   const Int_t firstY = (fDimension > 1) ? fYaxis.GetFirst() : 1;
   const Int_t lastY = (fDimension > 1) ? fYaxis.GetLast() : 1;
   const Int_t firstZ = (fDimension > 1) ? fZaxis.GetFirst() : 1;
   const Int_t lastZ = (fDimension > 1) ? fZaxis.GetLast() : 1;
 
   TH1* hintegrated = (TH1*) Clone(fName + suffix);
   hintegrated->Reset();
   Double_t sum = 0.;
   Double_t esum = 0;
   if (forward) { // Forward computation
      for (Int_t binz = firstZ; binz <= lastZ; ++binz) {
         for (Int_t biny = firstY; biny <= lastY; ++biny) {
            for (Int_t binx = firstX; binx <= lastX; ++binx) {
               const Int_t bin = hintegrated->GetBin(binx, biny, binz);
               sum += RetrieveBinContent(bin);
               hintegrated->AddBinContent(bin, sum);
               if (fSumw2.fN) {
                  esum += GetBinErrorSqUnchecked(bin);
                  hintegrated->fSumw2.fArray[bin] = esum;
               }
            }
         }
      }
   } else { // Backward computation
      for (Int_t binz = lastZ; binz >= firstZ; --binz) {
         for (Int_t biny = lastY; biny >= firstY; --biny) {
            for (Int_t binx = lastX; binx >= firstX; --binx) {
               const Int_t bin = hintegrated->GetBin(binx, biny, binz);
               sum += RetrieveBinContent(bin);
               hintegrated->AddBinContent(bin, sum);
               if (fSumw2.fN) {
                  esum += GetBinErrorSqUnchecked(bin);
                  hintegrated->fSumw2.fArray[bin] = esum;
               }
            }
         }
      }
   }
   return hintegrated;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Copy this histogram structure to newth1.
///
/// Note that this function does not copy the list of associated functions.
/// Use TObject::Clone to make a full copy of a histogram.
///
/// Note also that the histogram it will be created in gDirectory (if AddDirectoryStatus()=true)
/// or will not be added to any directory if  AddDirectoryStatus()=false
/// independently of the current directory stored in the original histogram
 
void TH1::Copy(TObject &obj) const
{
   if (((TH1&)obj).fDirectory) {
      // We are likely to change the hash value of this object
      // with TNamed::Copy, to keep things correct, we need to
      // clean up its existing entries.
      ((TH1&)obj).fDirectory->Remove(&obj);
      ((TH1&)obj).fDirectory = nullptr;
   }
   TNamed::Copy(obj);
   ((TH1&)obj).fDimension = fDimension;
   ((TH1&)obj).fNormFactor= fNormFactor;
   ((TH1&)obj).fNcells    = fNcells;
   ((TH1&)obj).fBarOffset = fBarOffset;
   ((TH1&)obj).fBarWidth  = fBarWidth;
   ((TH1&)obj).fOption    = fOption;
   ((TH1&)obj).fBinStatErrOpt = fBinStatErrOpt;
   ((TH1&)obj).fBufferSize= fBufferSize;
   // copy the Buffer
   // delete first a previously existing buffer
   if (((TH1&)obj).fBuffer != nullptr)  {
      delete [] ((TH1&)obj).fBuffer;
      ((TH1&)obj).fBuffer = nullptr;
   }
   if (fBuffer) {
      Double_t *buf = new Double_t[fBufferSize];
      for (Int_t i=0;i<fBufferSize;i++) buf[i] = fBuffer[i];
      // obj.fBuffer has been deleted before
      ((TH1&)obj).fBuffer    = buf;
   }
 
 
   TArray* a = dynamic_cast<TArray*>(&obj);
   if (a) a->Set(fNcells);
   for (Int_t i = 0; i < fNcells; i++) ((TH1&)obj).UpdateBinContent(i, RetrieveBinContent(i));
 
   ((TH1&)obj).fEntries   = fEntries;
 
   // which will call BufferEmpty(0) and set fBuffer[0] to a  Maybe one should call
   // assignment operator on the TArrayD
 
   ((TH1&)obj).fTsumw     = fTsumw;
   ((TH1&)obj).fTsumw2    = fTsumw2;
   ((TH1&)obj).fTsumwx    = fTsumwx;
   ((TH1&)obj).fTsumwx2   = fTsumwx2;
   ((TH1&)obj).fMaximum   = fMaximum;
   ((TH1&)obj).fMinimum   = fMinimum;
 
   TAttLine::Copy(((TH1&)obj));
   TAttFill::Copy(((TH1&)obj));
   TAttMarker::Copy(((TH1&)obj));
   fXaxis.Copy(((TH1&)obj).fXaxis);
   fYaxis.Copy(((TH1&)obj).fYaxis);
   fZaxis.Copy(((TH1&)obj).fZaxis);
   ((TH1&)obj).fXaxis.SetParent(&obj);
   ((TH1&)obj).fYaxis.SetParent(&obj);
   ((TH1&)obj).fZaxis.SetParent(&obj);
   fContour.Copy(((TH1&)obj).fContour);
   fSumw2.Copy(((TH1&)obj).fSumw2);
   //   fFunctions->Copy(((TH1&)obj).fFunctions);
   // when copying an histogram if the AddDirectoryStatus() is true it
   // will be added to gDirectory independently of the fDirectory stored.
   // and if the AddDirectoryStatus() is false it will not be added to
   // any directory (fDirectory = nullptr)
   if (fgAddDirectory && gDirectory) {
      gDirectory->Append(&obj);
      ((TH1&)obj).fFunctions->UseRWLock();
      ((TH1&)obj).fDirectory = gDirectory;
   } else
      ((TH1&)obj).fDirectory = nullptr;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Make a complete copy of the underlying object.  If 'newname' is set,
/// the copy's name will be set to that name.
 
TObject* TH1::Clone(const char* newname) const
{
   TH1* obj = (TH1*)IsA()->GetNew()(nullptr);
   Copy(*obj);
 
   // Now handle the parts that Copy doesn't do
   if(fFunctions) {
      // The Copy above might have published 'obj' to the ListOfCleanups.
      // Clone can call RecursiveRemove, for example via TCheckHashRecursiveRemoveConsistency
      // when dictionary information is initialized, so we need to
      // keep obj->fFunction valid during its execution and
      // protect the update with the write lock.
 
      // Reset stats parent - else cloning the stats will clone this histogram, too.
      auto oldstats = dynamic_cast<TVirtualPaveStats*>(fFunctions->FindObject("stats"));
      TObject *oldparent = nullptr;
      if (oldstats) {
         oldparent = oldstats->GetParent();
         oldstats->SetParent(nullptr);
      }
 
      auto newlist = (TList*)fFunctions->Clone();
 
      if (oldstats)
         oldstats->SetParent(oldparent);
      auto newstats = dynamic_cast<TVirtualPaveStats*>(obj->fFunctions->FindObject("stats"));
      if (newstats)
         newstats->SetParent(obj);
 
      auto oldlist = obj->fFunctions;
      {
         R__WRITE_LOCKGUARD(ROOT::gCoreMutex);
         obj->fFunctions = newlist;
      }
      delete oldlist;
   }
   if(newname && strlen(newname) ) {
      obj->SetName(newname);
   }
   return obj;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Perform the automatic addition of the histogram to the given directory
///
/// Note this function is called in place when the semantic requires
/// this object to be added to a directory (I.e. when being read from
/// a TKey or being Cloned)
 
void TH1::DirectoryAutoAdd(TDirectory *dir)
{
   Bool_t addStatus = TH1::AddDirectoryStatus();
   if (addStatus) {
      SetDirectory(dir);
      if (dir) {
         ResetBit(kCanDelete);
      }
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute distance from point px,py to a line.
///
///  Compute the closest distance of approach from point px,py to elements
///  of a histogram.
///  The distance is computed in pixels units.
///
///  #### Algorithm:
///  Currently, this simple model computes the distance from the mouse
///  to the histogram contour only.
 
Int_t TH1::DistancetoPrimitive(Int_t px, Int_t py)
{
   if (!fPainter) return 9999;
   return fPainter->DistancetoPrimitive(px,py);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Performs the operation: `this = this/(c1*f1)`
/// if errors are defined (see TH1::Sumw2), errors are also recalculated.
///
/// Only bins inside the function range are recomputed.
/// IMPORTANT NOTE: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Divide
///
/// The function return kFALSE if the divide operation failed
 
Bool_t TH1::Divide(TF1 *f1, Double_t c1)
{
   if (!f1) {
      Error("Divide","Attempt to divide by a non-existing function");
      return kFALSE;
   }
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   Int_t nx = GetNbinsX() + 2; // normal bins + uf / of
   Int_t ny = GetNbinsY() + 2;
   Int_t nz = GetNbinsZ() + 2;
   if (fDimension < 2) ny = 1;
   if (fDimension < 3) nz = 1;
 
 
   SetMinimum();
   SetMaximum();
 
   //   - Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu, w;
   Double_t xx[3];
   Double_t *params = nullptr;
   f1->InitArgs(xx,params);
   for (binz = 0; binz < nz; ++binz) {
      xx[2] = fZaxis.GetBinCenter(binz);
      for (biny = 0; biny < ny; ++biny) {
         xx[1] = fYaxis.GetBinCenter(biny);
         for (binx = 0; binx < nx; ++binx) {
            xx[0] = fXaxis.GetBinCenter(binx);
            if (!f1->IsInside(xx)) continue;
            TF1::RejectPoint(kFALSE);
            bin = binx + nx * (biny + ny * binz);
            cu  = c1 * f1->EvalPar(xx);
            if (TF1::RejectedPoint()) continue;
            if (cu) w = RetrieveBinContent(bin) / cu;
            else    w = 0;
            UpdateBinContent(bin, w);
            if (fSumw2.fN) {
               if (cu != 0) fSumw2.fArray[bin] = GetBinErrorSqUnchecked(bin) / (cu * cu);
               else         fSumw2.fArray[bin] = 0;
            }
         }
      }
   }
   ResetStats();
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Divide this histogram by h1.
///
/// `this = this/h1`
/// if errors are defined (see TH1::Sumw2), errors are also recalculated.
/// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
/// if not already set.
/// The resulting errors are calculated assuming uncorrelated histograms.
/// See the other TH1::Divide that gives the possibility to optionally
/// compute binomial errors.
///
/// IMPORTANT NOTE: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Scale
///
/// The function return kFALSE if the divide operation failed
 
Bool_t TH1::Divide(const TH1 *h1)
{
   if (!h1) {
      Error("Divide", "Input histogram passed does not exist (NULL).");
      return kFALSE;
   }
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   if (LoggedInconsistency("Divide", this, h1) >= kDifferentNumberOfBins) {
      return false;
   }
 
   //    Create Sumw2 if h1 has Sumw2 set
   if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
 
   //   - Loop on bins (including underflows/overflows)
   for (Int_t i = 0; i < fNcells; ++i) {
      Double_t c0 = RetrieveBinContent(i);
      Double_t c1 = h1->RetrieveBinContent(i);
      if (c1) UpdateBinContent(i, c0 / c1);
      else UpdateBinContent(i, 0);
 
      if(fSumw2.fN) {
         if (c1 == 0) { fSumw2.fArray[i] = 0; continue; }
         Double_t c1sq = c1 * c1;
         fSumw2.fArray[i] = (GetBinErrorSqUnchecked(i) * c1sq + h1->GetBinErrorSqUnchecked(i) * c0 * c0) / (c1sq * c1sq);
      }
   }
   ResetStats();
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Replace contents of this histogram by the division of h1 by h2.
///
/// `this = c1*h1/(c2*h2)`
///
/// If errors are defined (see TH1::Sumw2), errors are also recalculated
/// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
/// if not already set.
/// The resulting errors are calculated assuming uncorrelated histograms.
/// However, if option ="B" is specified, Binomial errors are computed.
/// In this case c1 and c2 do not make real sense and they are ignored.
///
/// IMPORTANT NOTE: If you intend to use the errors of this histogram later
/// you should call Sumw2 before making this operation.
/// This is particularly important if you fit the histogram after TH1::Divide
///
///  Please note also that in the binomial case errors are calculated using standard
///  binomial statistics, which means when b1 = b2, the error is zero.
///  If you prefer to have efficiency errors not going to zero when the efficiency is 1, you must
///  use the function TGraphAsymmErrors::BayesDivide, which will return an asymmetric and non-zero lower
///  error for the case b1=b2.
///
/// The function return kFALSE if the divide operation failed
 
Bool_t TH1::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
{
 
   TString opt = option;
   opt.ToLower();
   Bool_t binomial = kFALSE;
   if (opt.Contains("b")) binomial = kTRUE;
   if (!h1 || !h2) {
      Error("Divide", "At least one of the input histograms passed does not exist (NULL).");
      return kFALSE;
   }
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   if (LoggedInconsistency("Divide", this, h1) >= kDifferentNumberOfBins ||
       LoggedInconsistency("Divide", h1, h2) >= kDifferentNumberOfBins) {
      return false;
   }
 
   if (!c2) {
      Error("Divide","Coefficient of dividing histogram cannot be zero");
      return kFALSE;
   }
 
   //    Create Sumw2 if h1 or h2 have Sumw2 set, or if binomial errors are explicitly requested
   if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0 || binomial)) Sumw2();
 
   SetMinimum();
   SetMaximum();
 
   //   - Loop on bins (including underflows/overflows)
   for (Int_t i = 0; i < fNcells; ++i) {
      Double_t b1 = h1->RetrieveBinContent(i);
      Double_t b2 = h2->RetrieveBinContent(i);
      if (b2) UpdateBinContent(i, c1 * b1 / (c2 * b2));
      else UpdateBinContent(i, 0);
 
      if (fSumw2.fN) {
         if (b2 == 0) { fSumw2.fArray[i] = 0; continue; }
         Double_t b1sq = b1 * b1; Double_t b2sq = b2 * b2;
         Double_t c1sq = c1 * c1; Double_t c2sq = c2 * c2;
         Double_t e1sq = h1->GetBinErrorSqUnchecked(i);
         Double_t e2sq = h2->GetBinErrorSqUnchecked(i);
         if (binomial) {
            if (b1 != b2) {
               // in the case of binomial statistics c1 and c2 must be 1 otherwise it does not make sense
               // c1 and c2 are ignored
               //fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));//this is the formula in Hbook/Hoper1
               //fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/b2);     // old formula from G. Flucke
               // formula which works also for weighted histogram (see http://root-forum.cern.ch/viewtopic.php?t=3753 )
               fSumw2.fArray[i] = TMath::Abs( ( (1. - 2.* b1 / b2) * e1sq  + b1sq * e2sq / b2sq ) / b2sq );
            } else {
               //in case b1=b2 error is zero
               //use  TGraphAsymmErrors::BayesDivide for getting the asymmetric error not equal to zero
               fSumw2.fArray[i] = 0;
            }
         } else {
            fSumw2.fArray[i] = c1sq * c2sq * (e1sq * b2sq + e2sq * b1sq) / (c2sq * c2sq * b2sq * b2sq);
         }
      }
   }
   ResetStats();
   if (binomial)
      // in case of binomial division use denominator for number of entries
      SetEntries ( h2->GetEntries() );
 
   return kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Draw this histogram with options.
///
/// Histograms are drawn via the THistPainter class. Each histogram has
/// a pointer to its own painter (to be usable in a multithreaded program).
/// The same histogram can be drawn with different options in different pads.
/// When a histogram drawn in a pad is deleted, the histogram is
/// automatically removed from the pad or pads where it was drawn.
/// If a histogram is drawn in a pad, then filled again, the new status
/// of the histogram will be automatically shown in the pad next time
/// the pad is updated. One does not need to redraw the histogram.
/// To draw the current version of a histogram in a pad, one can use
/// `h->DrawCopy();`
/// This makes a clone of the histogram. Once the clone is drawn, the original
/// histogram may be modified or deleted without affecting the aspect of the
/// clone.
/// By default, TH1::Draw clears the current pad.
///
/// One can use TH1::SetMaximum and TH1::SetMinimum to force a particular
/// value for the maximum or the minimum scale on the plot.
///
/// TH1::UseCurrentStyle can be used to change all histogram graphics
/// attributes to correspond to the current selected style.
/// This function must be called for each histogram.
/// In case one reads and draws many histograms from a file, one can force
/// the histograms to inherit automatically the current graphics style
/// by calling before gROOT->ForceStyle();
///
/// See the THistPainter class for a description of all the drawing options.
 
void TH1::Draw(Option_t *option)
{
   TString opt1 = option; opt1.ToLower();
   TString opt2 = option;
   Int_t index  = opt1.Index("same");
 
   // Check if the string "same" is part of a TCutg name.
   if (index>=0) {
      Int_t indb = opt1.Index("[");
      if (indb>=0) {
         Int_t indk = opt1.Index("]");
         if (index>indb && index<indk) index = -1;
      }
   }
 
   // If there is no pad or an empty pad the "same" option is ignored.
   if (gPad) {
      if (!gPad->IsEditable()) gROOT->MakeDefCanvas();
      if (index>=0) {
         if (gPad->GetX1() == 0   && gPad->GetX2() == 1 &&
             gPad->GetY1() == 0   && gPad->GetY2() == 1 &&
             gPad->GetListOfPrimitives()->GetSize()==0) opt2.Remove(index,4);
      } else {
         //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();
      }
      gPad->IncrementPaletteColor(1, opt1);
   } else {
      if (index>=0) opt2.Remove(index,4);
   }
 
   AppendPad(opt2.Data());
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
/// By default a postfix "_copy" is added to the histogram name. Pass an empty postfix in case
/// you want to draw a histogram with the same name
///
/// See Draw for the list of options
 
TH1 *TH1::DrawCopy(Option_t *option, const char * name_postfix) const
{
   TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) gPad->Clear();
   TString newName;
   if (name_postfix) newName.Form("%s%s", GetName(), name_postfix);
   TH1 *newth1 = (TH1 *)Clone(newName.Data());
   newth1->SetDirectory(nullptr);
   newth1->SetBit(kCanDelete);
   if (gPad) gPad->IncrementPaletteColor(1, opt);
 
   newth1->AppendPad(option);
   return newth1;
}
 
////////////////////////////////////////////////////////////////////////////////
///  Draw a normalized copy of this histogram.
///
///  A clone of this histogram is normalized to norm and drawn with option.
///  A pointer to the normalized histogram is returned.
///  The contents of the histogram copy are scaled such that the new
///  sum of weights (excluding under and overflow) is equal to norm.
///  Note that the returned normalized histogram is not added to the list
///  of histograms in the current directory in memory.
///  It is the user's responsibility to delete this histogram.
///  The kCanDelete bit is set for the returned object. If a pad containing
///  this copy is cleared, the histogram will be automatically deleted.
///
///  See Draw for the list of options
 
TH1 *TH1::DrawNormalized(Option_t *option, Double_t norm) const
{
   Double_t sum = GetSumOfWeights();
   if (sum == 0) {
      Error("DrawNormalized","Sum of weights is null. Cannot normalize histogram: %s",GetName());
      return nullptr;
   }
   Bool_t addStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
   TH1 *h = (TH1*)Clone();
   h->SetBit(kCanDelete);
   // in case of drawing with error options - scale correctly the error
   TString opt(option); opt.ToUpper();
   if (fSumw2.fN == 0) {
      h->Sumw2();
      // do not use in this case the "Error option " for drawing which is enabled by default since the normalized histogram has now errors
      if (opt.IsNull() || opt == "SAME") opt += "HIST";
   }
   h->Scale(norm/sum);
   if (TMath::Abs(fMaximum+1111) > 1e-3) h->SetMaximum(fMaximum*norm/sum);
   if (TMath::Abs(fMinimum+1111) > 1e-3) h->SetMinimum(fMinimum*norm/sum);
   h->Draw(opt);
   TH1::AddDirectory(addStatus);
   return h;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Display a panel with all histogram drawing options.
///
///   See class TDrawPanelHist for example
 
void TH1::DrawPanel()
{
   if (!fPainter) {Draw(); if (gPad) gPad->Update();}
   if (fPainter) fPainter->DrawPanel();
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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 a value, distributed according to the Poisson
///    distribution, with f1 as the mean.
 
void TH1::Eval(TF1 *f1, Option_t *option)
{
   Double_t x[3];
   Int_t range, stat, add;
   if (!f1) return;
 
   TString opt = option;
   opt.ToLower();
   if (opt.Contains("a")) add   = 1;
   else                   add   = 0;
   if (opt.Contains("s")) stat  = 1;
   else                   stat  = 0;
   if (opt.Contains("r")) range = 1;
   else                   range = 0;
 
   // delete buffer if it is there since it will become invalid
   if (fBuffer) BufferEmpty(1);
 
   Int_t nbinsx  = fXaxis.GetNbins();
   Int_t nbinsy  = fYaxis.GetNbins();
   Int_t nbinsz  = fZaxis.GetNbins();
   if (!add) Reset();
 
   for (Int_t binz = 1; binz <= nbinsz; ++binz) {
      x[2]  = fZaxis.GetBinCenter(binz);
      for (Int_t biny = 1; biny <= nbinsy; ++biny) {
         x[1]  = fYaxis.GetBinCenter(biny);
         for (Int_t binx = 1; binx <= nbinsx; ++binx) {
            Int_t bin = GetBin(binx,biny,binz);
            x[0]  = fXaxis.GetBinCenter(binx);
            if (range && !f1->IsInside(x)) continue;
            Double_t fu = f1->Eval(x[0], x[1], x[2]);
            if (stat) fu = gRandom->PoissonD(fu);
            AddBinContent(bin, fu);
            if (fSumw2.fN) fSumw2.fArray[bin] += TMath::Abs(fu);
         }
      }
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
 
void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
   if (fPainter) fPainter->ExecuteEvent(event, px, py);
}
 
////////////////////////////////////////////////////////////////////////////////
/// This function allows to do discrete Fourier transforms of TH1 and TH2.
/// Available transform types and flags are described below.
///
/// To extract more information about the transform, use the function
/// TVirtualFFT::GetCurrentTransform() to get a pointer to the current
/// transform object.
///
/// \param[out] h_output histogram for the output. If a null pointer is passed, a new histogram is created
///          and returned, otherwise, the provided histogram is used and should be big enough
/// \param[in] option option parameters consists of 3 parts:
/// - option on what to return
///   - "RE" - returns a histogram of the real part of the output
///   - "IM" - returns a histogram of the imaginary part of the output
///   - "MAG"- returns a histogram of the magnitude of the output
///   - "PH" - returns a histogram of the phase of the output
/// - option of transform type
///   - "R2C"  - real to complex transforms - default
///   - "R2HC" - real to halfcomplex (special format of storing output data,
///     results the same as for R2C)
///   - "DHT" - discrete Hartley transform
///     real to real transforms (sine and cosine):
///   - "R2R_0", "R2R_1", "R2R_2", "R2R_3" - discrete cosine transforms of types I-IV
///   - "R2R_4", "R2R_5", "R2R_6", "R2R_7" - discrete sine transforms of types I-IV
///     To specify the type of each dimension of a 2-dimensional real to real
///     transform, use options of form "R2R_XX", for example, "R2R_02" for a transform,
///     which is of type "R2R_0" in 1st dimension and  "R2R_2" in the 2nd.
/// - option of transform flag
///   - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
///     performance
///   - "M" (from "measure")   - some time spend in finding the optimal way to do the transform
///   - "P" (from "patient")   - more time spend in finding the optimal way to do the transform
///   - "EX" (from "exhaustive") - the most optimal way is found
///     This option should be chosen depending on how many transforms of the same size and
///     type are going to be done. Planning is only done once, for the first transform of this
///     size and type. Default is "ES".
///
/// Examples of valid options: "Mag R2C M" "Re R2R_11" "Im R2C ES" "PH R2HC EX"
 
TH1* TH1::FFT(TH1* h_output, Option_t *option)
{
 
   Int_t ndim[3];
   ndim[0] = this->GetNbinsX();
   ndim[1] = this->GetNbinsY();
   ndim[2] = this->GetNbinsZ();
 
   TVirtualFFT *fft;
   TString opt = option;
   opt.ToUpper();
   if (!opt.Contains("2R")){
      if (!opt.Contains("2C") && !opt.Contains("2HC") && !opt.Contains("DHT")) {
         //no type specified, "R2C" by default
         opt.Append("R2C");
      }
      fft = TVirtualFFT::FFT(this->GetDimension(), ndim, opt.Data());
   }
   else {
      //find the kind of transform
      Int_t ind = opt.Index("R2R", 3);
      Int_t *kind = new Int_t[2];
      char t;
      t = opt[ind+4];
      kind[0] = atoi(&t);
      if (h_output->GetDimension()>1) {
         t = opt[ind+5];
         kind[1] = atoi(&t);
      }
      fft = TVirtualFFT::SineCosine(this->GetDimension(), ndim, kind, option);
      delete [] kind;
   }
 
   if (!fft) return nullptr;
   Int_t in=0;
   for (Int_t binx = 1; binx<=ndim[0]; binx++) {
      for (Int_t biny=1; biny<=ndim[1]; biny++) {
         for (Int_t binz=1; binz<=ndim[2]; binz++) {
            fft->SetPoint(in, this->GetBinContent(binx, biny, binz));
            in++;
         }
      }
   }
   fft->Transform();
   h_output = TransformHisto(fft, h_output, option);
   return h_output;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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 equal to or 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.
///
/// The function returns the corresponding bin number which has its content incremented by 1
 
Int_t TH1::Fill(Double_t x)
{
   if (fBuffer)  return BufferFill(x,1);
 
   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(x);
   if (bin <0) return -1;
   AddBinContent(bin);
   if (fSumw2.fN) ++fSumw2.fArray[bin];
   if (bin == 0 || bin > fXaxis.GetNbins()) {
      if (!GetStatOverflowsBehaviour()) return -1;
   }
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   return bin;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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 equal to or greater than the upper edge of last bin, the Overflow bin is incremented
///
/// If the weight is not equal to 1, the storage of the sum of squares of
/// weights is automatically triggered and the sum of the squares of weights is incremented
/// by \f$ w^2 \f$ in the bin corresponding to x.
///
/// The function returns the corresponding bin number which has its content incremented by w
 
Int_t TH1::Fill(Double_t x, Double_t w)
{
 
   if (fBuffer) return BufferFill(x,w);
 
   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(x);
   if (bin <0) return -1;
   if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW) )  Sumw2();   // must be called before AddBinContent
   if (fSumw2.fN)  fSumw2.fArray[bin] += w*w;
   AddBinContent(bin, w);
   if (bin == 0 || bin > fXaxis.GetNbins()) {
      if (!GetStatOverflowsBehaviour()) return -1;
   }
   Double_t z= w;
   fTsumw   += z;
   fTsumw2  += z*z;
   fTsumwx  += z*x;
   fTsumwx2 += z*x*x;
   return bin;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Increment bin with namex with a weight w
///
/// if x is less than the low-edge of the first bin, the Underflow bin is incremented
/// if x is equal to or greater than the upper edge of last bin, the Overflow bin is incremented
///
/// If the weight is not equal to 1, the storage of the sum of squares of
/// weights is automatically triggered and the sum of the squares of weights is incremented
/// by \f$ w^2 \f$ in the bin corresponding to x.
///
/// The function returns the corresponding bin number which has its content
/// incremented by w.
 
Int_t TH1::Fill(const char *namex, Double_t w)
{
   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(namex);
   if (bin <0) return -1;
   if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW))  Sumw2();
   if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
   AddBinContent(bin, w);
   if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
   Double_t z= w;
   fTsumw   += z;
   fTsumw2  += z*z;
   // this make sense if the histogram is not expanding (the x axis cannot be extended)
   if (!fXaxis.CanExtend() || !fXaxis.IsAlphanumeric()) {
      Double_t x = fXaxis.GetBinCenter(bin);
      fTsumwx  += z*x;
      fTsumwx2 += z*x*x;
   }
   return bin;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill this histogram with an array x and weights w.
///
/// \param[in] ntimes number of entries in arrays x and w (array size must be ntimes*stride)
/// \param[in] x array of values to be histogrammed
/// \param[in] w array of weighs
/// \param[in] stride step size through arrays x and w
///
/// If the weight is not equal to 1, the storage of the sum of squares of
/// weights is automatically triggered and the sum of the squares of weights is incremented
/// by \f$ w^2 \f$ in the bin corresponding to x.
/// if w is NULL each entry is assumed a weight=1
 
void TH1::FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
{
   //If a buffer is activated, fill buffer
   if (fBuffer) {
      ntimes *= stride;
      Int_t i = 0;
      for (i=0;i<ntimes;i+=stride) {
         if (!fBuffer) break;   // buffer can be deleted in BufferFill when is empty
         if (w) BufferFill(x[i],w[i]);
         else BufferFill(x[i], 1.);
      }
      // fill the remaining entries if the buffer has been deleted
      if (i < ntimes && !fBuffer) {
         auto weights = w ? &w[i] : nullptr;
         DoFillN((ntimes-i)/stride,&x[i],weights,stride);
      }
      return;
   }
   // call internal method
   DoFillN(ntimes, x, w, stride);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Internal method to fill histogram content from a vector
/// called directly by TH1::BufferEmpty
 
void TH1::DoFillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
{
   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 (bin <0) continue;
      if (w) ww = w[i];
      if (!fSumw2.fN && ww != 1.0 && !TestBit(TH1::kIsNotW))  Sumw2();
      if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
      AddBinContent(bin, ww);
      if (bin == 0 || bin > nbins) {
         if (!GetStatOverflowsBehaviour()) continue;
      }
      Double_t z= ww;
      fTsumw   += z;
      fTsumw2  += z*z;
      fTsumwx  += z*x[i];
      fTsumwx2 += z*x[i]*x[i];
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill histogram following distribution in function fname.
///
///  @param fname  : Function name used for filling the histogram
///  @param ntimes : number of times the histogram is filled
///  @param rng    : (optional) Random number generator used to sample
///
///
/// The distribution contained in the function fname (TF1) is integrated
/// over the channel contents for the bin range of this histogram.
/// It is normalized to 1.
///
/// Getting one random number implies:
///  - Generating a random number between 0 and 1 (say r1)
///  - Look in which bin in the normalized integral r1 corresponds to
///  - Fill histogram channel
///    ntimes random numbers are generated
///
/// One can also call TF1::GetRandom to get a random variate from a function.
 
void TH1::FillRandom(const char *fname, Int_t ntimes, TRandom * rng)
{
   Int_t bin, binx, ibin, loop;
   Double_t r1, x;
   //   - 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
 
   TAxis * xAxis = &fXaxis;
 
   // in case axis of histogram is not defined use the function axis
   if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
      Double_t xmin,xmax;
      f1->GetRange(xmin,xmax);
      Info("FillRandom","Using function axis and range [%g,%g]",xmin, xmax);
      xAxis = f1->GetHistogram()->GetXaxis();
   }
 
   Int_t first  = xAxis->GetFirst();
   Int_t last   = xAxis->GetLast();
   Int_t nbinsx = last-first+1;
 
   Double_t *integral = new Double_t[nbinsx+1];
   integral[0] = 0;
   for (binx=1;binx<=nbinsx;binx++) {
      Double_t fint = f1->Integral(xAxis->GetBinLowEdge(binx+first-1),xAxis->GetBinUpEdge(binx+first-1), 0.);
      integral[binx] = integral[binx-1] + fint;
   }
 
   //   - Normalize integral to 1
   if (integral[nbinsx] == 0 ) {
      delete [] integral;
      Error("FillRandom", "Integral = zero"); return;
   }
   for (bin=1;bin<=nbinsx;bin++)  integral[bin] /= integral[nbinsx];
 
   //   --------------Start main loop ntimes
   for (loop=0;loop<ntimes;loop++) {
      r1 = (rng) ? rng->Rndm() : gRandom->Rndm();
      ibin = TMath::BinarySearch(nbinsx,&integral[0],r1);
      //binx = 1 + ibin;
      //x    = xAxis->GetBinCenter(binx); //this is not OK when SetBuffer is used
      x    = xAxis->GetBinLowEdge(ibin+first)
             +xAxis->GetBinWidth(ibin+first)*(r1-integral[ibin])/(integral[ibin+1] - integral[ibin]);
      Fill(x);
   }
   delete [] integral;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill histogram following distribution in histogram h.
///
///  @param h      : Histogram  pointer used for sampling random number
///  @param ntimes : number of times the histogram is filled
///  @param rng    : (optional) Random number generator used for sampling
///
/// The distribution contained in the histogram h (TH1) is integrated
/// over the channel contents for the bin range of this histogram.
/// 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
///
/// SPECIAL CASE when the target histogram has the same binning as the source.
/// in this case we simply use a poisson distribution where
/// the mean value per bin = bincontent/integral.
 
void TH1::FillRandom(TH1 *h, Int_t ntimes, TRandom * rng)
{
   if (!h) { Error("FillRandom", "Null histogram"); return; }
   if (fDimension != h->GetDimension()) {
      Error("FillRandom", "Histograms with different dimensions"); return;
   }
   if (std::isnan(h->ComputeIntegral(true))) {
      Error("FillRandom", "Histograms contains negative bins, does not represent probabilities");
      return;
   }
 
   //in case the target histogram has the same binning and ntimes much greater
   //than the number of bins we can use a fast method
   Int_t first  = fXaxis.GetFirst();
   Int_t last   = fXaxis.GetLast();
   Int_t nbins = last-first+1;
   if (ntimes > 10*nbins) {
         auto inconsistency = CheckConsistency(this,h);
         if (inconsistency != kFullyConsistent) return; // do nothing
         Double_t sumw = h->Integral(first,last);
         if (sumw == 0) return;
         Double_t sumgen = 0;
         for (Int_t bin=first;bin<=last;bin++) {
            Double_t mean = h->RetrieveBinContent(bin)*ntimes/sumw;
            Double_t cont = (rng) ? rng->Poisson(mean) : gRandom->Poisson(mean);
            sumgen += cont;
            AddBinContent(bin,cont);
            if (fSumw2.fN) fSumw2.fArray[bin] += cont;
         }
 
         // fix for the fluctuations in the total number n
         // since we use Poisson instead of multinomial
         // add a correction to have ntimes as generated entries
         Int_t i;
         if (sumgen < ntimes) {
            // add missing entries
            for (i = Int_t(sumgen+0.5); i < ntimes; ++i)
            {
               Double_t x = h->GetRandom();
               Fill(x);
            }
         }
         else if (sumgen > ntimes) {
            // remove extra entries
            i =  Int_t(sumgen+0.5);
            while( i > ntimes) {
               Double_t x = h->GetRandom(rng);
               Int_t ibin = fXaxis.FindBin(x);
               Double_t y = RetrieveBinContent(ibin);
               // skip in case bin is empty
               if (y > 0) {
                  SetBinContent(ibin, y-1.);
                  i--;
               }
            }
         }
 
         ResetStats();
         return;
   }
   // case of different axis and not too large ntimes
 
   if (h->ComputeIntegral() ==0) return;
   Int_t loop;
   Double_t x;
   for (loop=0;loop<ntimes;loop++) {
      x = h->GetRandom();
      Fill(x);
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
/// This function tries to extend the axis if the given point belongs to an
/// under-/overflow bin AND if CanExtendAllAxes() is true.
///
/// See also TH1::GetBin, TAxis::FindBin and TAxis::FindFixBin
 
Int_t TH1::FindBin(Double_t x, Double_t y, Double_t z)
{
   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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
/// This function DOES NOT try to extend the axis if the given point belongs
/// to an under-/overflow bin.
///
/// See also TH1::GetBin, TAxis::FindBin and TAxis::FindFixBin
 
Int_t TH1::FindFixBin(Double_t x, Double_t y, Double_t z) const
{
   if (GetDimension() < 2) {
      return fXaxis.FindFixBin(x);
   }
   if (GetDimension() < 3) {
      Int_t nx   = fXaxis.GetNbins()+2;
      Int_t binx = fXaxis.FindFixBin(x);
      Int_t biny = fYaxis.FindFixBin(y);
      return  binx + nx*biny;
   }
   if (GetDimension() < 4) {
      Int_t nx   = fXaxis.GetNbins()+2;
      Int_t ny   = fYaxis.GetNbins()+2;
      Int_t binx = fXaxis.FindFixBin(x);
      Int_t biny = fYaxis.FindFixBin(y);
      Int_t binz = fZaxis.FindFixBin(z);
      return  binx + nx*(biny +ny*binz);
   }
   return -1;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Find first bin with content > threshold for axis (1=x, 2=y, 3=z)
/// if no bins with content > threshold is found the function returns -1.
/// The search will occur between the specified first and last bin. Specifying
/// the value of the last bin to search to less than zero will search until the
/// last defined bin.
 
Int_t TH1::FindFirstBinAbove(Double_t threshold, Int_t axis, Int_t firstBin, Int_t lastBin) const
{
   if (fBuffer) ((TH1*)this)->BufferEmpty();
 
   if (axis < 1 || (axis > 1 && GetDimension() == 1 ) ||
       ( axis > 2 && GetDimension() == 2 ) || ( axis > 3 && GetDimension() > 3 ) ) {
      Warning("FindFirstBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
      axis = 1;
   }
   if (firstBin < 1) {
      firstBin = 1;
   }
   Int_t nbinsx = fXaxis.GetNbins();
   Int_t nbinsy = (GetDimension() > 1 ) ? fYaxis.GetNbins() : 1;
   Int_t nbinsz = (GetDimension() > 2 ) ? fZaxis.GetNbins() : 1;
 
   if (axis == 1) {
      if (lastBin < 0 || lastBin > fXaxis.GetNbins()) {
         lastBin = fXaxis.GetNbins();
      }
      for (Int_t binx = firstBin; binx <= lastBin; binx++) {
         for (Int_t biny = 1; biny <= nbinsy; biny++) {
            for (Int_t binz = 1; binz <= nbinsz; binz++) {
               if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return binx;
            }
         }
      }
   }
   else if (axis == 2) {
      if (lastBin < 0 || lastBin > fYaxis.GetNbins()) {
         lastBin = fYaxis.GetNbins();
      }
      for (Int_t biny = firstBin; biny <= lastBin; biny++) {
         for (Int_t binx = 1; binx <= nbinsx; binx++) {
            for (Int_t binz = 1; binz <= nbinsz; binz++) {
               if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return biny;
           }
         }
      }
   }
   else if (axis == 3) {
      if (lastBin < 0 || lastBin > fZaxis.GetNbins()) {
         lastBin = fZaxis.GetNbins();
      }
      for (Int_t binz = firstBin; binz <= lastBin; binz++) {
         for (Int_t binx = 1; binx <= nbinsx; binx++) {
            for (Int_t biny = 1; biny <= nbinsy; biny++) {
               if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return binz;
            }
         }
      }
   }
 
   return -1;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Find last bin with content > threshold for axis (1=x, 2=y, 3=z)
/// if no bins with content > threshold is found the function returns -1.
/// The search will occur between the specified first and last bin. Specifying
/// the value of the last bin to search to less than zero will search until the
/// last defined bin.
 
Int_t TH1::FindLastBinAbove(Double_t threshold, Int_t axis, Int_t firstBin, Int_t lastBin) const
{
   if (fBuffer) ((TH1*)this)->BufferEmpty();
 
 
   if (axis < 1 || ( axis > 1 && GetDimension() == 1 ) ||
       ( axis > 2 && GetDimension() == 2 ) || ( axis > 3 && GetDimension() > 3) ) {
      Warning("FindFirstBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
      axis = 1;
   }
   if (firstBin < 1) {
      firstBin = 1;
   }
   Int_t nbinsx = fXaxis.GetNbins();
   Int_t nbinsy = (GetDimension() > 1 ) ? fYaxis.GetNbins() : 1;
   Int_t nbinsz = (GetDimension() > 2 ) ? fZaxis.GetNbins() : 1;
 
   if (axis == 1) {
      if (lastBin < 0 || lastBin > fXaxis.GetNbins()) {
         lastBin = fXaxis.GetNbins();
      }
      for (Int_t binx = lastBin; binx >= firstBin; binx--) {
         for (Int_t biny = 1; biny <= nbinsy; biny++) {
            for (Int_t binz = 1; binz <= nbinsz; binz++) {
               if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return binx;
            }
         }
      }
   }
   else if (axis == 2) {
      if (lastBin < 0 || lastBin > fYaxis.GetNbins()) {
         lastBin = fYaxis.GetNbins();
      }
      for (Int_t biny = lastBin; biny >= firstBin; biny--) {
         for (Int_t binx = 1; binx <= nbinsx; binx++) {
            for (Int_t binz = 1; binz <= nbinsz; binz++) {
               if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return biny;
            }
         }
      }
   }
   else if (axis == 3) {
      if (lastBin < 0 || lastBin > fZaxis.GetNbins()) {
         lastBin = fZaxis.GetNbins();
      }
      for (Int_t binz = lastBin; binz >= firstBin; binz--) {
         for (Int_t binx = 1; binx <= nbinsx; binx++) {
            for (Int_t biny = 1; biny <= nbinsy; biny++) {
               if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return binz;
            }
         }
      }
   }
 
   return -1;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Search object named name in the list of functions.
 
TObject *TH1::FindObject(const char *name) const
{
   if (fFunctions) return fFunctions->FindObject(name);
   return nullptr;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Search object obj in the list of functions.
 
TObject *TH1::FindObject(const TObject *obj) const
{
   if (fFunctions) return fFunctions->FindObject(obj);
   return nullptr;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fit histogram with function fname.
///
///
/// fname is the name of a function available in the global ROOT list of functions
/// `gROOT->GetListOfFunctions`
/// The list include any TF1 object created by the user plus some pre-defined functions
/// which are automatically created by ROOT the first time a pre-defined function is requested from `gROOT`
/// (i.e. when calling `gROOT->GetFunction(const char *name)`).
/// These pre-defined functions are:
///  - `gaus, gausn` where gausn is the normalized Gaussian
///  - `landau, landaun`
///  - `expo`
///  - `pol1,...9, chebyshev1,...9`.
///
/// For printing the list of all available functions do:
///
///       TF1::InitStandardFunctions();   // not needed if `gROOT->GetFunction` is called before
///       gROOT->GetListOfFunctions()->ls()
///
/// `fname` can also be a formula that is accepted by the linear fitter containing the special operator `++`,
/// representing linear components separated by `++` sign, for example `x++sin(x)` for fitting `[0]*x+[1]*sin(x)`
///
///  This function finds a pointer to the TF1 object with name `fname` and calls TH1::Fit(TF1 *, Option_t *, Option_t *,
///  Double_t, Double_t). See there for the fitting options and the details about fitting histograms
 
TFitResultPtr TH1::Fit(const char *fname ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
{
   char *linear;
   linear= (char*)strstr(fname, "++");
   Int_t ndim=GetDimension();
   if (linear){
      if (ndim<2){
         TF1 f1(fname, fname, xxmin, xxmax);
         return Fit(&f1,option,goption,xxmin,xxmax);
      }
      else if (ndim<3){
         TF2 f2(fname, fname);
         return Fit(&f2,option,goption,xxmin,xxmax);
      }
      else{
         TF3 f3(fname, fname);
         return Fit(&f3,option,goption,xxmin,xxmax);
      }
   }
   else{
      TF1 * f1 = (TF1*)gROOT->GetFunction(fname);
      if (!f1) { Printf("Unknown function: %s",fname); return -1; }
      return Fit(f1,option,goption,xxmin,xxmax);
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fit histogram with the function pointer f1.
///
/// \param[in] f1 pointer to the function object
/// \param[in] option string defining the fit options (see table below).
/// \param[in] goption specify a list of graphics options. See TH1::Draw for a complete list of these options.
/// \param[in] xxmin lower fitting range
/// \param[in] xxmax upper fitting range
/// \return A smart pointer to the TFitResult class
///
/// \anchor HFitOpt
/// ### Histogram Fitting Options
///
/// Here is the full list of fit options that can be given in the parameter `option`.
/// Several options can be used together by concatanating the strings without the need of any delimiters.
///
///  option | description
///  -------|------------
///   "L"  | Uses a log likelihood method (default is chi-square method). To be used when the histogram represents counts.
///   "WL" | Weighted log likelihood method. To be used when the histogram has been filled with weights different than 1. This is needed for getting correct parameter uncertainties for weighted fits.
///   "P"  | Uses Pearson chi-square method. Uses expected errors instead of the observed one (default case). The expected error is instead estimated from the square-root of the bin function value.
///   "MULTI" | Uses Loglikelihood method based on multi-nomial distribution. In this case the function must be normalized and one fits only the function shape.
///   "W"  | Fit using the chi-square method and ignoring the bin uncertainties and skip empty bins.
///   "WW" | Fit using the chi-square method and ignoring the bin uncertainties and include the empty bins.
///   "I"  | Uses the integral of function in the bin instead of the default bin center value.
///   "F"  | Uses the default minimizer (e.g. Minuit) when fitting a linear function (e.g. polN) instead of the linear fitter.
///   "U"  | Uses a user specified objective function (e.g. user providedlikelihood function) defined using `TVirtualFitter::SetFCN`
///   "E"  | Performs a better parameter errors estimation using the Minos technique for all fit parameters.
///   "M"  | Uses the IMPROVE algorithm (available only in TMinuit). This algorithm attempts improve the found local minimum by searching for a better one.
///   "S"  | The full result of the fit is returned in the `TFitResultPtr`. This is needed to get the covariance matrix of the fit. See `TFitResult` and the base class `ROOT::Math::FitResult`.
///   "Q"  | Quiet mode (minimum printing)
///   "V"  | Verbose mode (default is between Q and V)
///   "+"  | Adds this new fitted function to the list of fitted functions. By default, the previous function is deleted and only the last one is kept.
///   "N"  | Does not store the graphics function, does not draw the histogram with the function after fitting.
///   "0"  | Does not draw the histogram and the fitted function after fitting, but in contrast to option "N", it stores the fitted function in the histogram list of functions.
///   "R"  | Fit using a fitting range specified in the function range with `TF1::SetRange`.
///   "B"  | Use this option when you want to fix or set limits on one or more parameters and the fitting function is a predefined one (e.g gaus, expo,..), otherwise in case of pre-defined functions, some default initial values and limits will be used.
///   "C"  | In case of linear fitting, do no calculate the chisquare (saves CPU time).
///   "G"  | Uses the gradient implemented in `TF1::GradientPar` for the minimization. This allows to use Automatic Differentiation when it is supported by the provided TF1 function.
///   "WIDTH" | Scales the histogran bin content by the bin width (useful for variable bins histograms)
///   "SERIAL" | Runs in serial mode. By defult if ROOT is built with MT support and MT is enables, the fit is perfomed in multi-thread     - "E"  Perform better Errors estimation using Minos technique
///   "MULTITHREAD" | Forces usage of multi-thread execution whenever possible
///
/// The default fitting of an histogram (when no option is given) is perfomed as following:
///   - a chi-square fit (see below Chi-square Fits) computed using the bin histogram errors and excluding bins with zero errors (empty bins);
///   - the full range of the histogram is used;
///   - the default Minimizer with its default configuration is used (see below Minimizer Configuration) except for linear function;
///   - for linear functions (`polN`, `chenbyshev` or formula expressions combined using operator `++`) a linear minimization is used.
///   - only the status of the fit is returned;
///   - the fit is performed in Multithread whenever is enabled in ROOT;
///   - only the last fitted function is saved in the histogram;
///   - the histogram is drawn after fitting overalyed with the resulting fitting function
///
/// \anchor HFitMinimizer
/// ### Minimizer Configuration
///
/// The Fit is perfomed using  the default Minimizer, defined in the `ROOT::Math::MinimizerOptions` class.
/// It is possible to change the default minimizer and its configuration parameters by calling these static functions before fitting (before calling `TH1::Fit`):
///  - `ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerName, minimizerAgorithm)` for changing the minmizer and/or the corresponding algorithm.
///     For example `ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLMultiMin","BFGS");` will set the usage of the BFGS algorithm of the GSL multi-dimensional minimization
///     The current defaults are ("Minuit","Migrad").
///     See the documentation of the `ROOT::Math::MinimizerOptions` for the available minimizers in ROOT  and their corresponding algorithms.
///  - `ROOT::Math::MinimizerOptions::SetDefaultTolerance` for setting a different tolerance value for the minimization.
///  - `ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls` for setting the maximum number of function calls.
///  - `ROOT::Math::MinimizerOptions::SetDefaultPrintLevel` for changing the minimizer print level from level=0 (minimal printing) to level=3 maximum printing
///
/// Other options are possible depending on the Minimizer used, see the corresponding documentation.
/// The default minimizer can be also set in the resource file in etc/system.rootrc. For example
///
/// ~~~ {.cpp}
///     Root.Fitter:      Minuit2
/// ~~~
///
/// \anchor HFitChi2
/// ### Chi-square Fits
///
/// By default a chi-square (least-square) fit is performed on the histogram. The so-called modified least-square method
/// is used where the residual for each bin is computed using as error the observed value (the bin error) returned by `TH1::GetBinError`
///
/// \f[
///      Chi2 = \sum_{i}{ \left(\frac{y(i) - f(x(i) | p )}{e(i)} \right)^2 }
/// \f]
///
/// where `y(i)` is the bin content for each bin `i`, `x(i)` is the bin center and `e(i)` is the bin error (`sqrt(y(i)` for
/// an un-weighted histogram). Bins with zero errors are excluded from the fit. See also later the note on the treatment
/// of empty bins. When using option "I" the residual is computed not using the function value at the bin center, `f(x(i)|p)`,
/// but the integral of the function in the bin,   Integral{ f(x|p)dx }, divided by the bin volume.
/// When using option `P` (Pearson chi2), the expected error computed as `e(i) = sqrt(f(x(i)|p))` is used.
/// In this case empty bins are considered in the fit.
/// Both chi-square methods should not be used when the bin content represent counts, especially in case of low bin statistics,
/// because they could return a biased result.
///
/// \anchor HFitNLL
/// ### Likelihood Fits
///
/// When using option "L" a likelihood fit is used instead of the default chi-square fit.
/// The likelihood is built assuming a Poisson probability density function for each bin.
/// The negative log-likelihood to be minimized is
///
/// \f[
///       NLL = - \sum_{i}{ \log {\mathrm P} ( y(i) | f(x(i) | p ) ) }
/// \f]
/// where `P(y|f)` is the Poisson distribution of observing a count `y(i)` in the bin when the expected count is `f(x(i)|p)`.
/// The exact likelihood used is the Poisson likelihood described in this paper:
/// S. Baker and R. D. Cousins, “Clarification of the use of chi-square and likelihood functions in fits to histograms,”
/// Nucl. Instrum. Meth. 221 (1984) 437.
///
/// \f[
///       NLL = \sum_{i}{( f(x(i) | p ) + y(i)\log(y(i)/ f(x(i) | p )) - y(i)) }
/// \f]
/// By using this formulation, `2*NLL` can be interpreted as the chi-square resulting from the fit.
///
/// This method should be always used when the bin content represents counts (i.e. errors are sqrt(N) ).
/// The likelihood method has the advantage of treating correctly bins with low statistics. In case of high
/// statistics/bin the distribution of the bin content becomes a normal distribution and the likelihood and the chi2 fit
/// give the same result.
///
/// The likelihood method, although a bit slower, it is therefore the recommended method,
/// when the histogram represent counts (Poisson statistics), where the chi-square methods may
/// give incorrect results, especially in case of low statistics.
/// In case of a weighted histogram, it is possible to perform also a likelihood fit by using the
/// option "WL". Note a weighted histogram is a histogram which has been filled with weights and it
/// has the information on the sum of the weight square for each bin ( TH1::Sumw2() has been called).
/// The bin error for a weighted histogram is the square root of the sum of the weight square.
///
/// \anchor HFitRes
/// ### Fit Result
///
/// The function returns a TFitResultPtr which can hold a  pointer to a TFitResult object.
/// By default the TFitResultPtr contains only the status of the fit which is return by an
/// automatic conversion of the TFitResultPtr to an integer. One can write in this case directly:
///
/// ~~~ {.cpp}
///     Int_t fitStatus =  h->Fit(myFunc);
/// ~~~
///
/// If the option "S" is instead used, TFitResultPtr behaves as a smart
/// pointer to the TFitResult object. This is useful for retrieving the full result information from the fit, such as the covariance matrix,
/// as shown in this example code:
///
/// ~~~ {.cpp}
///     TFitResultPtr r = h->Fit(myFunc,"S");
///     TMatrixDSym cov = r->GetCovarianceMatrix();  //  to access the covariance matrix
///     Double_t chi2   = r->Chi2(); // to retrieve the fit chi2
///     Double_t par0   = r->Parameter(0); // retrieve the value for the parameter 0
///     Double_t err0   = r->ParError(0); // retrieve the error for the parameter 0
///     r->Print("V");     // print full information of fit including covariance matrix
///     r->Write();        // store the result in a file
/// ~~~
///
/// The fit parameters, error and chi-square (but not covariance matrix) can be retrieved also
/// directly from the fitted function that is passed to this call.
/// Given a pointer to an associated fitted function `myfunc`, one can retrieve the function/fit
/// parameters with calls such as:
///
/// ~~~ {.cpp}
///     Double_t chi2 = myfunc->GetChisquare();
///     Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
///     Double_t err0 = myfunc->GetParError(0);  //error on first parameter
/// ~~~
///
/// ##### Associated functions
///
/// One or more object ( can be added to the list
/// of functions (fFunctions) associated to each histogram.
/// When TH1::Fit is invoked, the fitted function is added to the histogram list of functions (fFunctions).
/// If the histogram is made persistent, the list of associated functions is also persistent.
/// Given a histogram h, one can retrieve an associated function with:
///
/// ~~~ {.cpp}
///  TF1 *myfunc = h->GetFunction("myfunc");
/// ~~~
/// or by quering directly the list obtained by calling `TH1::GetListOfFunctions`.
///
/// \anchor HFitStatus
/// ### Fit status
///
/// The status of the fit is obtained converting the TFitResultPtr to an integer
/// independently if the fit option "S" is used or not:
///
/// ~~~ {.cpp}
///     TFitResultPtr r = h->Fit(myFunc,opt);
///     Int_t fitStatus = r;
/// ~~~
///
/// - `status = 0` : the fit has been performed successfully (i.e no error occurred).
/// - `status < 0` : there is an error not connected with the minimization procedure, for example  when a wrong function is used.
/// - `status > 0` : return status from Minimizer, depends on used Minimizer. For example for TMinuit and Minuit2 we have:
///     - `status =  migradStatus + 10*minosStatus + 100*hesseStatus + 1000*improveStatus`.
///     TMinuit returns 0 (for migrad, minos, hesse or improve) in case of success and 4 in case of error (see the documentation of TMinuit::mnexcm). For example, for an error
///     only in Minos but not in Migrad a fitStatus of 40 will be returned.
///     Minuit2 returns 0 in case of success and different values in migrad,minos or
///     hesse depending on the error. See in this case the documentation of
///     Minuit2Minimizer::Minimize for the migrad return status, Minuit2Minimizer::GetMinosError for the
///     minos return status and Minuit2Minimizer::Hesse for the hesse return status.
///     If other minimizers are used see their specific documentation for the status code returned.
///     For example in the case of Fumili, see TFumili::Minimize.
///
/// \anchor HFitRange
/// ### Fitting in a range
///
/// In order to fit in a sub-range of the histogram you have two options:
/// -  pass to this function the lower (`xxmin`) and upper (`xxmax`) values for the fitting range;
/// -  define a specific range in the fitted function and use the fitting option "R".
///    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:
///
/// ~~~ {.cpp}
///      TF1 *f1 = new TF1("f1", "gaus", 1, 3);
///      histo->Fit("f1", "R");
/// ~~~
///
/// The fitting range is also limited by the histogram range defined using TAxis::SetRange
/// or TAxis::SetRangeUser. Therefore the fitting range is the smallest range between the
/// histogram one and the one defined by one of the two previous options described above.
///
/// \anchor HFitInitial
/// ### Setting initial conditions
///
/// Parameters must be initialized before invoking the Fit function.
/// The setting of the parameter initial values is automatic for the
/// predefined functions such as poln, expo, gaus, landau. One can however disable
/// this automatic computation by using the option "B".
/// Note that if a predefined function is defined with an argument,
/// eg, gaus(0), expo(1), you must specify the initial values for
/// the parameters.
/// You can specify boundary limits for some or all parameters via
///
/// ~~~ {.cpp}
///      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:
///
/// ~~~ {.cpp}
///      func->SetParameters(0, 3.1, 1.e-6, -8, 0, 100);
///      func->SetParLimits(3, -10, -4);
///      func->FixParameter(4, 0);
///      func->SetParLimits(5, 1, 1);
/// ~~~
///
/// With this setup, parameters 0->2 can vary freely
/// Parameter 3 has boundaries [-10,-4] with initial value -8
/// Parameter 4 is fixed to 0
/// Parameter 5 is fixed to 100.
/// When the lower limit and upper limit are equal, the parameter is fixed.
/// However to fix a parameter to 0, one must call the FixParameter function.
///
/// \anchor HFitStatBox
/// ### Fit Statistics Box
///
///   The statistics box can display the result of the fit.
///   You can change the statistics box to display the fit parameters with
///   the TStyle::SetOptFit(mode) method. This mode has four digits.
///   mode = pcev  (default = 0111)
///
///       v = 1;  print name/values of parameters
///       e = 1;  print errors (if e=1, v must be 1)
///       c = 1;  print Chisquare/Number of degrees of freedom
///       p = 1;  print Probability
///
///   For example: gStyle->SetOptFit(1011);
///   prints the fit probability, parameter names/values, and errors.
///   You can change the position of the statistics box with these lines
///   (where g is a pointer to the TGraph):
///
///       TPaveStats *st = (TPaveStats*)g->GetListOfFunctions()->FindObject("stats");
///       st->SetX1NDC(newx1); //new x start position
///       st->SetX2NDC(newx2); //new x end position
///
/// \anchor HFitExtra
/// ### Additional Notes on Fitting
///
/// #### Fitting a histogram of dimension N with a function of dimension N-1
///
/// It is possible to fit a TH2 with a TF1 or a TH3 with a TF2.
/// In this case the chi-square is computed from the squared error distance between the function values and the bin centers weighted by the bin content.
/// For correct error scaling, the obtained parameter error are corrected as in the case when the
/// option "W" is used.
///
/// #### User defined objective functions
///
/// By default when fitting a chi square function is used for fitting. When option "L" is used
/// a Poisson likelihood function is used. Using option "MULTI" a multinomial likelihood fit is used.
/// Thes functions are defined in the header Fit/Chi2Func.h or Fit/PoissonLikelihoodFCN and they
/// are implemented using the routines FitUtil::EvaluateChi2 or FitUtil::EvaluatePoissonLogL in
/// the file math/mathcore/src/FitUtil.cxx.
/// It is possible to specify a user defined fitting function, using option "U" and
/// calling the following functions:
///
/// ~~~ {.cpp}
///      TVirtualFitter::Fitter(myhist)->SetFCN(MyFittingFunction);
/// ~~~
///
/// where MyFittingFunction is of type:
///
/// ~~~ {.cpp}
///      extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
/// ~~~
///
/// #### Note on treatment of empty bins
///
/// Empty bins, which have the content equal to zero AND error equal to zero,
/// are excluded by default from the chi-square fit, but they are considered in the likelihood fit.
/// since they affect the likelihood if the function value in these bins is not negligible.
/// Note that if the histogram is having bins with zero content and non zero-errors they are considered as
/// any other bins in the fit. Instead bins with zero error and non-zero content are by default excluded in the chi-squared fit.
/// In general, one should not fit a histogram with non-empty bins and zero errors.
///
/// If the bin errors are not known, one should use the fit option "W", which gives a weight=1 for each bin (it is an unweighted least-square
/// fit). When using option "WW" the empty bins will be also considered in the chi-square fit with an error of 1.
/// Note that in this fitting case (option "W" or "WW") the resulting fitted parameter errors
/// are corrected by the obtained chi2 value using this scaling expression:
/// `errorp *= sqrt(chisquare/(ndf-1))` as it is done when fitting a TGraph with
/// no point errors.
///
/// #### Excluding points
///
/// You can use TF1::RejectPoint inside your fitting function to exclude some points
/// within a certain range from the fit. See the tutorial `fit/fitExclude.C`.
///
///
/// #### Warning when using the option "0"
///
/// When selecting the option "0", the fitted function is added to
/// the list of functions of the histogram, but it is not drawn when the histogram is drawn.
/// You can undo this behaviour resetting its corresponding bit in the TF1 object as following:
///
/// ~~~ {.cpp}
///     h.Fit("myFunction", "0"); // fit, store function but do not draw
///     h.Draw(); // function is not drawn
///     h.GetFunction("myFunction")->ResetBit(TF1::kNotDraw);
///     h.Draw();  // function is visible again
/// ~~~
///
 
TFitResultPtr TH1::Fit(TF1 *f1 ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
{
   // implementation of Fit method is in file hist/src/HFitImpl.cxx
   Foption_t fitOption;
   ROOT::Fit::FitOptionsMake(ROOT::Fit::EFitObjectType::kHistogram,option,fitOption);
 
   // create range and minimizer options with default values
   ROOT::Fit::DataRange range(xxmin,xxmax);
   ROOT::Math::MinimizerOptions minOption;
 
   // need to empty the buffer before
   // (t.b.d. do a ML unbinned fit with buffer data)
   if (fBuffer) BufferEmpty();
 
   return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Display a panel with all histogram fit options.
///
/// See class TFitPanel for example
 
void TH1::FitPanel()
{
   if (!gPad)
      gROOT->MakeDefCanvas();
 
   if (!gPad) {
      Error("FitPanel", "Unable to create a default canvas");
      return;
   }
 
 
   // use plugin manager to create instance of TFitEditor
   TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
   if (handler && handler->LoadPlugin() != -1) {
      if (handler->ExecPlugin(2, gPad, this) == 0)
         Error("FitPanel", "Unable to create the FitPanel");
   }
   else
         Error("FitPanel", "Unable to find the FitPanel plug-in");
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return a histogram containing the asymmetry of this histogram with h2,
/// where the asymmetry is defined as:
///
/// ~~~ {.cpp}
/// Asymmetry = (h1 - h2)/(h1 + h2)  where h1 = this
/// ~~~
///
/// works for 1D, 2D, etc. histograms
/// c2 is an optional argument that gives a relative weight between the two
/// histograms, and dc2 is the error on this weight.  This is useful, for example,
/// when forming an asymmetry between two histograms from 2 different data sets that
/// need to be normalized to each other in some way.  The function calculates
/// the errors assuming Poisson statistics on h1 and h2 (that is, dh = sqrt(h)).
///
/// example:  assuming 'h1' and 'h2' are already filled
///
/// ~~~ {.cpp}
/// h3 = h1->GetAsymmetry(h2)
/// ~~~
///
/// then 'h3' is created and filled with the asymmetry between 'h1' and 'h2';
/// h1 and h2 are left intact.
///
/// Note that it is the user's responsibility to manage the created histogram.
/// The name of the returned histogram will be `Asymmetry_nameOfh1-nameOfh2`
///
/// code proposed by Jason Seely (seely@mit.edu) and adapted by R.Brun
///
/// clone the histograms so top and bottom will have the
/// correct dimensions:
/// Sumw2 just makes sure the errors will be computed properly
/// when we form sums and ratios below.
 
TH1 *TH1::GetAsymmetry(TH1* h2, Double_t c2, Double_t dc2)
{
   TH1 *h1 = this;
   TString name =  TString::Format("Asymmetry_%s-%s",h1->GetName(),h2->GetName() );
   TH1 *asym   = (TH1*)Clone(name);
 
   // set also the title
   TString title = TString::Format("(%s - %s)/(%s+%s)",h1->GetName(),h2->GetName(),h1->GetName(),h2->GetName() );
   asym->SetTitle(title);
 
   asym->Sumw2();
   Bool_t addStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
   TH1 *top    = (TH1*)asym->Clone();
   TH1 *bottom = (TH1*)asym->Clone();
   TH1::AddDirectory(addStatus);
 
   // form the top and bottom of the asymmetry, and then divide:
   top->Add(h1,h2,1,-c2);
   bottom->Add(h1,h2,1,c2);
   asym->Divide(top,bottom);
 
   Int_t   xmax = asym->GetNbinsX();
   Int_t   ymax = asym->GetNbinsY();
   Int_t   zmax = asym->GetNbinsZ();
 
   if (h1->fBuffer) h1->BufferEmpty(1);
   if (h2->fBuffer) h2->BufferEmpty(1);
   if (bottom->fBuffer) bottom->BufferEmpty(1);
 
   // now loop over bins to calculate the correct errors
   // the reason this error calculation looks complex is because of c2
   for(Int_t i=1; i<= xmax; i++){
      for(Int_t j=1; j<= ymax; j++){
         for(Int_t k=1; k<= zmax; k++){
            Int_t bin = GetBin(i, j, k);
            // here some bin contents are written into variables to make the error
            // calculation a little more legible:
            Double_t a   = h1->RetrieveBinContent(bin);
            Double_t b   = h2->RetrieveBinContent(bin);
            Double_t bot = bottom->RetrieveBinContent(bin);
 
            // make sure there are some events, if not, then the errors are set = 0
            // automatically.
            //if(bot < 1){} was changed to the next line from recommendation of Jason Seely (28 Nov 2005)
            if(bot < 1e-6){}
            else{
               // computation of errors by Christos Leonidopoulos
               Double_t dasq  = h1->GetBinErrorSqUnchecked(bin);
               Double_t dbsq  = h2->GetBinErrorSqUnchecked(bin);
               Double_t error = 2*TMath::Sqrt(a*a*c2*c2*dbsq + c2*c2*b*b*dasq+a*a*b*b*dc2*dc2)/(bot*bot);
               asym->SetBinError(i,j,k,error);
            }
         }
      }
   }
   delete top;
   delete bottom;
 
   return asym;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Static function
/// return the default buffer size for automatic histograms
/// the parameter fgBufferSize may be changed via SetDefaultBufferSize
 
Int_t TH1::GetDefaultBufferSize()
{
   return fgBufferSize;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
/// see TH1::SetDefaultSumw2.
 
Bool_t TH1::GetDefaultSumw2()
{
   return fgDefaultSumw2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the current number of entries.
 
Double_t TH1::GetEntries() const
{
   if (fBuffer) {
      Int_t nentries = (Int_t) fBuffer[0];
      if (nentries > 0) return nentries;
   }
 
   return fEntries;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Number of effective entries of the histogram.
///
/// \f[
/// neff = \frac{(\sum Weights )^2}{(\sum Weight^2 )}
/// \f]
///
/// In case of an unweighted histogram this number is equivalent to the
/// number of entries of the histogram.
/// For a weighted histogram, this number corresponds to the hypothetical number of unweighted entries
/// a histogram would need to have the same statistical power as this weighted histogram.
/// Note: The underflow/overflow are included if one has set the TH1::StatOverFlows flag
/// and if the statistics has been computed at filling time.
/// If a range is set in the histogram the number is computed from the given range.
 
Double_t TH1::GetEffectiveEntries() const
{
   Stat_t s[kNstat];
   this->GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
   return (s[1] ? s[0]*s[0]/s[1] : TMath::Abs(s[0]) );
}
 
////////////////////////////////////////////////////////////////////////////////
/// Set highlight (enable/disable) mode for the histogram
/// by default highlight mode is disable
 
void TH1::SetHighlight(Bool_t set)
{
   if (IsHighlight() == set)
      return;
   if (fDimension > 2) {
      Info("SetHighlight", "Supported only 1-D or 2-D histograms");
      return;
   }
 
   SetBit(kIsHighlight, set);
 
   if (fPainter)
      fPainter->SetHighlight();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Redefines TObject::GetObjectInfo.
/// Displays the histogram info (bin number, contents, integral up to bin
/// corresponding to cursor position px,py
 
char *TH1::GetObjectInfo(Int_t px, Int_t py) const
{
   return ((TH1*)this)->GetPainter()->GetObjectInfo(px,py);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return pointer to painter.
/// If painter does not exist, it is created
 
TVirtualHistPainter *TH1::GetPainter(Option_t *option)
{
   if (!fPainter) {
      TString opt = option;
      opt.ToLower();
      if (opt.Contains("gl") || gStyle->GetCanvasPreferGL()) {
         //try to create TGLHistPainter
         TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TGLHistPainter");
 
         if (handler && handler->LoadPlugin() != -1)
            fPainter = reinterpret_cast<TVirtualHistPainter *>(handler->ExecPlugin(1, this));
      }
   }
 
   if (!fPainter) fPainter = TVirtualHistPainter::HistPainter(this);
 
   return fPainter;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute Quantiles for this histogram
/// Quantile x_q of a probability distribution Function F is defined as
///
/// ~~~ {.cpp}
///      F(x_q) = q with 0 <= q <= 1.
/// ~~~
///
/// For instance the median x_0.5 of a distribution is defined as that value
/// of the random variable for which the distribution function equals 0.5:
///
/// ~~~ {.cpp}
///      F(x_0.5) = Probability(x < x_0.5) = 0.5
/// ~~~
///
/// code from Eddy Offermann, Renaissance
///
/// \param[in] nprobSum maximum size of array q and size of array probSum (if given)
/// \param[in] probSum array of positions where quantiles will be computed.
///   - if probSum is null, probSum will be computed internally and will
///     have a size = number of bins + 1 in h. it will correspond to the
///     quantiles calculated at the lowest edge of the histogram (quantile=0) and
///     all the upper edges of the bins.
///   - if probSum is not null, it is assumed to contain at least nprobSum values.
/// \param[out] q array q filled with nq quantiles
/// \return value nq (<=nprobSum) with the number of quantiles computed
///
/// Note that the Integral of the histogram is automatically recomputed
/// if the number of entries is different of the number of entries when
/// the integral was computed last time. In case you do not use the Fill
/// functions to fill your histogram, but SetBinContent, you must call
/// TH1::ComputeIntegral before calling this function.
///
/// Getting quantiles q from two histograms and storing results in a TGraph,
/// a so-called QQ-plot
///
/// ~~~ {.cpp}
/// TGraph *gr = new TGraph(nprob);
/// h1->GetQuantiles(nprob,gr->GetX());
/// h2->GetQuantiles(nprob,gr->GetY());
/// gr->Draw("alp");
/// ~~~
///
/// Example:
///
/// ~~~ {.cpp}
/// void quantiles() {
///    // demo for quantiles
///    const Int_t nq = 20;
///    TH1F *h = new TH1F("h","demo quantiles",100,-3,3);
///    h->FillRandom("gaus",5000);
///
///    Double_t xq[nq];  // position where to compute the quantiles in [0,1]
///    Double_t yq[nq];  // array to contain the quantiles
///    for (Int_t i=0;i<nq;i++) xq[i] = Float_t(i+1)/nq;
///    h->GetQuantiles(nq,yq,xq);
///
///    //show the original histogram in the top pad
///    TCanvas *c1 = new TCanvas("c1","demo quantiles",10,10,700,900);
///    c1->Divide(1,2);
///    c1->cd(1);
///    h->Draw();
///
///    // show the quantiles in the bottom pad
///    c1->cd(2);
///    gPad->SetGrid();
///    TGraph *gr = new TGraph(nq,xq,yq);
///    gr->SetMarkerStyle(21);
///    gr->Draw("alp");
/// }
/// ~~~
 
Int_t TH1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
{
   if (GetDimension() > 1) {
      Error("GetQuantiles","Only available for 1-d histograms");
      return 0;
   }
 
   const Int_t nbins = GetXaxis()->GetNbins();
   if (!fIntegral) ComputeIntegral();
   if (fIntegral[nbins+1] != fEntries) ComputeIntegral();
 
   Int_t i, ibin;
   Double_t *prob = (Double_t*)probSum;
   Int_t nq = nprobSum;
   if (probSum == nullptr) {
      nq = nbins+1;
      prob = new Double_t[nq];
      prob[0] = 0;
      for (i=1;i<nq;i++) {
         prob[i] = fIntegral[i]/fIntegral[nbins];
      }
   }
 
   for (i = 0; i < nq; i++) {
      ibin = TMath::BinarySearch(nbins,fIntegral,prob[i]);
      while (ibin < nbins-1 && fIntegral[ibin+1] == prob[i]) {
         if (fIntegral[ibin+2] == prob[i]) ibin++;
         else break;
      }
      q[i] = GetBinLowEdge(ibin+1);
      const Double_t dint = fIntegral[ibin+1]-fIntegral[ibin];
      if (dint > 0) q[i] += GetBinWidth(ibin+1)*(prob[i]-fIntegral[ibin])/dint;
   }
 
   if (!probSum) delete [] prob;
   return nq;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Decode string choptin and fill fitOption structure.
 
Int_t TH1::FitOptionsMake(Option_t *choptin, Foption_t &fitOption)
{
   ROOT::Fit::FitOptionsMake(ROOT::Fit::EFitObjectType::kHistogram, choptin,fitOption);
   return 1;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute Initial values of parameters for a gaussian.
 
void H1InitGaus()
{
   Double_t allcha, sumx, sumx2, x, val, stddev, mean;
   Int_t bin;
   const Double_t sqrtpi = 2.506628;
 
   //   - Compute mean value and StdDev of the histogram in the given range
   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   TH1 *curHist = (TH1*)hFitter->GetObjectFit();
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   Double_t valmax  = curHist->GetBinContent(hxfirst);
   Double_t binwidx = curHist->GetBinWidth(hxfirst);
   allcha = sumx = sumx2 = 0;
   for (bin=hxfirst;bin<=hxlast;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;
   stddev  = sumx2/allcha - mean*mean;
   if (stddev > 0) stddev  = TMath::Sqrt(stddev);
   else         stddev  = 0;
   if (stddev == 0) stddev = binwidx*(hxlast-hxfirst+1)/4;
   //if the distribution is really gaussian, the best approximation
   //is binwidx*allcha/(sqrtpi*stddev)
   //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*stddev));
 
   //In case the mean value is outside the histo limits and
   //the StdDev is bigger than the range, we take
   //  mean = center of bins
   //  stddev  = half range
   Double_t xmin = curHist->GetXaxis()->GetXmin();
   Double_t xmax = curHist->GetXaxis()->GetXmax();
   if ((mean < xmin || mean > xmax) && stddev > (xmax-xmin)) {
      mean = 0.5*(xmax+xmin);
      stddev  = 0.5*(xmax-xmin);
   }
   TF1 *f1 = (TF1*)hFitter->GetUserFunc();
   f1->SetParameter(0,constant);
   f1->SetParameter(1,mean);
   f1->SetParameter(2,stddev);
   f1->SetParLimits(2,0,10*stddev);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute Initial values of parameters for an exponential.
 
void H1InitExpo()
{
   Double_t constant, slope;
   Int_t ifail;
   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   Int_t nchanx  = hxlast - hxfirst + 1;
 
   H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
 
   TF1 *f1 = (TF1*)hFitter->GetUserFunc();
   f1->SetParameter(0,constant);
   f1->SetParameter(1,slope);
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute Initial values of parameters for a polynom.
 
void H1InitPolynom()
{
   Double_t fitpar[25];
 
   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   TF1 *f1 = (TF1*)hFitter->GetUserFunc();
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   Int_t nchanx  = hxlast - hxfirst + 1;
   Int_t npar    = f1->GetNpar();
 
   if (nchanx <=1 || npar == 1) {
      TH1 *curHist = (TH1*)hFitter->GetObjectFit();
      fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
   } else {
      H1LeastSquareFit( nchanx, npar, fitpar);
   }
   for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Least squares lpolynomial fitting without weights.
///
/// \param[in] n number of points to fit
/// \param[in] m number of parameters
/// \param[in] a array of parameters
///
///  based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
///  (E.Keil.  revised by B.Schorr, 23.10.1981.)
 
void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a)
{
   const Double_t zero = 0.;
   const Double_t one = 1.;
   const Int_t idim = 20;
 
   Double_t  b[400]        /* was [20][20] */;
   Int_t i, k, l, ifail;
   Double_t power;
   Double_t da[20], xk, yk;
 
   if (m <= 2) {
      H1LeastSquareLinearFit(n, a[0], a[1], ifail);
      return;
   }
   if (m > idim || m > n) return;
   b[0]  = Double_t(n);
   da[0] = zero;
   for (l = 2; l <= m; ++l) {
      b[l-1]           = zero;
      b[m + l*20 - 21] = zero;
      da[l-1]          = zero;
   }
   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   TH1 *curHist  = (TH1*)hFitter->GetObjectFit();
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   for (k = hxfirst; k <= hxlast; ++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];
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.)
 
void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
{
   Double_t xbar, ybar, x2bar;
   Int_t i, n;
   Double_t xybar;
   Double_t fn, xk, yk;
   Double_t det;
 
   n     = TMath::Abs(ndata);
   ifail = -2;
   xbar  = ybar  = x2bar = xybar = 0;
   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   TH1 *curHist  = (TH1*)hFitter->GetObjectFit();
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   for (i = hxfirst; i <= hxlast; ++i) {
      xk = curHist->GetBinCenter(i);
      yk = curHist->GetBinContent(i);
      if (ndata < 0) {
         if (yk <= 0) yk = 1e-9;
         yk = TMath::Log(yk);
      }
      xbar  += xk;
      ybar  += yk;
      x2bar += xk*xk;
      xybar += xk*yk;
   }
   fn    = Double_t(n);
   det   = fn*x2bar - xbar*xbar;
   ifail = -1;
   if (det <= 0) {
      a0 = ybar/fn;
      a1 = 0;
      return;
   }
   ifail = 0;
   a0 = (x2bar*ybar - xbar*xybar) / det;
   a1 = (fn*xybar - xbar*ybar) / det;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Extracted from CERN Program library routine DSEQN.
///
/// Translated to C++ by Rene Brun
 
void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
{
   Int_t a_dim1, a_offset, b_dim1, b_offset;
   Int_t nmjp1, i, j, l;
   Int_t im1, jp1, nm1, nmi;
   Double_t s1, s21, s22;
   const Double_t one = 1.;
 
   /* Parameter adjustments */
   b_dim1 = idim;
   b_offset = b_dim1 + 1;
   b -= b_offset;
   a_dim1 = idim;
   a_offset = a_dim1 + 1;
   a -= a_offset;
 
   if (idim < n) return;
 
   ifail = 0;
   for (j = 1; j <= n; ++j) {
      if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
      a[j + j*a_dim1] = one / a[j + j*a_dim1];
      if (j == n) continue;
      jp1 = j + 1;
      for (l = jp1; l <= n; ++l) {
         a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
         s1 = -a[l + (j+1)*a_dim1];
         for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
         a[l + (j+1)*a_dim1] = -s1;
      }
   }
   if (k <= 0) return;
 
   for (l = 1; l <= k; ++l) {
      b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
   }
   if (n == 1) return;
   for (l = 1; l <= k; ++l) {
      for (i = 2; i <= n; ++i) {
         im1 = i - 1;
         s21 = -b[i + l*b_dim1];
         for (j = 1; j <= im1; ++j) {
            s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
         }
         b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
      }
      nm1 = n - 1;
      for (i = 1; i <= nm1; ++i) {
         nmi = n - i;
         s22 = -b[nmi + l*b_dim1];
         for (j = 1; j <= i; ++j) {
            nmjp1 = n - j + 1;
            s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
         }
         b[nmi + l*b_dim1] = -s22;
      }
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return Global bin number corresponding to binx,y,z.
///
/// 2-D and 3-D histograms are represented with a one dimensional
/// structure.
/// This has the advantage that all existing functions, such as
/// GetBinContent, GetBinError, GetBinFunction work for all dimensions.
///
/// In case of a TH1x, returns binx directly.
/// see TH1::GetBinXYZ for the inverse transformation.
///
/// Convention for numbering bins
///
/// For all histogram types: nbins, xlow, xup
///
///  - bin = 0;       underflow bin
///  - bin = 1;       first bin with low-edge xlow INCLUDED
///  - bin = nbins;   last bin with upper-edge xup EXCLUDED
///  - bin = nbins+1; overflow bin
///
/// In case of 2-D or 3-D histograms, a "global bin" number is defined.
/// For example, assuming a 3-D histogram with binx,biny,binz, the function
///
/// ~~~ {.cpp}
///     Int_t bin = h->GetBin(binx,biny,binz);
/// ~~~
///
/// returns a global/linearized bin number. This global bin is useful
/// to access the bin information independently of the dimension.
 
Int_t TH1::GetBin(Int_t binx, Int_t, Int_t) const
{
   Int_t ofx = fXaxis.GetNbins() + 1; // overflow bin
   if (binx < 0) binx = 0;
   if (binx > ofx) binx = ofx;
 
   return binx;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return binx, biny, binz corresponding to the global bin number globalbin
/// see TH1::GetBin function above
 
void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
{
   Int_t nx  = fXaxis.GetNbins()+2;
   Int_t ny  = fYaxis.GetNbins()+2;
 
   if (GetDimension() == 1) {
      binx = binglobal%nx;
      biny = 0;
      binz = 0;
      return;
   }
   if (GetDimension() == 2) {
      binx = binglobal%nx;
      biny = ((binglobal-binx)/nx)%ny;
      binz = 0;
      return;
   }
   if (GetDimension() == 3) {
      binx = binglobal%nx;
      biny = ((binglobal-binx)/nx)%ny;
      binz = ((binglobal-binx)/nx -biny)/ny;
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
///
/// @param rng (optional) Random number generator pointer used (default is gRandom)
///
/// The integral is automatically recomputed if the number of entries
/// is not the same then when the integral was computed.
/// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.
/// If the histogram has a bin with negative content a NaN is returned
 
Double_t TH1::GetRandom(TRandom * rng) const
{
   if (fDimension > 1) {
      Error("GetRandom","Function only valid for 1-d histograms");
      return 0;
   }
   Int_t nbinsx = GetNbinsX();
   Double_t integral = 0;
   // compute integral checking that all bins have positive content (see ROOT-5894)
   if (fIntegral) {
      if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true);
      else  integral = fIntegral[nbinsx];
   } else {
      integral = ((TH1*)this)->ComputeIntegral(true);
   }
   if (integral == 0) return 0;
   // return a NaN in case some bins have negative content
   if (integral == TMath::QuietNaN() ) return TMath::QuietNaN();
 
   Double_t r1 = (rng) ? rng->Rndm() : gRandom->Rndm();
   Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1);
   Double_t x = GetBinLowEdge(ibin+1);
   if (r1 > fIntegral[ibin]) x +=
      GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
   return x;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return content of bin number bin.
///
/// Implemented in TH1C,S,F,D
///
///  Convention for numbering bins
///
///  For all histogram types: nbins, xlow, xup
///
///  - bin = 0;       underflow bin
///  - bin = 1;       first bin with low-edge xlow INCLUDED
///  - bin = nbins;   last bin with upper-edge xup EXCLUDED
///  - bin = nbins+1; overflow bin
///
///  In case of 2-D or 3-D histograms, a "global bin" number is defined.
///  For example, assuming a 3-D histogram with binx,biny,binz, the function
///
/// ~~~ {.cpp}
///    Int_t bin = h->GetBin(binx,biny,binz);
/// ~~~
///
///  returns a global/linearized bin number. This global bin is useful
///  to access the bin information independently of the dimension.
 
Double_t TH1::GetBinContent(Int_t bin) const
{
   if (fBuffer) const_cast<TH1*>(this)->BufferEmpty();
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
 
   return RetrieveBinContent(bin);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute first binx in the range [firstx,lastx] for which
/// diff = abs(bin_content-c) <= maxdiff
///
/// In case several bins in the specified range with diff=0 are found
/// the first bin found is returned in binx.
/// In case several bins in the specified range satisfy diff <=maxdiff
/// the bin with the smallest difference is returned in binx.
/// In all cases the function returns the smallest difference.
///
/// NOTE1: if firstx <= 0, firstx is set to bin 1
///    if (lastx < firstx then firstx is set to the number of bins
///    ie if firstx=0 and lastx=0 (default) the search is on all bins.
///
/// NOTE2: if maxdiff=0 (default), the first bin with content=c is returned.
 
Double_t TH1::GetBinWithContent(Double_t c, Int_t &binx, Int_t firstx, Int_t lastx,Double_t maxdiff) const
{
   if (fDimension > 1) {
      binx = 0;
      Error("GetBinWithContent","function is only valid for 1-D histograms");
      return 0;
   }
 
   if (fBuffer) ((TH1*)this)->BufferEmpty();
 
   if (firstx <= 0) firstx = 1;
   if (lastx < firstx) lastx = fXaxis.GetNbins();
   Int_t binminx = 0;
   Double_t diff, curmax = 1.e240;
   for (Int_t i=firstx;i<=lastx;i++) {
      diff = TMath::Abs(RetrieveBinContent(i)-c);
      if (diff <= 0) {binx = i; return diff;}
      if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i;}
   }
   binx = binminx;
   return curmax;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Given a point x, approximates the value via linear interpolation
/// based on the two nearest bin centers
///
/// Andy Mastbaum 10/21/08
 
Double_t TH1::Interpolate(Double_t x) const
{
   if (fBuffer) ((TH1*)this)->BufferEmpty();
 
   Int_t xbin = fXaxis.FindFixBin(x);
   Double_t x0,x1,y0,y1;
 
   if(x<=GetBinCenter(1)) {
      return RetrieveBinContent(1);
   } else if(x>=GetBinCenter(GetNbinsX())) {
      return RetrieveBinContent(GetNbinsX());
   } else {
      if(x<=GetBinCenter(xbin)) {
         y0 = RetrieveBinContent(xbin-1);
         x0 = GetBinCenter(xbin-1);
         y1 = RetrieveBinContent(xbin);
         x1 = GetBinCenter(xbin);
      } else {
         y0 = RetrieveBinContent(xbin);
         x0 = GetBinCenter(xbin);
         y1 = RetrieveBinContent(xbin+1);
         x1 = GetBinCenter(xbin+1);
      }
      return y0 + (x-x0)*((y1-y0)/(x1-x0));
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// 2d Interpolation. Not yet implemented.
 
Double_t TH1::Interpolate(Double_t, Double_t) const
{
   Error("Interpolate","This function must be called with 1 argument for a TH1");
   return 0;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 3d Interpolation. Not yet implemented.
 
Double_t TH1::Interpolate(Double_t, Double_t, Double_t) const
{
   Error("Interpolate","This function must be called with 1 argument for a TH1");
   return 0;
}
 
///////////////////////////////////////////////////////////////////////////////
/// Check if a histogram is empty
///  (this is a protected method used mainly by TH1Merger )
 
Bool_t TH1::IsEmpty() const
{
   // if fTsumw or fentries are not zero histogram is not empty
   // need to use GetEntries() instead of fEntries in case of bugger histograms
   // so we will flash the buffer
   if (fTsumw != 0) return kFALSE;
   if (GetEntries() != 0) return kFALSE;
   // case fTSumw == 0 amd entries are also zero
   // this should not really happening, but if one sets content by hand
   // it can happen. a call to ResetStats() should be done in such cases
   double sumw = 0;
   for (int i = 0; i< GetNcells(); ++i) sumw += RetrieveBinContent(i);
   return (sumw != 0) ? kFALSE : kTRUE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return true if the bin is overflow.
 
Bool_t TH1::IsBinOverflow(Int_t bin, Int_t iaxis) const
{
   Int_t binx, biny, binz;
   GetBinXYZ(bin, binx, biny, binz);
 
   if (iaxis == 0) {
      if ( fDimension == 1 )
         return binx >= GetNbinsX() + 1;
      if ( fDimension == 2 )
         return (binx >= GetNbinsX() + 1) ||
            (biny >= GetNbinsY() + 1);
      if ( fDimension == 3 )
         return (binx >= GetNbinsX() + 1) ||
            (biny >= GetNbinsY() + 1) ||
            (binz >= GetNbinsZ() + 1);
      return kFALSE;
   }
   if (iaxis == 1)
      return binx >= GetNbinsX() + 1;
   if (iaxis == 2)
      return biny >= GetNbinsY() + 1;
   if (iaxis == 3)
      return binz >= GetNbinsZ() + 1;
 
   Error("IsBinOverflow","Invalid axis value");
   return kFALSE;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return