include/root/THnBase.h

0001 // @(#)root/hist:$Id$
0002 // Author: Axel Naumann (2011-12-20)
0003
0004 /*************************************************************************
0005  * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOT_THnBase
0013 #define ROOT_THnBase
0014
0015 /*************************************************************************
0016
0017  THnBase: Common base class for n-dimensional histogramming.
0018  Defines interfaces and algorithms.
0019
0020 *************************************************************************/
0021
0022
0023 #include "TNamed.h"
0024 #include "TMath.h"
0025 #include "TFitResultPtr.h"
0026 #include "TObjArray.h"
0027 #include "TArrayD.h"
0028 #include <vector>
0029
0030
0031 class TAxis;
0032 class TH1;
0033 class TH1D;
0034 class TH2D;
0035 class TH3D;
0036 class TF1;
0037 class THnIter;
0038
0039 namespace ROOT {
0040 namespace Internal {
0041    class THnBaseBinIter;
0042 }
0043 }
0044
0045 class THnBase: public TNamed {
0046 protected:
0047    Int_t      fNdimensions;  ///<  Number of dimensions
0048    TObjArray  fAxes;         ///<  Axes of the histogram
0049    TObjArray  fBrowsables;   ///<! Browser-helpers for each axis
0050    Double_t   fEntries;      ///<  Number of entries, spread over chunks
0051    Double_t   fTsumw;        ///<  Total sum of weights
0052    Double_t   fTsumw2;       ///<  Total sum of weights squared; -1 if no errors are calculated
0053    TArrayD    fTsumwx;       ///<  Total sum of weight*X for each dimension
0054    TArrayD    fTsumwx2;      ///<  Total sum of weight*X*X for each dimension
0055    std::vector<Double_t> fIntegral; ///<! vector with bin weight sums
0056    enum {
0057       kNoInt,
0058       kValidInt,
0059       kInvalidInt
0060    } fIntegralStatus;        ///<! status of integral
0061
0062  protected:
0063     THnBase() : fNdimensions(0), fEntries(0), fTsumw(0), fTsumw2(-1.), fIntegral(), fIntegralStatus(kNoInt) {}
0064
0065     THnBase(const char *name, const char *title, Int_t dim, const Int_t *nbins, const Double_t *xmin,
0066             const Double_t *xmax);
0067
0068     THnBase(const char* name, const char* title, const std::vector<TAxis>& axes);
0069
0070     THnBase(const char *name, const char *title, Int_t dim, const Int_t *nbins,
0071             const std::vector<std::vector<double>> &xbins);
0072
0073     THnBase(const THnBase &other);
0074
0075     THnBase &operator=(const THnBase &other);
0076
0077     THnBase(THnBase &&other);
0078
0079     THnBase &operator=(THnBase &&other);
0080
0081     void UpdateXStat(const Double_t *x, Double_t w = 1.)
0082     {
0083        if (GetCalculateErrors()) {
0084           for (Int_t d = 0; d < fNdimensions; ++d) {
0085              const Double_t xd = x[d];
0086              fTsumwx[d] += w * xd;
0087              fTsumwx2[d] += w * xd * xd;
0088           }
0089        }
0090     }
0091
0092    /// Increment the statistics due to filled weight "w",
0093    void FillBinBase(Double_t w) {
0094       fEntries += 1;
0095       if (GetCalculateErrors()) {
0096          fTsumw += w;
0097          fTsumw2 += w*w;
0098       }
0099       fIntegralStatus = kInvalidInt;
0100    }
0101
0102    virtual void InitStorage(Int_t* nbins, Int_t chunkSize) = 0;
0103    void Init(const char* name, const char* title,
0104              const TObjArray* axes, Bool_t keepTargetAxis,
0105              Int_t chunkSize = 1024 * 16);
0106    THnBase* CloneEmpty(const char* name, const char* title,
0107                        const TObjArray* axes, Bool_t keepTargetAxis) const;
0108    virtual void Reserve(Long64_t /*nbins*/) {}
0109    virtual void SetFilledBins(Long64_t /*nbins*/) {};
0110
0111    Bool_t CheckConsistency(const THnBase *h, const char *tag) const;
0112    TH1* CreateHist(const char* name, const char* title,
0113                    const TObjArray* axes, Bool_t keepTargetAxis) const;
0114    TObject* ProjectionAny(Int_t ndim, const Int_t* dim,
0115                           Bool_t wantNDim, Option_t* option = "") const;
0116    Bool_t PrintBin(Long64_t idx, Int_t* coord, Option_t* options) const;
0117    void AddInternal(const THnBase* h, Double_t c, Bool_t rebinned);
0118    THnBase* RebinBase(Int_t group) const;
0119    THnBase* RebinBase(const Int_t* group) const;
0120    void ResetBase(Option_t *option= "");
0121    void SetTitleImpl(const char *title, bool overrideAxesTitle);
0122
0123    static THnBase* CreateHnAny(const char* name, const char* title,
0124                                const TH1* h1, Bool_t sparse,
0125                                Int_t chunkSize = 1024 * 16);
0126    static THnBase* CreateHnAny(const char* name, const char* title,
0127                                const THnBase* hn, Bool_t sparse,
0128                                Int_t chunkSize = 1024 * 16);
0129
0130  public:
0131    ~THnBase() override;
0132
0133    TObjArray* GetListOfAxes() { return &fAxes; }
0134    const TObjArray* GetListOfAxes() const { return &fAxes; }
0135    TAxis* GetAxis(Int_t dim) const { return (TAxis*)fAxes[dim]; }
0136
0137    TFitResultPtr Fit(TF1 *f1 ,Option_t *option = "", Option_t *goption = "");
0138    TList* GetListOfFunctions() { return nullptr; }
0139
0140    virtual ROOT::Internal::THnBaseBinIter* CreateIter(Bool_t respectAxisRange) const = 0;
0141
0142    virtual Long64_t GetNbins() const = 0;
0143    Double_t GetEntries() const { return fEntries; }
0144    Double_t GetWeightSum() const { return fTsumw; }
0145    Int_t    GetNdimensions() const { return fNdimensions; }
0146    Bool_t   GetCalculateErrors() const { return fTsumw2 >= 0.; }
0147
0148    /// Calculate errors (or not if "calc" == kFALSE)
0149    void     CalculateErrors(Bool_t calc = kTRUE) {
0150       if (calc) Sumw2();
0151       else fTsumw2 = -1.;
0152    }
0153
0154    Long64_t Fill(const Double_t *x, Double_t w = 1.) {
0155       UpdateXStat(x, w);
0156       Long64_t bin = GetBin(x, kTRUE /*alloc*/);
0157       FillBin(bin, w);
0158       return bin;
0159    }
0160    Long64_t Fill(const char* name[], Double_t w = 1.) {
0161       Long64_t bin = GetBin(name, kTRUE /*alloc*/);
0162       FillBin(bin, w);
0163       return bin;
0164    }
0165
0166    /// Fill with the provided variadic arguments.
0167    /// The number of arguments must be equal to the number of histogram dimensions or, for weighted fills, to the
0168    /// number of dimensions + 1; in the latter case, the last function argument is used as weight.
0169    /// A separate `firstval` argument is needed so the compiler does not pick this overload instead of the non-templated
0170    /// Fill overloads
0171    template <typename... MoreTypes>
0172    Long64_t Fill(Double_t firstval, MoreTypes... morevals)
0173    {
0174       const std::array<double, 1 + sizeof...(morevals)> x{firstval, static_cast<double>(morevals)...};
0175       if (Int_t(x.size()) == GetNdimensions()) {
0176          // without weight
0177          return Fill(x.data());
0178       } else if (Int_t(x.size()) == (GetNdimensions() + 1)) {
0179          // with weight
0180          return Fill(x.data(), x.back());
0181       } else {
0182          Error("Fill", "Wrong number of arguments for number of histogram axes.");
0183       }
0184
0185       return -1;
0186    }
0187
0188    virtual void FillBin(Long64_t bin, Double_t w) = 0;
0189
0190    void SetBinEdges(Int_t idim, const Double_t* bins);
0191    Bool_t IsInRange(Int_t *coord) const;
0192    Double_t GetBinError(const Int_t *idx) const { return GetBinError(GetBin(idx)); }
0193    Double_t GetBinError(Long64_t linidx) const { return TMath::Sqrt(GetBinError2(linidx)); }
0194    void SetBinError(const Int_t* idx, Double_t e) { SetBinError(GetBin(idx), e); }
0195    void SetBinError(Long64_t bin, Double_t e) { SetBinError2(bin, e*e); }
0196    void AddBinContent(const Int_t* x, Double_t v = 1.) { AddBinContent(GetBin(x), v); }
0197    void SetEntries(Double_t entries) { fEntries = entries; }
0198    void SetTitle(const char *title) override;
0199
0200    std::vector<Double_t> GetBinCenter(const std::vector<Int_t> &idx) const;
0201
0202    Double_t GetBinContent(const Int_t *idx) const { return GetBinContent(GetBin(idx)); } // intentionally non-virtual
0203    virtual Double_t GetBinContent(Long64_t bin, Int_t* idx = nullptr) const = 0;
0204    virtual Double_t GetBinError2(Long64_t linidx) const = 0;
0205    virtual Long64_t GetBin(const Int_t* idx) const = 0;
0206    virtual Long64_t GetBin(const Double_t* x) const = 0;
0207    virtual Long64_t GetBin(const char* name[]) const = 0;
0208    virtual Long64_t GetBin(const Int_t* idx, Bool_t /*allocate*/ = kTRUE) = 0;
0209    virtual Long64_t GetBin(const Double_t* x, Bool_t /*allocate*/ = kTRUE) = 0;
0210    virtual Long64_t GetBin(const char* name[], Bool_t /*allocate*/ = kTRUE) = 0;
0211
0212    void SetBinContent(const Int_t* idx, Double_t v) { SetBinContent(GetBin(idx), v); } // intentionally non-virtual
0213    virtual void SetBinContent(Long64_t bin, Double_t v) = 0;
0214    virtual void SetBinError2(Long64_t bin, Double_t e2) = 0;
0215    virtual void AddBinError2(Long64_t bin, Double_t e2) = 0;
0216    virtual void AddBinContent(Long64_t bin, Double_t v = 1.) = 0;
0217
0218    Double_t GetSumw() const  { return fTsumw; }
0219    Double_t GetSumw2() const { return fTsumw2; }
0220    Double_t GetSumwx(Int_t dim) const  { return fTsumwx[dim]; }
0221    Double_t GetSumwx2(Int_t dim) const { return fTsumwx2[dim]; }
0222
0223    /// Project all bins into a 1-dimensional histogram,
0224    /// keeping only axis "xDim".
0225    /// If "option" contains:
0226    ///  - "E" errors will be calculated.
0227    ///  - "A" ranges of the taget axes will be ignored.
0228    ///  - "O" original axis range of the taget axes will be
0229    ///    kept, but only bins inside the selected range
0230    ///    will be filled.
0231    TH1D*    Projection(Int_t xDim, Option_t* option = "") const {
0232       return (TH1D*) ProjectionAny(1, &xDim, false, option);
0233    }
0234
0235    /// Project all bins into a 2-dimensional histogram,
0236    /// keeping only axes "xDim" and "yDim".
0237    ///
0238    /// WARNING: just like TH3::Project3D("yx") and TTree::Draw("y:x"),
0239    /// Projection(y,x) uses the first argument to define the y-axis and the
0240    /// second for the x-axis!
0241    ///
0242    /// If "option" contains "E" errors will be calculated.
0243    ///                      "A" ranges of the taget axes will be ignored.
0244    TH2D*    Projection(Int_t yDim, Int_t xDim, Option_t* option = "") const {
0245       const Int_t dim[2] = {xDim, yDim};
0246       return (TH2D*) ProjectionAny(2, dim, false, option);
0247    }
0248
0249    /// Project all bins into a 3-dimensional histogram,
0250    /// keeping only axes "xDim", "yDim", and "zDim".
0251    /// If "option" contains:
0252    ///  - "E" errors will be calculated.
0253    ///  - "A" ranges of the taget axes will be ignored.
0254    ///  - "O" original axis range of the taget axes will be
0255    ///    kept, but only bins inside the selected range
0256    ///    will be filled.
0257    TH3D*    Projection(Int_t xDim, Int_t yDim, Int_t zDim, Option_t* option = "") const {
0258       const Int_t dim[3] = {xDim, yDim, zDim};
0259       return (TH3D*) ProjectionAny(3, dim, false, option);
0260    }
0261
0262    /// @brief Project histogram onto the axes specified in dim.
0263    /// @param ndim Size of the dim array.
0264    /// @param dim Indices of the axes to project on.
0265    /// @param option - "E" errors will be calculated.
0266    ///               - "A" ranges of the target axes will be ignored.
0267    ///               - "O" original axis range of the target axes will be kept, but only bins inside the selected range
0268    ///               will be filled.
0269    /// @return A new THnBase derived class, such as THN or THNSparse.
0270    THnBase* ProjectionND(Int_t ndim, const Int_t* dim,
0271                          Option_t* option = "") const {
0272       return (THnBase*)ProjectionAny(ndim, dim, kTRUE /*wantNDim*/, option);
0273    }
0274
0275    Long64_t   Merge(TCollection* list);
0276
0277    void Scale(Double_t c);
0278    void Add(const THnBase* h, Double_t c=1.);
0279    void Add(const TH1* hist, Double_t c=1.);
0280    void Multiply(const THnBase* h);
0281    void Multiply(TF1* f, Double_t c = 1.);
0282    void Divide(const THnBase* h);
0283    void Divide(const THnBase* h1, const THnBase* h2, Double_t c1 = 1., Double_t c2 = 1., Option_t* option="");
0284    void RebinnedAdd(const THnBase* h, Double_t c=1.);
0285
0286    virtual void Reset(Option_t* option = "") = 0;
0287    virtual void Sumw2() = 0;
0288
0289    Double_t ComputeIntegral();
0290    void GetRandom(Double_t *rand, Bool_t subBinRandom = kTRUE);
0291    Double_t Integral(Bool_t respectAxisRange) const;
0292
0293    void Print(Option_t* option = "") const override;
0294    void PrintEntries(Long64_t from = 0, Long64_t howmany = -1, Option_t* options = nullptr) const;
0295    void PrintBin(Int_t* coord, Option_t* options) const {
0296       PrintBin(-1, coord, options);
0297    }
0298    void PrintBin(Long64_t idx, Option_t* options) const;
0299
0300    void Browse(TBrowser *b) override;
0301    Bool_t IsFolder() const override { return kTRUE; }
0302
0303    //void Draw(Option_t* option = "");
0304
0305    ClassDefOverride(THnBase, 1); // Common base for n-dimensional histogram
0306
0307    friend class THnIter;
0308 };
0309
0310 namespace ROOT {
0311 namespace Internal {
0312    // Helper class for browsing THnBase objects
0313    class THnBaseBrowsable: public TNamed {
0314    public:
0315       THnBaseBrowsable(THnBase* hist, Int_t axis);
0316       ~THnBaseBrowsable() override;
0317       void Browse(TBrowser *b) override;
0318       Bool_t IsFolder() const override { return kFALSE; }
0319
0320    private:
0321       THnBase* fHist; // Original histogram
0322       Int_t    fAxis; // Axis to visualize
0323       TH1*     fProj; // Projection result
0324       ClassDefOverride(THnBaseBrowsable, 0); // Browser-helper for THnBase
0325    };
0326
0327    // Base class for iterating over THnBase bins
0328    class THnBaseBinIter {
0329    public:
0330       THnBaseBinIter(Bool_t respectAxisRange):
0331          fRespectAxisRange(respectAxisRange), fHaveSkippedBin(kFALSE) {}
0332       virtual ~THnBaseBinIter();
0333       Bool_t HaveSkippedBin() const { return fHaveSkippedBin; }
0334       Bool_t RespectsAxisRange() const { return fRespectAxisRange; }
0335
0336       virtual Int_t GetCoord(Int_t dim) const = 0;
0337       virtual Long64_t Next(Int_t* coord = nullptr) = 0;
0338
0339    protected:
0340       Bool_t fRespectAxisRange;
0341       Bool_t fHaveSkippedBin;
0342    };
0343 }
0344 }
0345
0346 class THnIter: public TObject {
0347 public:
0348    THnIter(const THnBase* hist, Bool_t respectAxisRange = kFALSE):
0349       fIter(hist->CreateIter(respectAxisRange)) {}
0350    ~THnIter() override;
0351
0352    /// Return the next bin's index.
0353    /// If provided, set coord to that bin's coordinates (bin indexes).
0354    /// I.e. coord must point to Int_t[hist->GetNdimensions()]
0355    /// Returns -1 when all bins have been visited.
0356    Long64_t Next(Int_t* coord = nullptr) {
0357       return fIter->Next(coord);
0358    }
0359
0360    Int_t GetCoord(Int_t dim) const { return fIter->GetCoord(dim); }
0361    Bool_t HaveSkippedBin() const { return fIter->HaveSkippedBin(); }
0362    Bool_t RespectsAxisRange() const { return fIter->RespectsAxisRange(); }
0363
0364 private:
0365    ROOT::Internal::THnBaseBinIter* fIter;
0366    ClassDefOverride(THnIter, 0); //Iterator over bins of a THnBase.
0367 };
0368
0369 #endif //  ROOT_THnBase