diff --git a/hist/hist/inc/HFitInterface.h b/hist/hist/inc/HFitInterface.h
index f917d0e4b84bdbe527cbff58547956be5620f1ff..de941065ea15e95dc625dcceeb3eb2b731e703b4 100644
--- a/hist/hist/inc/HFitInterface.h
+++ b/hist/hist/inc/HFitInterface.h
@@ -19,7 +19,7 @@
class TH1;
-class THnSparse;
+class THnBase;
class TF1;
class TF2;
class TGraph;
@@ -73,9 +73,9 @@ namespace ROOT {
TFitResultPtr FitObject(TGraph2D * gr, TF1 *f1 , Foption_t & option , const ROOT::Math::MinimizerOptions & moption, const char *goption, ROOT::Fit::DataRange & range);
/**
- fitting function for a THnSparse (called from THnSparse::Fit)
+ fitting function for a THn / THnSparse (called from THnBase::Fit)
*/
- TFitResultPtr FitObject(THnSparse * s1, TF1 *f1, Foption_t & option, const ROOT::Math::MinimizerOptions & moption, const char *goption, ROOT::Fit::DataRange & range);
+ TFitResultPtr FitObject(THnBase * s1, TF1 *f1, Foption_t & option, const ROOT::Math::MinimizerOptions & moption, const char *goption, ROOT::Fit::DataRange & range);
#endif
/**
@@ -100,16 +100,16 @@ namespace ROOT {
void FillData ( SparseData & dv, const TH1 * hist, TF1 * func = 0);
/**
- fill the data vector from a THnSparse. Pass also the TF1 function which is
+ fill the data vector from a THnBase. Pass also the TF1 function which is
needed in case of integral option and to reject points rejected by the function
*/
- void FillData ( SparseData & dv, const THnSparse * hist, TF1 * func = 0);
+ void FillData ( SparseData & dv, const THnBase * hist, TF1 * func = 0);
/**
- fill the data vector from a THnSparse. Pass also the TF1 function which is
+ fill the data vector from a THnBase. Pass also the TF1 function which is
needed in case of integral option and to reject points rejected by the function
*/
- void FillData ( BinData & dv, const THnSparse * hist, TF1 * func = 0);
+ void FillData ( BinData & dv, const THnBase * hist, TF1 * func = 0);
/**
fill the data vector from a TGraph2D. Pass also the TF1 function which is
diff --git a/hist/hist/inc/LinkDef.h b/hist/hist/inc/LinkDef.h
index 4fc8beda7fef07e44eac0c2f0f3f07722861ca99..cf820e8c60143c1227c57e587ecb61281f415e75 100644
--- a/hist/hist/inc/LinkDef.h
+++ b/hist/hist/inc/LinkDef.h
@@ -59,6 +59,62 @@
#pragma link C++ class TH3S-;
#pragma link C++ class TH3I+;
#pragma link C++ class THLimitsFinder+;
+#pragma link C++ class THnBase+;
+#pragma link C++ class THnIter+;
+#pragma link C++ class TNDArray+;
+#pragma link C++ class TNDArrayT<Float_t>+;
+//#pragma link C++ class TNDArrayT<Float16_t>+;
+#pragma link C++ class TNDArrayT<Double_t>+;
+//#pragma link C++ class TNDArrayT<Double32_t>+;
+#pragma link C++ class TNDArrayT<Long64_t>+;
+#pragma link C++ class TNDArrayT<Long_t>+;
+#pragma link C++ class TNDArrayT<Int_t>+;
+#pragma link C++ class TNDArrayT<Short_t>+;
+#pragma link C++ class TNDArrayT<Char_t>+;
+#pragma link C++ class TNDArrayT<ULong64_t>+;
+#pragma link C++ class TNDArrayT<ULong_t>+;
+#pragma link C++ class TNDArrayT<UInt_t>+;
+#pragma link C++ class TNDArrayT<UShort_t>+;
+#pragma link C++ class TNDArrayRef<Float_t>+;
+//#pragma link C++ class TNDArrayRef<Float16_t>+;
+#pragma link C++ class TNDArrayRef<Double_t>+;
+//#pragma link C++ class TNDArrayRef<Double32_t>+;
+#pragma link C++ class TNDArrayRef<Long64_t>+;
+#pragma link C++ class TNDArrayRef<Long_t>+;
+#pragma link C++ class TNDArrayRef<Int_t>+;
+#pragma link C++ class TNDArrayRef<Short_t>+;
+#pragma link C++ class TNDArrayRef<Char_t>+;
+#pragma link C++ class TNDArrayRef<ULong64_t>+;
+#pragma link C++ class TNDArrayRef<ULong_t>+;
+#pragma link C++ class TNDArrayRef<UInt_t>+;
+#pragma link C++ class TNDArrayRef<UShort_t>+;
+#pragma link C++ class TNDArrayRef<const Float_t>+;
+//#pragma link C++ class TNDArrayRef<const Float16_t>+;
+#pragma link C++ class TNDArrayRef<const Double_t>+;
+//#pragma link C++ class TNDArrayRef<const Double32_t>+;
+#pragma link C++ class TNDArrayRef<const Long64_t>+;
+#pragma link C++ class TNDArrayRef<const Long_t>+;
+#pragma link C++ class TNDArrayRef<const Int_t>+;
+#pragma link C++ class TNDArrayRef<const Short_t>+;
+#pragma link C++ class TNDArrayRef<const Char_t>+;
+#pragma link C++ class TNDArrayRef<const ULong64_t>+;
+#pragma link C++ class TNDArrayRef<const ULong_t>+;
+#pragma link C++ class TNDArrayRef<const UInt_t>+;
+#pragma link C++ class TNDArrayRef<const UShort_t>+;
+#pragma link C++ class THn+;
+#pragma link C++ class THnT<Float_t>+;
+//#pragma link C++ class THnT<Float16_t>+;
+#pragma link C++ class THnT<Double_t>+;
+//#pragma link C++ class THnT<Double32_t>+;
+#pragma link C++ class THnT<Long64_t>+;
+#pragma link C++ class THnT<Long_t>+;
+#pragma link C++ class THnT<Int_t>+;
+#pragma link C++ class THnT<Short_t>+;
+#pragma link C++ class THnT<Char_t>+;
+#pragma link C++ class THnT<ULong64_t>+;
+#pragma link C++ class THnT<ULong_t>+;
+#pragma link C++ class THnT<UInt_t>+;
+#pragma link C++ class THnT<UShort_t>+;
#pragma link C++ class THnSparse+;
#pragma link C++ class THnSparseT<TArrayD>+;
#pragma link C++ class THnSparseT<TArrayF>+;
@@ -221,7 +277,7 @@
#pragma link C++ function R__H(Int_t);
#pragma link C++ function R__H(const char*);
-#pragma link C++ class ROOT::THnSparseBrowsable;
+#pragma link C++ class ROOT::THnBaseBrowsable;
#pragma link C++ class ROOT::Math::WrappedTF1;
#pragma link C++ class ROOT::Math::WrappedMultiTF1;
diff --git a/hist/hist/inc/THn.h b/hist/hist/inc/THn.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4fd6556e6a8f4fc5946c37cfb60048bfb632c20
--- /dev/null
+++ b/hist/hist/inc/THn.h
@@ -0,0 +1,220 @@
+// @(#)root/hist:$Id$
+// Author: Axel Naumann, Nov 2011
+
+/*************************************************************************
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
+ * All rights reserved. *
+ * *
+ * For the licensing terms see $ROOTSYS/LICENSE. *
+ * For the list of contributors see $ROOTSYS/README/CREDITS. *
+ *************************************************************************/
+
+#ifndef ROOT_THN
+#define ROOT_THN
+
+#ifndef ROOT_THnBase
+#include "THnBase.h"
+#endif
+
+#ifndef ROOT_TNDArray
+#include "TNDArray.h"
+#endif
+
+#ifndef ROOT_TArrayD
+#include "TArrayD.h"
+#endif
+
+#ifndef ROOT_TObjArray
+#include "TObjArray.h"
+#endif
+
+#ifndef ROOT_TAxis
+#include "TAxis.h"
+#endif
+
+#ifndef ROOT_TMath
+#include "TMath.h"
+#endif
+
+class TH1;
+class TH1D;
+class TH2D;
+class TH3D;
+class THnSparse;
+class TF1;
+
+class THn: public THnBase {
+private:
+ THn(const THn&); // Not implemented
+ THn& operator=(const THn&); // Not implemented
+
+protected:
+ void FillBin(Long64_t bin, Double_t w) {
+ // Increment the bin content of "bin" by "w",
+ // return the bin index.
+ GetArray().AddAt(bin, w);
+ if (GetCalculateErrors()) {
+ fSumw2.AddAt(bin, w * w);
+ }
+ FillBinBase(w);
+ }
+
+ THnBase* CloneEmpty(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const;
+
+public:
+ THn(): fCoordBuf() {}
+ THn(const char* name, const char* title, Int_t dim, const Int_t* nbins,
+ const Double_t* xmin, const Double_t* xmax);
+ virtual ~THn();
+
+ static THn* CreateHn(const char* name, const char* title, const TH1* h1) {
+ return (THn*) CreateHnAny(name, title, h1, kFALSE /*THn*/, -1);
+ }
+
+ ROOT::THnBaseBinIter* CreateIter(Bool_t respectAxisRange) const;
+ Long64_t GetNbins() const { return GetArray().GetNbins(); }
+
+ Long64_t GetBin(const Int_t* idx) const {
+ return GetArray().GetBin(idx);
+ }
+ Long64_t GetBin(const Double_t* x) const {
+ for (Int_t d = 0; d < fNdimensions; ++d) {
+ fCoordBuf[d] = GetAxis(d)->FindFixBin(x[d]);
+ }
+ return GetArray().GetBin(fCoordBuf);
+ }
+ Long64_t GetBin(const char* name[]) const {
+ for (Int_t d = 0; d < fNdimensions; ++d) {
+ fCoordBuf[d] = GetAxis(d)->FindBin(name[d]);
+ }
+ return GetArray().GetBin(fCoordBuf);
+ }
+
+ Long64_t GetBin(const Int_t* idx, Bool_t /*allocate*/ = kTRUE) {
+ return const_cast<const THn*>(this)->GetBin(idx);
+ }
+ Long64_t GetBin(const Double_t* x, Bool_t /*allocate*/ = kTRUE) {
+ return const_cast<const THn*>(this)->GetBin(x);
+ }
+ Long64_t GetBin(const char* name[], Bool_t /*allocate*/ = kTRUE) {
+ return const_cast<const THn*>(this)->GetBin(name);
+ }
+
+ void SetBinContent(const Int_t* idx, Double_t v) {
+ // Forwards to THnBase::SetBinContent().
+ // Non-virtual, CINT-compatible replacement of a using declaration.
+ THnBase::SetBinContent(idx, v);
+ }
+ void SetBinContent(Long64_t bin, Double_t v) {
+ GetArray().SetAsDouble(bin, v);
+ }
+ void SetBinError2(Long64_t bin, Double_t e2) {
+ if (!GetCalculateErrors()) Sumw2();
+ fSumw2.At(bin) = e2;
+ }
+ void AddBinContent(Long64_t bin, Double_t v = 1.) {
+ GetArray().AddAt(bin, v);
+ }
+ void AddBinError2(Long64_t bin, Double_t e2) {
+ fSumw2.At(bin) += e2;
+ }
+ Double_t GetBinContent(const Int_t *idx) const {
+ // Forwards to THnBase::GetBinContent() overload.
+ // Non-virtual, CINT-compatible replacement of a using declaration.
+ return THnBase::GetBinContent(idx);
+ }
+ Double_t GetBinContent(Long64_t bin, Int_t* idx = 0) const {
+ // Get the content of bin, and set its index if idx is != 0.
+ if (idx) {
+ const TNDArray& arr = GetArray();
+ Long64_t prevCellSize = arr.GetNbins();
+ for (Int_t i = 0; i < GetNdimensions(); ++i) {
+ Long64_t cellSize = arr.GetCellSize(i);
+ idx[i] = (bin % prevCellSize) / cellSize;
+ prevCellSize = cellSize;
+ }
+ }
+ return GetArray().AtAsDouble(bin);
+ }
+ Double_t GetBinError2(Long64_t linidx) const {
+ return GetCalculateErrors() ? fSumw2.At(linidx) : GetBinContent(linidx);
+ }
+
+ virtual const TNDArray& GetArray() const = 0;
+ virtual TNDArray& GetArray() = 0;
+
+ void Sumw2();
+
+ TH1D* Projection(Int_t xDim, Option_t* option = "") const {
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(xDim, option);
+ }
+
+ TH2D* Projection(Int_t yDim, Int_t xDim,
+ Option_t* option = "") const {
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(yDim, xDim, option);
+ }
+
+ TH3D* Projection(Int_t xDim, Int_t yDim, Int_t zDim,
+ Option_t* option = "") const {
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(xDim, yDim, zDim, option);
+ }
+
+ THn* Projection(Int_t ndim, const Int_t* dim,
+ Option_t* option = "") const {
+ return (THn*) ProjectionND(ndim, dim, option);
+ }
+
+ THn* Rebin(Int_t group) const {
+ return (THn*) RebinBase(group);
+ }
+ THn* Rebin(const Int_t* group) const {
+ return (THn*) RebinBase(group);
+ }
+
+ void Reset(Option_t* option = "");
+
+protected:
+ TNDArrayT<Double_t> fSumw2; // bin error, lazy allocation happens in TNDArrayT
+ mutable Int_t* fCoordBuf; //! Temporary buffer
+
+ ClassDef(THn, 1); //Base class for multi-dimensional histogram
+};
+
+template <typename T>
+class THnT: public THn {
+public:
+ THnT() {}
+
+ THnT(const char* name, const char* title,
+ Int_t dim, const Int_t* nbins,
+ const Double_t* xmin, const Double_t* xmax):
+ THn(name, title, dim, nbins, xmin, xmax),
+ fArray(dim, nbins, true) {}
+
+ const TNDArray& GetArray() const { return fArray; }
+ TNDArray& GetArray() { return fArray; }
+
+protected:
+ TNDArrayT<T> fArray; // bin content
+ ClassDef(THnT, 1); // multi-dimensional histogram with templated storage
+};
+
+typedef THnT<Float_t> THnF;
+typedef THnT<Double_t> THnD;
+typedef THnT<Char_t> THnC;
+typedef THnT<Short_t> THnS;
+typedef THnT<Int_t> THnI;
+typedef THnT<Long_t> THnL;
+typedef THnT<Long64_t> THnL64;
+
+#endif // ROOT_THN
diff --git a/hist/hist/inc/THnBase.h b/hist/hist/inc/THnBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..576a03f76055fa685cacfc0f01a7da863184387b
--- /dev/null
+++ b/hist/hist/inc/THnBase.h
@@ -0,0 +1,324 @@
+// @(#)root/hist:$Id$
+// Author: Axel Naumann (2011-12-20)
+
+/*************************************************************************
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
+ * All rights reserved. *
+ * *
+ * For the licensing terms see $ROOTSYS/LICENSE. *
+ * For the list of contributors see $ROOTSYS/README/CREDITS. *
+ *************************************************************************/
+
+#ifndef ROOT_THnBase
+#define ROOT_THnBase
+
+/*************************************************************************
+
+ THnBase: Common base class for n-dimensional histogramming.
+ Defines interfaces and algorithms.
+
+*************************************************************************/
+
+
+#ifndef ROOT_TNamed
+#include "TNamed.h"
+#endif
+#ifndef ROOT_TMath
+#include "TMath.h"
+#endif
+#ifndef ROOT_TFitResultPtr
+#include "TFitResultPtr.h"
+#endif
+#ifndef ROOT_TObjArray
+#include "TObjArray.h"
+#endif
+#ifndef ROOT_TArrayD
+#include "TArrayD.h"
+#endif
+
+class TAxis;
+class TH1;
+class TH1D;
+class TH2D;
+class TH3D;
+class TF1;
+class THnIter;
+
+namespace ROOT {
+ class THnBaseBinIter;
+}
+
+class THnBase: public TNamed {
+protected:
+ Int_t fNdimensions; // number of dimensions
+ TObjArray fAxes; // axes of the histogram
+ TObjArray fBrowsables; //! browser-helpers for each axis
+ Double_t fEntries; // number of entries, spread over chunks
+ Double_t fTsumw; // total sum of weights
+ Double_t fTsumw2; // total sum of weights squared; -1 if no errors are calculated
+ TArrayD fTsumwx; // total sum of weight*X for each dimension
+ TArrayD fTsumwx2; // total sum of weight*X*X for each dimension
+ Double_t *fIntegral; //! array with bin weight sums
+ enum {
+ kNoInt,
+ kValidInt,
+ kInvalidInt
+ } fIntegralStatus; //! status of integral
+
+private:
+ THnBase(const THnBase&); // Not implemented
+ THnBase& operator=(const THnBase&); // Not implemented
+
+ protected:
+ THnBase():
+ fNdimensions(0), fEntries(0),
+ fTsumw(0), fTsumw2(-1.), fIntegral(0), fIntegralStatus(kNoInt)
+ {}
+
+ THnBase(const char* name, const char* title, Int_t dim,
+ const Int_t* nbins, const Double_t* xmin, const Double_t* xmax);
+
+ void UpdateXStat(const Double_t *x, Double_t w = 1.) {
+ if (GetCalculateErrors()) {
+ for (Int_t d = 0; d < fNdimensions; ++d) {
+ const Double_t xd = x[d];
+ fTsumwx[d] += w * xd;
+ fTsumwx2[d] += w * xd * xd;
+ }
+ }
+ }
+
+ virtual void FillBin(Long64_t bin, Double_t w) = 0;
+ void FillBinBase(Double_t w) {
+ // Increment the statistics due to filled weight "w",
+ fEntries += 1;
+ if (GetCalculateErrors()) {
+ fTsumw += w;
+ fTsumw2 += w*w;
+ }
+ fIntegralStatus = kInvalidInt;
+ }
+
+ virtual THnBase* CloneEmpty(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const = 0;
+ virtual void Reserve(Long64_t /*nbins*/) {}
+ virtual void SetFilledBins(Long64_t /*nbins*/) {};
+
+ Bool_t CheckConsistency(const THnBase *h, const char *tag) const;
+ TH1* CreateHist(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const;
+ TObject* ProjectionAny(Int_t ndim, const Int_t* dim,
+ Bool_t wantNDim, Option_t* option = "") const;
+ Bool_t PrintBin(Long64_t idx, Int_t* coord, Option_t* options) const;
+ void AddInternal(const THnBase* h, Double_t c, Bool_t rebinned);
+ THnBase* RebinBase(Int_t group) const;
+ THnBase* RebinBase(const Int_t* group) const;
+ void ResetBase(Option_t *option= "");
+
+ static THnBase* CreateHnAny(const char* name, const char* title,
+ const TH1* h1, Bool_t sparse,
+ Int_t ChunkSize = 1024 * 16);
+
+ public:
+ virtual ~THnBase();
+
+ TObjArray* GetListOfAxes() { return &fAxes; }
+ const TObjArray* GetListOfAxes() const { return &fAxes; }
+ TAxis* GetAxis(Int_t dim) const { return (TAxis*)fAxes[dim]; }
+
+ TFitResultPtr Fit(TF1 *f1 ,Option_t *option = "", Option_t *goption = "");
+ TList* GetListOfFunctions() { return 0; }
+
+ virtual ROOT::THnBaseBinIter* CreateIter(Bool_t respectAxisRange) const = 0;
+
+ virtual Long64_t GetNbins() const = 0;
+ Double_t GetEntries() const { return fEntries; }
+ Double_t GetWeightSum() const { return fTsumw; }
+ Int_t GetNdimensions() const { return fNdimensions; }
+ Bool_t GetCalculateErrors() const { return fTsumw2 >= 0.; }
+ void CalculateErrors(Bool_t calc = kTRUE) {
+ // Calculate errors (or not if "calc" == kFALSE)
+ if (calc) Sumw2();
+ else fTsumw2 = -1.;
+ }
+
+ Long64_t Fill(const Double_t *x, Double_t w = 1.) {
+ UpdateXStat(x, w);
+ Long64_t bin = GetBin(x, kTRUE /*alloc*/);
+ FillBin(bin, w);
+ return bin;
+ }
+ Long64_t Fill(const char* name[], Double_t w = 1.) {
+ Long64_t bin = GetBin(name, kTRUE /*alloc*/);
+ FillBin(bin, w);
+ return bin;
+ }
+ void SetBinEdges(Int_t idim, const Double_t* bins);
+ Bool_t IsInRange(Int_t *coord) const;
+ Double_t GetBinError(const Int_t *idx) const { return GetBinError(GetBin(idx)); }
+ Double_t GetBinError(Long64_t linidx) const { return TMath::Sqrt(GetBinError2(linidx)); }
+ void SetBinError(const Int_t* idx, Double_t e) { SetBinError(GetBin(idx), e); }
+ void SetBinError(Long64_t bin, Double_t e) { SetBinError2(bin, e*e); }
+ void AddBinContent(const Int_t* x, Double_t v = 1.) { AddBinContent(GetBin(x), v); }
+ void SetEntries(Double_t entries) { fEntries = entries; }
+ void SetTitle(const char *title);
+
+ Double_t GetBinContent(const Int_t *idx) const { return GetBinContent(GetBin(idx)); } // intentionally non-virtual
+ virtual Double_t GetBinContent(Long64_t bin, Int_t* idx = 0) const = 0;
+ virtual Double_t GetBinError2(Long64_t linidx) const = 0;
+ virtual Long64_t GetBin(const Int_t* idx) const = 0;
+ virtual Long64_t GetBin(const Double_t* x) const = 0;
+ virtual Long64_t GetBin(const char* name[]) const = 0;
+ virtual Long64_t GetBin(const Int_t* idx, Bool_t /*allocate*/ = kTRUE) = 0;
+ virtual Long64_t GetBin(const Double_t* x, Bool_t /*allocate*/ = kTRUE) = 0;
+ virtual Long64_t GetBin(const char* name[], Bool_t /*allocate*/ = kTRUE) = 0;
+
+ void SetBinContent(const Int_t* idx, Double_t v) { SetBinContent(GetBin(idx), v); } // intentionally non-virtual
+ virtual void SetBinContent(Long64_t bin, Double_t v) = 0;
+ virtual void SetBinError2(Long64_t bin, Double_t e2) = 0;
+ virtual void AddBinError2(Long64_t bin, Double_t e2) = 0;
+ virtual void AddBinContent(Long64_t bin, Double_t v = 1.) = 0;
+
+ Double_t GetSumw() const { return fTsumw; }
+ Double_t GetSumw2() const { return fTsumw2; }
+ Double_t GetSumwx(Int_t dim) const { return fTsumwx[dim]; }
+ Double_t GetSumwx2(Int_t dim) const { return fTsumwx2[dim]; }
+
+ TH1D* Projection(Int_t xDim, Option_t* option = "") const {
+ // Project all bins into a 1-dimensional histogram,
+ // keeping only axis "xDim".
+ // If "option" contains "E" errors will be calculated.
+ // "A" ranges of the taget axes will be ignored.
+ // "O" original axis range of the taget axes will be
+ // kept, but only bins inside the selected range
+ // will be filled.
+ return (TH1D*) ProjectionAny(1, &xDim, false, option);
+ }
+
+ TH2D* Projection(Int_t yDim, Int_t xDim,
+ Option_t* option = "") const {
+ // Project all bins into a 2-dimensional histogram,
+ // keeping only axes "xDim" and "yDim".
+ //
+ // WARNING: just like TH3::Project3D("yx") and TTree::Draw("y:x"),
+ // Projection(y,x) uses the first argument to define the y-axis and the
+ // second for the x-axis!
+ //
+ // If "option" contains "E" errors will be calculated.
+ // "A" ranges of the taget axes will be ignored.
+
+ const Int_t dim[2] = {xDim, yDim};
+ return (TH2D*) ProjectionAny(2, dim, false, option);
+ }
+
+ TH3D* Projection(Int_t xDim, Int_t yDim, Int_t zDim,
+ Option_t* option = "") const {
+ // Project all bins into a 3-dimensional histogram,
+ // keeping only axes "xDim", "yDim", and "zDim".
+ // If "option" contains "E" errors will be calculated.
+ // "A" ranges of the taget axes will be ignored.
+ // "O" original axis range of the taget axes will be
+ // kept, but only bins inside the selected range
+ // will be filled.
+
+ const Int_t dim[3] = {xDim, yDim, zDim};
+ return (TH3D*) ProjectionAny(3, dim, false, option);
+ }
+
+ THnBase* ProjectionND(Int_t ndim, const Int_t* dim,
+ Option_t* option = "") const {
+ return (THnBase*)ProjectionAny(ndim, dim, kTRUE /*wantNDim*/, option);
+ }
+
+ Long64_t Merge(TCollection* list);
+
+ void Scale(Double_t c);
+ void Add(const THnBase* h, Double_t c=1.);
+ void Add(const TH1* hist, Double_t c=1.);
+ void Multiply(const THnBase* h);
+ void Multiply(TF1* f, Double_t c = 1.);
+ void Divide(const THnBase* h);
+ void Divide(const THnBase* h1, const THnBase* h2, Double_t c1 = 1., Double_t c2 = 1., Option_t* option="");
+ void RebinnedAdd(const THnBase* h, Double_t c=1.);
+
+ virtual void Reset(Option_t* option = "") = 0;
+ virtual void Sumw2() = 0;
+
+ Double_t ComputeIntegral();
+ void GetRandom(Double_t *rand, Bool_t subBinRandom = kTRUE);
+
+ void Print(Option_t* option = "") const;
+ void PrintEntries(Long64_t from = 0, Long64_t howmany = -1, Option_t* options = 0) const;
+ void PrintBin(Int_t* coord, Option_t* options) const {
+ PrintBin(-1, coord, options);
+ }
+ void PrintBin(Long64_t idx, Option_t* options) const;
+
+ void Browse(TBrowser *b);
+ Bool_t IsFolder() const { return kTRUE; }
+
+ //void Draw(Option_t* option = "");
+
+ ClassDef(THnBase, 1); // Common base for n-dimensional histogram
+
+ friend class THnIter;
+};
+
+namespace ROOT {
+ // Helper class for browing THnBase objects
+ class THnBaseBrowsable: public TNamed {
+ public:
+ THnBaseBrowsable(THnBase* hist, Int_t axis);
+ ~THnBaseBrowsable();
+ void Browse(TBrowser *b);
+ Bool_t IsFolder() const { return kFALSE; }
+
+ private:
+ THnBase* fHist; // Original histogram
+ Int_t fAxis; // Axis to visualize
+ TH1* fProj; // Projection result
+ ClassDef(THnBaseBrowsable, 0); // Browser-helper for THnBase
+ };
+
+ // Base class for iterating over THnBase bins
+ class THnBaseBinIter {
+ public:
+ THnBaseBinIter(Bool_t respectAxisRange):
+ fRespectAxisRange(respectAxisRange), fHaveSkippedBin(kFALSE) {}
+ virtual ~THnBaseBinIter();
+ Bool_t HaveSkippedBin() const { return fHaveSkippedBin; }
+ Bool_t RespectsAxisRange() const { return fRespectAxisRange; }
+
+ virtual Int_t GetCoord(Int_t dim) const = 0;
+ virtual Long64_t Next(Int_t* coord = 0) = 0;
+
+ protected:
+ Bool_t fRespectAxisRange;
+ Bool_t fHaveSkippedBin;
+ };
+}
+
+class THnIter: public TObject {
+public:
+ THnIter(const THnBase* hist, Bool_t respectAxisRange = kFALSE):
+ fIter(hist->CreateIter(respectAxisRange)) {}
+ virtual ~THnIter();
+
+ Long64_t Next(Int_t* coord = 0) {
+ // Return the next bin's index.
+ // If provided, set coord to that bin's coordinates (bin indexes).
+ // I.e. coord must point to Int_t[hist->GetNdimensions()]
+ // Returns -1 when all bins have been visited.
+ return fIter->Next(coord);
+ }
+
+ Int_t GetCoord(Int_t dim) const { return fIter->GetCoord(dim); }
+ Bool_t HaveSkippedBin() const { return fIter->HaveSkippedBin(); }
+ Bool_t RespectsAxisRange() const { return fIter->RespectsAxisRange(); }
+
+private:
+ ROOT::THnBaseBinIter* fIter;
+ ClassDef(THnIter, 0); //Iterator over bins of a THnBase.
+};
+
+#endif // ROOT_THnBase
diff --git a/hist/hist/inc/THnSparse.h b/hist/hist/inc/THnSparse.h
index 9171107763fcabde607059499cf8a1f816fffe59..5ffa1d378266192bb653c27d7bce09b95d91c21c 100644
--- a/hist/hist/inc/THnSparse.h
+++ b/hist/hist/inc/THnSparse.h
@@ -2,7 +2,7 @@
// Author: Axel Naumann (2007-09-11)
/*************************************************************************
- * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers. *
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
@@ -20,21 +20,12 @@
*************************************************************************/
+#ifndef ROOT_THnBase
+#include "THnBase.h"
+#endif
#ifndef ROOT_TExMap
#include "TExMap.h"
#endif
-#ifndef ROOT_TNamed
-#include "TNamed.h"
-#endif
-#ifndef ROOT_TObjArray
-#include "TObjArray.h"
-#endif
-#ifndef ROOT_TArrayD
-#include "TArrayD.h"
-#endif
-#ifndef ROOT_TFitResultPtr
-#include "TFitResultPtr.h"
-#endif
#ifndef ROOT_THnSparse_Internal
#include "THnSparse_Internal.h"
#endif
@@ -56,40 +47,16 @@
#include "TArrayC.h"
#endif
-class TAxis;
-class TCollection;
-class TH1;
-class TH1D;
-class TH2D;
-class TH3D;
-class TF1;
-
-
-
class THnSparseCompactBinCoord;
-class THnSparse: public TNamed {
+class THnSparse: public THnBase {
private:
- Int_t fNdimensions; // number of dimensions
Int_t fChunkSize; // number of entries for each chunk
Long64_t fFilledBins; // number of filled bins
- TObjArray fAxes; // axes of the histogram
- TObjArray fBrowsables; //! browser-helpers for each axis
TObjArray fBinContent; // array of THnSparseArrayChunk
TExMap fBins; //! filled bins
TExMap fBinsContinued;//! filled bins for non-unique hashes, containing pairs of (bin index 0, bin index 1)
- Double_t fEntries; // number of entries, spread over chunks
- Double_t fTsumw; // total sum of weights
- Double_t fTsumw2; // total sum of weights squared; -1 if no errors are calculated
- TArrayD fTsumwx; // total sum of weight*X for each dimension
- TArrayD fTsumwx2; // total sum of weight*X*X for each dimension
THnSparseCompactBinCoord *fCompactCoord; //! compact coordinate
- Double_t *fIntegral; //! array with bin weight sums
- enum {
- kNoInt,
- kValidInt,
- kInvalidInt
- } fIntegralStatus; //! status of integral
THnSparse(const THnSparse&); // Not implemented
THnSparse& operator=(const THnSparse&); // Not implemented
@@ -100,146 +67,109 @@ class THnSparse: public TNamed {
THnSparse(const char* name, const char* title, Int_t dim,
const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
Int_t chunksize);
- Int_t GetChunkSize() const { return fChunkSize; }
THnSparseCompactBinCoord* GetCompactCoord() const;
THnSparseArrayChunk* GetChunk(Int_t idx) const {
return (THnSparseArrayChunk*) fBinContent[idx]; }
THnSparseArrayChunk* AddChunk();
+ void Reserve(Long64_t nbins);
void FillExMap();
virtual TArray* GenerateArray() const = 0;
Long64_t GetBinIndexForCurrentBin(Bool_t allocate);
- Long64_t Fill(Long64_t bin, Double_t w) {
+ void FillBin(Long64_t bin, Double_t w) {
// Increment the bin content of "bin" by "w",
// return the bin index.
- fEntries += 1;
- if (GetCalculateErrors()) {
- fTsumw += w;
- fTsumw2 += w*w;
- }
- fIntegralStatus = kInvalidInt;
THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
chunk->AddBinContent(bin % fChunkSize, w);
- return bin;
+ FillBinBase(w);
}
- THnSparse* CloneEmpty(const char* name, const char* title,
- const TObjArray* axes, Int_t chunksize,
- Bool_t keepTargetAxis) const;
-
- Bool_t CheckConsistency(const THnSparse *h, const char *tag) const;
- Bool_t IsInRange(Int_t *coord) const;
- TH1* CreateHist(const char* name, const char* title,
- const TObjArray* axes, Bool_t keepTargetAxis) const;
- TObject* ProjectionAny(Int_t ndim, const Int_t* dim,
- Bool_t wantSparse, Option_t* option = "") const;
- Bool_t PrintBin(Long64_t idx, Int_t* coord, Option_t* options) const;
- void AddInternal(const THnSparse* h, Double_t c, Bool_t rebinned);
+ THnBase* CloneEmpty(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const;
public:
virtual ~THnSparse();
static THnSparse* CreateSparse(const char* name, const char* title,
- const TH1* h1, Int_t ChunkSize = 1024 * 16);
-
- Int_t GetNChunks() const { return fBinContent.GetEntriesFast(); }
- TObjArray* GetListOfAxes() { return &fAxes; }
- TAxis* GetAxis(Int_t dim) const { return (TAxis*)fAxes[dim]; }
-
- Long64_t Fill(const Double_t *x, Double_t w = 1.) {
- if (GetCalculateErrors()) {
- for (Int_t d = 0; d < fNdimensions; ++d) {
- const Double_t xd = x[d];
- fTsumwx[d] += w * xd;
- fTsumwx2[d] += w * xd * xd;
- }
- }
- return Fill(GetBin(x), w);
- }
- Long64_t Fill(const char* name[], Double_t w = 1.) {
- return Fill(GetBin(name), w);
+ const TH1* h1, Int_t chunkSize = 1024 * 16) {
+ return (THnSparse*) CreateHnAny(name, title, h1, kTRUE /*sparse*/,
+ chunkSize);
}
- TFitResultPtr Fit(TF1 *f1 ,Option_t *option = "", Option_t *goption = "");
- TList* GetListOfFunctions() { return 0; }
+ Int_t GetChunkSize() const { return fChunkSize; }
+ Int_t GetNChunks() const { return fBinContent.GetEntriesFast(); }
+
+ ROOT::THnBaseBinIter* CreateIter(Bool_t respectAxisRange) const;
- Double_t GetEntries() const { return fEntries; }
- Double_t GetWeightSum() const { return fTsumw; }
Long64_t GetNbins() const { return fFilledBins; }
- Int_t GetNdimensions() const { return fNdimensions; }
- Bool_t GetCalculateErrors() const { return fTsumw2 >= 0.; }
- void CalculateErrors(Bool_t calc = kTRUE) {
- // Calculate errors (or not if "calc" == kFALSE)
- if (calc) Sumw2();
- else fTsumw2 = -1.;
- }
+ void SetFilledBins(Long64_t nbins) { fFilledBins = nbins; }
+ Long64_t GetBin(const Int_t* idx) const { return const_cast<THnSparse*>(this)->GetBin(idx, kFALSE); }
+ Long64_t GetBin(const Double_t* x) const { return const_cast<THnSparse*>(this)->GetBin(x, kFALSE); }
+ Long64_t GetBin(const char* name[]) const { return const_cast<THnSparse*>(this)->GetBin(name, kFALSE); }
Long64_t GetBin(const Int_t* idx, Bool_t allocate = kTRUE);
Long64_t GetBin(const Double_t* x, Bool_t allocate = kTRUE);
Long64_t GetBin(const char* name[], Bool_t allocate = kTRUE);
- void SetBinEdges(Int_t idim, const Double_t* bins);
- void SetBinContent(const Int_t* x, Double_t v);
+ void SetBinContent(const Int_t* idx, Double_t v) {
+ // Forwards to THnBase::SetBinContent().
+ // Non-virtual, CINT-compatible replacement of a using declaration.
+ THnBase::SetBinContent(idx, v);
+ }
void SetBinContent(Long64_t bin, Double_t v);
- void SetBinError(const Int_t* x, Double_t e);
- void SetBinError(Long64_t bin, Double_t e);
- void AddBinContent(const Int_t* x, Double_t v = 1.);
+ void SetBinError2(Long64_t bin, Double_t e2);
void AddBinContent(Long64_t bin, Double_t v = 1.);
- void SetEntries(Double_t entries) { fEntries = entries; }
- void SetTitle(const char *title);
+ void AddBinError2(Long64_t bin, Double_t e2);
- Double_t GetBinContent(const Int_t *idx) const;
+ Double_t GetBinContent(const Int_t *idx) const {
+ // Forwards to THnBase::GetBinContent() overload.
+ // Non-virtual, CINT-compatible replacement of a using declaration.
+ return THnBase::GetBinContent(idx);
+ }
Double_t GetBinContent(Long64_t bin, Int_t* idx = 0) const;
- Double_t GetBinError(const Int_t *idx) const;
- Double_t GetBinError(Long64_t linidx) const;
+ Double_t GetBinError2(Long64_t linidx) const;
Double_t GetSparseFractionBins() const;
Double_t GetSparseFractionMem() const;
- Double_t GetSumw() const { return fTsumw; }
- Double_t GetSumw2() const { return fTsumw2; }
- Double_t GetSumwx(Int_t dim) const { return fTsumwx[dim]; }
- Double_t GetSumwx2(Int_t dim) const { return fTsumwx2[dim]; }
+ TH1D* Projection(Int_t xDim, Option_t* option = "") const{
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(xDim, option);
+ }
- TH1D* Projection(Int_t xDim, Option_t* option = "") const;
TH2D* Projection(Int_t yDim, Int_t xDim,
- Option_t* option = "") const;
- TH3D* Projection(Int_t xDim, Int_t yDim, Int_t zDim,
- Option_t* option = "") const;
- THnSparse* Projection(Int_t ndim, const Int_t* dim,
- Option_t* option = "") const;
+ Option_t* option = "") const {
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(yDim, xDim, option);
+ }
- THnSparse* Rebin(Int_t group) const;
- THnSparse* Rebin(const Int_t* group) const;
+ TH3D* Projection(Int_t xDim, Int_t yDim, Int_t zDim,
+ Option_t* option = "") const {
+ // Forwards to THnBase::Projection().
+ // Non-virtual, as a CINT-compatible replacement of a using
+ // declaration.
+ return THnBase::Projection(xDim, yDim, zDim, option);
+ }
- Long64_t Merge(TCollection* list);
+ THnSparse* Projection(Int_t ndim, const Int_t* dim,
+ Option_t* option = "") const {
+ return (THnSparse*) ProjectionND(ndim, dim, option);
+ }
- void Scale(Double_t c);
- void Add(const THnSparse* h, Double_t c=1.);
- void Multiply(const THnSparse* h);
- void Multiply(TF1* f, Double_t c = 1.);
- void Divide(const THnSparse* h);
- void Divide(const THnSparse* h1, const THnSparse* h2, Double_t c1 = 1., Double_t c2 = 1., Option_t* option="");
- void RebinnedAdd(const THnSparse* h, Double_t c=1.);
+ THnSparse* Rebin(Int_t group) const {
+ return (THnSparse*) RebinBase(group);
+ }
+ THnSparse* Rebin(const Int_t* group) const {
+ return (THnSparse*) RebinBase(group);
+ }
void Reset(Option_t* option = "");
void Sumw2();
- Double_t ComputeIntegral();
- void GetRandom(Double_t *rand, Bool_t subBinRandom = kTRUE);
-
- void Print(Option_t* option = "") const;
- void PrintEntries(Long64_t from = 0, Long64_t howmany = -1, Option_t* options = 0) const;
- void PrintBin(Int_t* coord, Option_t* options) const {
- PrintBin(-1, coord, options);
- }
- void PrintBin(Long64_t idx, Option_t* options) const;
-
- void Browse(TBrowser *b);
- Bool_t IsFolder() const { return kTRUE; }
-
- //void Draw(Option_t* option = "");
-
- ClassDef(THnSparse, 2); // Interfaces of sparse n-dimensional histogram
+ ClassDef(THnSparse, 3); // Interfaces of sparse n-dimensional histogram
};
@@ -257,7 +187,7 @@ class THnSparse: public TNamed {
// Typedefs exist for template parematers with ROOT's generic types:
//
// Templated name Typedef Bin content type
-// THnSparseT<TArrayC> THnSparseC Char_r
+// THnSparseT<TArrayC> THnSparseC Char_t
// THnSparseT<TArrayS> THnSparseS Short_t
// THnSparseT<TArrayI> THnSparseI Int_t
// THnSparseT<TArrayL> THnSparseL Long_t
@@ -295,4 +225,5 @@ typedef THnSparseT<TArrayI> THnSparseI;
typedef THnSparseT<TArrayS> THnSparseS;
typedef THnSparseT<TArrayC> THnSparseC;
+
#endif // ROOT_THnSparse
diff --git a/hist/hist/inc/THnSparse_Internal.h b/hist/hist/inc/THnSparse_Internal.h
index fd2e47040f2a0b8465ea143ff9e9b5e24e39cb9f..b3e6922ca91d4c0ad115059feb2a497d806fbead 100644
--- a/hist/hist/inc/THnSparse_Internal.h
+++ b/hist/hist/inc/THnSparse_Internal.h
@@ -2,7 +2,7 @@
// Author: Axel Naumann (2007-09-11)
/*************************************************************************
- * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers. *
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
@@ -65,21 +65,5 @@ class THnSparseArrayChunk: public TObject {
ClassDef(THnSparseArrayChunk, 1); // chunks of linearized bins
};
-
-namespace ROOT {
- class THnSparseBrowsable: public TNamed {
- public:
- THnSparseBrowsable(THnSparse* hist, Int_t axis);
- ~THnSparseBrowsable();
- void Browse(TBrowser *b);
- Bool_t IsFolder() const { return kFALSE; }
-
- private:
- THnSparse* fHist; // Original histogram
- Int_t fAxis; // Axis to visualize
- TH1* fProj; // Projection result
- ClassDef(THnSparseBrowsable, 0); // Browser-helper for THnSparse
- };
-}
#endif // ROOT_THnSparse_Internal
diff --git a/hist/hist/inc/TNDArray.h b/hist/hist/inc/TNDArray.h
new file mode 100644
index 0000000000000000000000000000000000000000..6cfaaf0bc0b851aa83e2b3fd81c846a2a85a4096
--- /dev/null
+++ b/hist/hist/inc/TNDArray.h
@@ -0,0 +1,200 @@
+// @(#)root/hist:$Id$
+// Author: Axel Naumann, Nov 2011
+
+/*************************************************************************
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
+ * All rights reserved. *
+ * *
+ * For the licensing terms see $ROOTSYS/LICENSE. *
+ * For the list of contributors see $ROOTSYS/README/CREDITS. *
+ *************************************************************************/
+
+#ifndef ROOT_TNDArray
+#define ROOT_TNDArray
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+#ifndef ROOT_TError
+#include "TError.h"
+#endif
+
+//////////////////////////////////////////////////////////////////////////
+// //
+// TNDArray //
+// //
+// N-Dim array class. //
+// //
+// Storage layout: //
+// Assume 3 dimensions, array sizes 2, 4 and 3 i.e. 24 bins: //
+// Data is stored as [0,0,0], [0,0,1], [0,0,2], [0,1,0],... //
+// //
+// fSizes stores the combined size of each bin in a dimension, i.e. in //
+// above example it would contain 24, 12, 3, 1. //
+// //
+// Storage is allocated lazily, only when data is written to the array. //
+// //
+// TNDArrayRef gives access to a sub-dimension, e.g. arr[0][1] in above //
+// three-dimensional example, up to an element with conversion operator //
+// to double: double value = arr[0][1][2]; //
+// //
+//////////////////////////////////////////////////////////////////////////
+
+// Array layout:
+// nbins[0] = 2, nbins[1] = 4, nbins[2] = 3 => 24 bins
+//
+// fSizes: 24, 12, 3 [, 1
+
+class TNDArray: public TObject {
+public:
+ TNDArray(): fNdimPlusOne(), fSizes() {}
+
+ TNDArray(Int_t ndim, const Int_t* nbins, bool addOverflow = false):
+ fNdimPlusOne(), fSizes() {
+ TNDArray::Init(ndim, nbins, addOverflow);
+ }
+ ~TNDArray() {
+ delete[] fSizes;
+ }
+
+ virtual void Init(Int_t ndim, const Int_t* nbins, bool addOverflow = false) {
+ // Calculate fSize based on ndim dimensions, nbins for each dimension,
+ // possibly adding over- and underflow bin to each dimensions' nbins.
+ delete[] fSizes;
+ fNdimPlusOne = ndim + 1;
+ fSizes = new Long64_t[ndim + 1];
+ Int_t overBins = addOverflow ? 2 : 0;
+ fSizes[ndim] = 1;
+ for (Int_t i = 0; i < ndim; ++i) {
+ fSizes[ndim - i - 1] = fSizes[ndim - i] * (nbins[ndim - i - 1] + overBins);
+ }
+ }
+
+ virtual void Reset(Option_t* option = "") = 0;
+
+ Int_t GetNdimensions() const { return fNdimPlusOne - 1; }
+ Long64_t GetNbins() const { return fSizes[0]; }
+ Long64_t GetCellSize(Int_t dim) const { return fSizes[dim + 1]; }
+
+ Long64_t GetBin(const Int_t* idx) const {
+ // Get the linear bin number for each dimension's bin index
+ Long64_t bin = idx[fNdimPlusOne - 2];
+ for (Int_t d = 0; d < fNdimPlusOne - 2; ++d) {
+ bin += fSizes[d + 1] * idx[d];
+ }
+ return bin;
+ }
+
+ virtual Double_t AtAsDouble(ULong64_t linidx) const = 0;
+ virtual void SetAsDouble(ULong64_t linidx, Double_t value) = 0;
+ virtual void AddAt(ULong64_t linidx, Double_t value) = 0;
+
+protected:
+ Int_t fNdimPlusOne; // Number of dimensions plus one
+ Long64_t* fSizes; //[fNdimPlusOne] bin count
+ ClassDef(TNDArray, 1); //Base for n-dimensional array
+};
+
+template <typename T>
+class TNDArrayRef {
+public:
+ TNDArrayRef(T* data, const Long64_t* sizes):
+ fData(data), fSizes(sizes) {}
+
+ TNDArrayRef<T> operator[] (Int_t idx) const {
+ if (!fData) return TNDArrayRef<T>(0, 0);
+ R__ASSERT(idx < fSizes[-1] / fSizes[0] && "index out of range!");
+ return TNDArrayRef<T>(fData + idx * fSizes[0], (fSizes[0] == 1) ? 0 : (fSizes + 1));
+ }
+ operator T() const {
+ if (!fData) return T();
+ R__ASSERT(fSizes == 0 && "Element operator can only be used on non-array element. Missing an operator[] level?");
+ return *fData;
+ }
+
+private:
+ T* fData; // pointer into TNDArray's fData
+ const Long64_t* fSizes; // pointer into TNDArray's fSizes
+ ClassDefNV(TNDArrayRef, 0); // subdimension of a TNDArray
+};
+
+template <typename T>
+class TNDArrayT: public TNDArray {
+public:
+ TNDArrayT(): fNumData(), fData() {}
+
+ TNDArrayT(Int_t ndim, const Int_t* nbins, bool addOverflow = false):
+ TNDArray(ndim, nbins, addOverflow),
+ fNumData(), fData() {
+ fNumData = fSizes[0];
+ }
+ ~TNDArrayT() {
+ delete[] fData;
+ }
+
+ void Init(Int_t ndim, const Int_t* nbins, bool addOverflow = false) {
+ delete[] fData;
+ fData = 0;
+ TNDArray::Init(ndim, nbins, addOverflow);
+ fNumData = fSizes[0];
+ }
+
+ void Reset(Option_t* /*option*/ = "") {
+ // Reset the content
+
+ // Use placement-new with value initialization:
+ if (fData) {
+ new (fData) T[fNumData]();
+ }
+ }
+
+#ifndef __CINT__
+ TNDArrayRef<const T> operator[] (Int_t idx) const {
+ if (!fData) return TNDArrayRef<const T>(0, 0);
+ R__ASSERT(idx < fSizes[0] / fSizes[1] && "index out of range!");
+ return TNDArrayRef<const T>(fData + idx * fSizes[1], fSizes + 2);
+ }
+
+ TNDArrayRef<T> operator[] (Int_t idx) {
+ if (!fData) return TNDArrayRef<T>(0, 0);
+ R__ASSERT(idx < fSizes[0] / fSizes[1] && "index out of range!");
+ return TNDArrayRef<T>(fData + idx * fSizes[1], fSizes + 2);
+ }
+#endif // __CINT__
+
+ T At(const Int_t* idx) const {
+ return At(GetBin(idx));
+ }
+ T& At(const Int_t* idx) {
+ return At(GetBin(idx));
+ }
+ T At(ULong64_t linidx) const {
+ if (!fData) return T();
+ return fData[linidx];
+ }
+ T& At(ULong64_t linidx) {
+ if (!fData) fData = new T[fNumData]();
+ return fData[linidx];
+ }
+
+ Double_t AtAsDouble(ULong64_t linidx) const {
+ if (!fData) return 0.;
+ return fData[linidx];
+ }
+ void SetAsDouble(ULong64_t linidx, Double_t value) {
+ if (!fData) fData = new T[fNumData]();
+ fData[linidx] = (T) value;
+ }
+ void AddAt(ULong64_t linidx, Double_t value) {
+ if (!fData) fData = new T[fNumData]();
+ fData[linidx] += (T) value;
+ }
+
+protected:
+ int fNumData; // number of bins, product of fSizes
+ T* fData; //[fNumData] data
+ ClassDef(TNDArray, 1); // N-dimensional array
+};
+
+
+#endif // ROOT_TNDArray
diff --git a/hist/hist/src/HFitImpl.cxx b/hist/hist/src/HFitImpl.cxx
index 1042083257ee3692a3314552a9bf337837d3dad2..e20273acf5c568cf84581e97abf230a5e8ec2897 100644
--- a/hist/hist/src/HFitImpl.cxx
+++ b/hist/hist/src/HFitImpl.cxx
@@ -11,7 +11,7 @@
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TGraph2D.h"
-#include "THnSparse.h"
+#include "THnBase.h"
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
@@ -48,7 +48,7 @@ namespace HFit {
int GetDimension(const TGraph * ) { return 1; }
int GetDimension(const TMultiGraph * ) { return 1; }
int GetDimension(const TGraph2D * ) { return 2; }
- int GetDimension(const THnSparse * s1) { return s1->GetNdimensions(); }
+ int GetDimension(const THnBase * s1) { return s1->GetNdimensions(); }
int CheckFitFunction(const TF1 * f1, int hdim);
@@ -61,7 +61,7 @@ namespace HFit {
void GetDrawingRange(TGraph * gr, ROOT::Fit::DataRange & range);
void GetDrawingRange(TMultiGraph * mg, ROOT::Fit::DataRange & range);
void GetDrawingRange(TGraph2D * gr, ROOT::Fit::DataRange & range);
- void GetDrawingRange(THnSparse * s, ROOT::Fit::DataRange & range);
+ void GetDrawingRange(THnBase * s, ROOT::Fit::DataRange & range);
template <class FitObject>
@@ -501,7 +501,7 @@ void HFit::GetDrawingRange(TGraph2D * gr, ROOT::Fit::DataRange & range) {
if (h1) HFit::GetDrawingRange(h1, range);
}
-void HFit::GetDrawingRange(THnSparse * s1, ROOT::Fit::DataRange & range) {
+void HFit::GetDrawingRange(THnBase * s1, ROOT::Fit::DataRange & range) {
// get range from histogram and update the DataRange class
// if a ranges already exist in that dimension use that one
@@ -839,7 +839,7 @@ TFitResultPtr ROOT::Fit::FitObject(TGraph2D * gr, TF1 *f1 , Foption_t & foption
return HFit::Fit(gr,f1,foption,moption,goption,range);
}
-TFitResultPtr ROOT::Fit::FitObject(THnSparse * s1, TF1 *f1 , Foption_t & foption , const ROOT::Math::MinimizerOptions & moption, const char *goption, ROOT::Fit::DataRange & range) {
+TFitResultPtr ROOT::Fit::FitObject(THnBase * s1, TF1 *f1 , Foption_t & foption , const ROOT::Math::MinimizerOptions & moption, const char *goption, ROOT::Fit::DataRange & range) {
// sparse histogram fitting
return HFit::Fit(s1,f1,foption,moption,goption,range);
}
diff --git a/hist/hist/src/HFitInterface.cxx b/hist/hist/src/HFitInterface.cxx
index 708aa172960e7ac60c8731e33f7eb12e551be55a..5ec548d0be33ff89dedb2b0d21b04e1dc7cdce81 100644
--- a/hist/hist/src/HFitInterface.cxx
+++ b/hist/hist/src/HFitInterface.cxx
@@ -23,7 +23,7 @@
#include <cmath>
#include "TH1.h"
-#include "THnSparse.h"
+#include "THnBase.h"
#include "TF1.h"
#include "TGraph2D.h"
#include "TGraph.h"
@@ -700,7 +700,7 @@ void FillData(SparseData & dv, const TH1 * h1, TF1 * /*func*/)
}
}
-void FillData(SparseData & dv, const THnSparse * h1, TF1 * /*func*/)
+void FillData(SparseData & dv, const THnBase * h1, TF1 * /*func*/)
{
const int dim = h1->GetNdimensions();
std::vector<double> min(dim);
@@ -743,9 +743,9 @@ void FillData(SparseData & dv, const THnSparse * h1, TF1 * /*func*/)
}
}
-void FillData(BinData & dv, const THnSparse * s1, TF1 * func)
+void FillData(BinData & dv, const THnBase * s1, TF1 * func)
{
- // Fill the Range of the THnSparse
+ // Fill the Range of the THnBase
unsigned int const ndim = s1->GetNdimensions();
std::vector<double> xmin(ndim);
std::vector<double> xmax(ndim);
@@ -767,7 +767,7 @@ void FillData(BinData & dv, const THnSparse * s1, TF1 * func)
ROOT::Fit::SparseData d(ndim, &xmin[0], &xmax[0]);
ROOT::Fit::FillData(d, s1, func);
-// cout << "FillData(BinData & dv, const THnSparse * s1, TF1 * func) (1)" << endl;
+// cout << "FillData(BinData & dv, const THnBase * s1, TF1 * func) (1)" << endl;
// Create the bin data from the sparse data
d.GetBinDataIntegral(dv);
diff --git a/hist/hist/src/THn.cxx b/hist/hist/src/THn.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..332db753f21868f83b679f63f50b1e328f102b1d
--- /dev/null
+++ b/hist/hist/src/THn.cxx
@@ -0,0 +1,193 @@
+// @(#)root/hist:$Id$
+// Author: Axel Naumann (2011-12-13)
+
+/*************************************************************************
+ * Copyright (C) 1995-2012, 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 "THn.h"
+
+#include "TClass.h"
+
+namespace {
+ struct CounterRange_t {
+ Int_t i;
+ Int_t first;
+ Int_t last;
+ Int_t len;
+ Long64_t cellSize;
+ };
+
+ class THnBinIter: public ROOT::THnBaseBinIter {
+ public:
+ THnBinIter(Int_t dim, const TObjArray* axes, const TNDArray* arr,
+ Bool_t respectAxisRange);
+ ~THnBinIter() { delete [] fCounter; }
+
+ Long64_t Next(Int_t* coord = 0);
+ Int_t GetCoord(Int_t dim) const { return fCounter[dim].i; }
+
+ Int_t fNdimensions;
+ Long64_t fIndex;
+ const TNDArray* fArray;
+ CounterRange_t* fCounter;
+ };
+
+
+ //______________________________________________________________________________
+ THnBinIter::THnBinIter(Int_t dim, const TObjArray* axes,
+ const TNDArray* arr, Bool_t respectAxisRange):
+ ROOT::THnBaseBinIter(respectAxisRange),
+ fNdimensions(dim), fIndex(-1), fArray(arr) {
+ fCounter = new CounterRange_t[dim]();
+ for (Int_t i = 0; i < dim; ++i) {
+ TAxis *axis = (TAxis*) axes->At(i);
+ fCounter[i].len = axis->GetNbins() + 2;
+ fCounter[i].cellSize = arr->GetCellSize(i);
+ if (!respectAxisRange || !axis->TestBit(TAxis::kAxisRange)) {
+ fCounter[i].first = 0;
+ fCounter[i].last = fCounter[i].len - 1;
+ fCounter[i].i = 0;
+ continue;
+ }
+ fHaveSkippedBin = kTRUE;
+ Int_t min = axis->GetFirst();
+ Int_t max = axis->GetLast();
+ if (min == 0 && max == 0) {
+ // special case where TAxis::SetBit(kAxisRange) and
+ // over- and underflow bins are de-selected.
+ // first and last are == 0 due to axis12->SetRange(1, axis12->GetNbins());
+ min = 1;
+ max = axis->GetNbins();
+ }
+ fCounter[i].first = min;
+ fCounter[i].last = max;
+ fCounter[i].i = min;
+ fIndex += fCounter[i].first * fCounter[i].cellSize;
+ }
+ // First Next() will increment it:
+ --fCounter[dim - 1].i;
+ }
+
+ //______________________________________________________________________________
+ Long64_t THnBinIter::Next(Int_t* coord /*= 0*/) {
+ if (fNdimensions < 0) return -1; // end
+ ++fCounter[fNdimensions - 1].i;
+ ++fIndex;
+ // Wrap around if needed
+ for (Int_t d = fNdimensions - 1; d > 0 && fCounter[d].i > fCounter[d].last; --d) {
+ // We skip last + 1..size and 0..first - 1, adjust fIndex
+ Int_t skippedCells = fCounter[d].len - (fCounter[d].last + 1);
+ skippedCells += fCounter[d].first;
+ fIndex += skippedCells * fCounter[d].cellSize;
+ fCounter[d].i = fCounter[d].first;
+ ++fCounter[d - 1].i;
+ }
+ if (fCounter[0].i > fCounter[0].last) {
+ fNdimensions = -1;
+ return -1;
+ }
+ if (coord) {
+ for (Int_t d = 0; d < fNdimensions; ++d) {
+ coord[d] = fCounter[d].i;
+ }
+ }
+ return fIndex;
+ }
+} // unnamed namespce
+
+
+
+//______________________________________________________________________________
+THn::THn(const char* name, const char* title,
+ Int_t dim, const Int_t* nbins,
+ const Double_t* xmin, const Double_t* xmax):
+ THnBase(name, title, dim, nbins, xmin, xmax),
+ fSumw2(dim, nbins, kTRUE /*overflow*/) {
+ // Construct a THn.
+ fCoordBuf = new Int_t[dim];
+}
+
+//______________________________________________________________________________
+THn::~THn()
+{
+ // Destruct a THn
+ delete [] fCoordBuf;
+}
+
+
+//______________________________________________________________________________
+ROOT::THnBaseBinIter* THn::CreateIter(Bool_t respectAxisRange) const
+{
+ // Create an iterator over all bins. Public interface is THnIter.
+ return new THnBinIter(GetNdimensions(), GetListOfAxes(), &GetArray(),
+ respectAxisRange);
+}
+
+//______________________________________________________________________________
+void THn::Sumw2() {
+ // Enable calculation of errors
+ if (!GetCalculateErrors()) {
+ fTsumw2 = 0.;
+ }
+}
+
+
+//______________________________________________________________________________
+THnBase* THn::CloneEmpty(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const
+{
+ // Create a new THn object that is of the same type as *this,
+ // but with dimensions and bins given by axes.
+ // If keepTargetAxis is true, the axes will keep their original xmin / xmax,
+ // else they will be restricted to the range selected (first / last).
+
+ THn* ret = (THn*)IsA()->New();
+ ret->SetNameTitle(name, title);
+
+ TIter iAxis(axes);
+ const TAxis* axis = 0;
+ Int_t pos = 0;
+ Int_t *nbins = new Int_t[axes->GetEntriesFast()];
+ while ((axis = (TAxis*)iAxis())) {
+ TAxis* reqaxis = (TAxis*)axis->Clone();
+ if (!keepTargetAxis && axis->TestBit(TAxis::kAxisRange)) {
+ Int_t binFirst = axis->GetFirst();
+ Int_t binLast = axis->GetLast();
+ Int_t nBins = binLast - binFirst + 1;
+ if (axis->GetXbins()->GetSize()) {
+ // non-uniform bins:
+ reqaxis->Set(nBins, axis->GetXbins()->GetArray() + binFirst - 1);
+ } else {
+ // uniform bins:
+ reqaxis->Set(nBins, axis->GetBinLowEdge(binFirst), axis->GetBinUpEdge(binLast));
+ }
+ reqaxis->ResetBit(TAxis::kAxisRange);
+ }
+
+ nbins[pos] = reqaxis->GetNbins();
+ ret->fAxes.AddAtAndExpand(reqaxis->Clone(), pos++);
+ }
+ ret->fAxes.SetOwner();
+
+ ret->fNdimensions = axes->GetEntriesFast();
+ ret->fCoordBuf = new Int_t[ret->fNdimensions];
+ ret->GetArray().Init(ret->fNdimensions, nbins, true /*addOverflow*/);
+
+ delete [] nbins;
+
+ return ret;
+
+}
+
+//______________________________________________________________________________
+void THn::Reset(Option_t* option /*= ""*/)
+{
+ // Reset the contents of a THn.
+ GetArray().Reset(option);
+ fSumw2.Reset(option);
+}
diff --git a/hist/hist/src/THnBase.cxx b/hist/hist/src/THnBase.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9379e2e44e5e2309570e81fbd28eac4c1c55f960
--- /dev/null
+++ b/hist/hist/src/THnBase.cxx
@@ -0,0 +1,1338 @@
+// @(#)root/hist:$Id$
+// Author: Axel Naumann (2011-12-20)
+
+/*************************************************************************
+ * Copyright (C) 1995-2012, 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 "THnBase.h"
+
+#include "TAxis.h"
+#include "TBrowser.h"
+#include "TError.h"
+#include "TClass.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TH3D.h"
+#include "THn.h"
+#include "THnSparse.h"
+#include "TMath.h"
+#include "TRandom.h"
+#include "TVirtualPad.h"
+
+#include "HFitInterface.h"
+#include "Fit/SparseData.h"
+#include "Math/MinimizerOptions.h"
+#include "Math/WrappedMultiTF1.h"
+
+
+//______________________________________________________________________________
+//
+// Multidimensional histogram base.
+//
+// Defines common functionality and interfaces for THn, THnSparse.
+//
+
+ClassImp(THnBase);
+
+//______________________________________________________________________________
+THnBase::THnBase(const char* name, const char* title, Int_t dim,
+ const Int_t* nbins, const Double_t* xmin, const Double_t* xmax):
+ TNamed(name, title), fNdimensions(dim), fAxes(dim), fEntries(0), fTsumw(0),
+ fTsumw2(-1.), fTsumwx(dim), fTsumwx2(dim),
+ fIntegral(0), fIntegralStatus(kNoInt)
+{
+ // Construct a THnBase with "dim" dimensions,
+ // "nbins" holds the number of bins for each dimension;
+ // "xmin" and "xmax" the minimal and maximal value for each dimension.
+ // The arrays "xmin" and "xmax" can be NULL; in that case SetBinEdges()
+ // must be called for each dimension.
+
+ for (Int_t i = 0; i < fNdimensions; ++i) {
+ TAxis* axis = new TAxis(nbins[i], xmin ? xmin[i] : 0., xmax ? xmax[i] : 1.);
+ fAxes.AddAtAndExpand(axis, i);
+ }
+ SetTitle(title);
+ fAxes.SetOwner();
+}
+
+//______________________________________________________________________________
+THnBase::~THnBase() {
+ // Destruct a THnBase
+ if (fIntegralStatus != kNoInt) delete [] fIntegral;
+}
+
+//______________________________________________________________________________
+TH1* THnBase::CreateHist(const char* name, const char* title,
+ const TObjArray* axes,
+ Bool_t keepTargetAxis ) const {
+ // Create an empty histogram with name and title with a given
+ // set of axes. Create a TH1D/TH2D/TH3D, depending on the number
+ // of elements in axes.
+
+ const int ndim = axes->GetSize();
+
+ TH1* hist = 0;
+ // create hist with dummy axes, we fix them later.
+ if (ndim == 1)
+ hist = new TH1D(name, title, 1, 0., 1.);
+ else if (ndim == 2)
+ hist = new TH2D(name, title, 1, 0., 1., 1, 0., 1.);
+ else if (ndim == 3)
+ hist = new TH3D(name, title, 1, 0., 1., 1, 0., 1., 1, 0., 1.);
+ else {
+ Error("CreateHist", "Cannot create histogram %s with %d dimensions!", name, ndim);
+ return 0;
+ }
+
+ TAxis* hax[3] = {hist->GetXaxis(), hist->GetYaxis(), hist->GetZaxis()};
+ for (Int_t d = 0; d < ndim; ++d) {
+ TAxis* reqaxis = (TAxis*)(*axes)[d];
+ hax[d]->SetTitle(reqaxis->GetTitle());
+ if (!keepTargetAxis && reqaxis->TestBit(TAxis::kAxisRange)) {
+ Int_t binFirst = reqaxis->GetFirst();
+ if (binFirst == 0) binFirst = 1;
+ Int_t binLast = reqaxis->GetLast();
+ Int_t nBins = binLast - binFirst + 1;
+ if (reqaxis->GetXbins()->GetSize()) {
+ // non-uniform bins:
+ hax[d]->Set(nBins, reqaxis->GetXbins()->GetArray() + binFirst - 1);
+ } else {
+ // uniform bins:
+ hax[d]->Set(nBins, reqaxis->GetBinLowEdge(binFirst), reqaxis->GetBinUpEdge(binLast));
+ }
+ } else {
+ if (reqaxis->GetXbins()->GetSize()) {
+ // non-uniform bins:
+ hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXbins()->GetArray());
+ } else {
+ // uniform bins:
+ hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXmin(), reqaxis->GetXmax());
+ }
+ }
+ }
+
+ hist->Rebuild();
+
+ return hist;
+}
+
+//______________________________________________________________________________
+THnBase* THnBase::CreateHnAny(const char* name, const char* title,
+ const TH1* h, Bool_t sparse, Int_t chunkSize)
+{
+ // Create a THn / THnSparse object from a histogram deriving from TH1.
+
+ // Get the dimension of the TH1
+ int ndim = h->GetDimension();
+
+ // Axis properties
+ int nbins[3] = {0,0,0};
+ double minRange[3] = {0.,0.,0.};
+ double maxRange[3] = {0.,0.,0.};
+ TAxis* axis[3] = { h->GetXaxis(), h->GetYaxis(), h->GetZaxis() };
+ for (int i = 0; i < ndim; ++i) {
+ nbins[i] = axis[i]->GetNbins();
+ minRange[i] = axis[i]->GetXmin();
+ maxRange[i] = axis[i]->GetXmax();
+ }
+
+ // Create the corresponding THnSparse, depending on the storage
+ // type of the TH1. The class name will be "TH??\0" where the first
+ // ? is 1,2 or 3 and the second ? indicates the stograge as C, S,
+ // I, F or D.
+ THnBase* s = 0;
+ const char* cname( h->ClassName() );
+ if (cname[0] == 'T' && cname[1] == 'H'
+ && cname[2] >= '1' && cname[2] <= '3' && cname[4] == 0) {
+
+#define R__THNBCASE(TAG) \
+ if (sparse) { \
+ s = new _NAME2_(THnSparse,TAG)(name, title, ndim, nbins, \
+ minRange, maxRange, chunkSize); \
+ } else { \
+ s = new _NAME2_(THn,TAG)(name, title, ndim, nbins, \
+ minRange, maxRange); \
+ } \
+ break;
+
+ switch (cname[3]) {
+ case 'F': R__THNBCASE(F);
+ case 'D': R__THNBCASE(D);
+ case 'I': R__THNBCASE(I);
+ case 'S': R__THNBCASE(S);
+ case 'C': R__THNBCASE(C);
+ }
+#undef R__THNBCASE
+ }
+ if (!s) {
+ ::Warning("THnSparse::CreateHnAny", "Unknown Type of Histogram");
+ return 0;
+ }
+
+ for (int i = 0; i < ndim; ++i) {
+ s->GetAxis(i)->SetTitle(axis[i]->GetTitle());
+ }
+
+ // Get the array to know the number of entries of the TH1
+ const TArray *array = dynamic_cast<const TArray*>(h);
+ if (!array) {
+ ::Warning("THnSparse::CreateHnAny", "Unknown Type of Histogram");
+ return 0;
+ }
+
+ s->Add(h);
+ return s;
+}
+
+
+//______________________________________________________________________________
+void THnBase::Add(const TH1* hist, Double_t c /*=1.*/)
+{
+ // Fill the THnBase with the bins of hist that have content
+ // or error != 0.
+ Long64_t nbins = hist->GetNbinsX() + 2;
+ if (hist->GetDimension() >= 2) nbins *= hist->GetNbinsY() + 2;
+ if (hist->GetDimension() >= 3) nbins *= hist->GetNbinsZ() + 2;
+ int x[3] = {0,0,0};
+ for (int i = 0; i < nbins; ++i) {
+ double value = hist->GetBinContent(i);
+ double error = hist->GetBinError(i);
+ if (!value && !error) continue;
+ hist->GetBinXYZ(i, x[0], x[1], x[2]);
+ SetBinContent(x, value * c);
+ SetBinError(x, error * c);
+ }
+}
+
+//______________________________________________________________________________
+TFitResultPtr THnBase::Fit(TF1 *f ,Option_t *option ,Option_t *goption)
+{
+// Fit a THnSparse with function f
+//
+// since the data is sparse by default a likelihood fit is performed
+// merging all the regions with empty bins for betetr performance efficiency
+//
+// Since the THnSparse is not drawn no graphics options are passed
+// Here is the list of possible options
+//
+// = "I" Use integral of function in bin instead of value at bin center
+// = "X" Use chi2 method (default is log-likelihood method)
+// = "U" Use a User specified fitting algorithm (via SetFCN)
+// = "Q" Quiet mode (minimum printing)
+// = "V" Verbose mode (default is between Q and V)
+// = "E" Perform better Errors estimation using Minos technique
+// = "B" Use this option when you want to fix one or more parameters
+// and the fitting function is like "gaus", "expo", "poln", "landau".
+// = "M" More. Improve fit results
+// = "R" Use the Range specified in the function range
+
+
+ Foption_t fitOption;
+
+ if (!TH1::FitOptionsMake(option,fitOption)) return 0;
+
+ // The function used to fit cannot be stored in a THnSparse. It
+ // cannot be drawn either. Perhaps in the future.
+ fitOption.Nostore = true;
+ // Use likelihood fit if not specified
+ if (!fitOption.Chi2) fitOption.Like = true;
+ // create range and minimizer options with default values
+ ROOT::Fit::DataRange range(GetNdimensions());
+ for ( int i = 0; i < GetNdimensions(); ++i ) {
+ TAxis *axis = GetAxis(i);
+ range.AddRange(i, axis->GetXmin(), axis->GetXmax());
+ }
+ ROOT::Math::MinimizerOptions minOption;
+
+ return ROOT::Fit::FitObject(this, f , fitOption , minOption, goption, range);
+}
+
+//______________________________________________________________________________
+void THnBase::GetRandom(Double_t *rand, Bool_t subBinRandom /* = kTRUE */)
+{
+ // Generate an n-dimensional random tuple based on the histogrammed
+ // distribution. If subBinRandom, the returned tuple will be additionally
+ // randomly distributed within the randomized bin, using a flat
+ // distribution.
+
+ // check whether the integral array is valid
+ if (fIntegralStatus != kValidInt)
+ ComputeIntegral();
+
+ // generate a random bin
+ Double_t p = gRandom->Rndm();
+ Long64_t idx = TMath::BinarySearch(GetNbins() + 1, fIntegral, p);
+ const Int_t nStaticBins = 40;
+ Int_t bin[nStaticBins];
+ Int_t* pBin = bin;
+ if (GetNdimensions() > nStaticBins) {
+ pBin = new Int_t[GetNdimensions()];
+ }
+ GetBinContent(idx, pBin);
+
+ // convert bin coordinates to real values
+ for (Int_t i = 0; i < fNdimensions; i++) {
+ rand[i] = GetAxis(i)->GetBinCenter(pBin[i]);
+
+ // randomize the vector withing a bin
+ if (subBinRandom)
+ rand[i] += (gRandom->Rndm() - 0.5) * GetAxis(i)->GetBinWidth(pBin[i]);
+ }
+ if (pBin != bin) {
+ delete [] pBin;
+ }
+
+ return;
+}
+
+//______________________________________________________________________________
+Bool_t THnBase::IsInRange(Int_t *coord) const
+{
+ // Check whether bin coord is in range, as defined by TAxis::SetRange().
+ // Currently, TAxis::SetRange() does not allow to select all but over- and
+ // underflow bins (it instead resets the axis to "no range selected").
+ // Instead, simply call
+ // TAxis* axis12 = hsparse.GetAxis(12);
+ // axis12->SetRange(1, axis12->GetNbins());
+ // axis12->SetBit(TAxis::kAxisRange);
+ // to deselect the under- and overflow bins in the 12th dimension.
+
+ Int_t min = 0;
+ Int_t max = 0;
+ for (Int_t i = 0; i < fNdimensions; ++i) {
+ TAxis *axis = GetAxis(i);
+ if (!axis->TestBit(TAxis::kAxisRange)) continue;
+ min = axis->GetFirst();
+ max = axis->GetLast();
+ if (min == 0 && max == 0) {
+ // special case where TAxis::SetBit(kAxisRange) and
+ // over- and underflow bins are de-selected.
+ // first and last are == 0 due to axis12->SetRange(1, axis12->GetNbins());
+ min = 1;
+ max = axis->GetNbins();
+ }
+ if (coord[i] < min || coord[i] > max)
+ return kFALSE;
+ }
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+TObject* THnBase::ProjectionAny(Int_t ndim, const Int_t* dim,
+ Bool_t wantNDim,
+ Option_t* option /*= ""*/) const
+{
+ // Project all bins into a ndim-dimensional THn / THnSparse (whatever
+ // *this is) or if (ndim < 4 and !wantNDim) a TH1/2/3 histogram,
+ // keeping only axes in dim (specifying ndim dimensions).
+ // If "option" contains "E" errors will be calculated.
+ // "A" ranges of the taget axes will be ignored.
+ // "O" original axis range of the taget axes will be
+ // kept, but only bins inside the selected range
+ // will be filled.
+
+ TString name(GetName());
+ name +="_proj";
+
+ for (Int_t d = 0; d < ndim; ++d) {
+ name += "_";
+ name += dim[d];
+ }
+
+ TString title(GetTitle());
+ Ssiz_t posInsert = title.First(';');
+ if (posInsert == kNPOS) {
+ title += " projection ";
+ for (Int_t d = 0; d < ndim; ++d)
+ title += GetAxis(dim[d])->GetTitle();
+ } else {
+ for (Int_t d = ndim - 1; d >= 0; --d) {
+ title.Insert(posInsert, GetAxis(d)->GetTitle());
+ if (d)
+ title.Insert(posInsert, ", ");
+ }
+ title.Insert(posInsert, " projection ");
+ }
+
+ TObjArray newaxes(ndim);
+ for (Int_t d = 0; d < ndim; ++d) {
+ newaxes.AddAt(GetAxis(dim[d]),d);
+ }
+
+ THnBase* hn = 0;
+ TH1* hist = 0;
+ TObject* ret = 0;
+
+ Bool_t* hadRange = 0;
+ Bool_t ignoreTargetRange = (option && (strchr(option, 'A') || strchr(option, 'a')));
+ Bool_t keepTargetAxis = ignoreTargetRange || (option && (strchr(option, 'O') || strchr(option, 'o')));
+ if (ignoreTargetRange) {
+ hadRange = new Bool_t[ndim];
+ for (Int_t d = 0; d < ndim; ++d){
+ TAxis *axis = GetAxis(dim[d]);
+ hadRange[d] = axis->TestBit(TAxis::kAxisRange);
+ axis->SetBit(TAxis::kAxisRange, kFALSE);
+ }
+ }
+
+ if (wantNDim)
+ ret = hn = CloneEmpty(name, title, &newaxes, keepTargetAxis);
+ else
+ ret = hist = CreateHist(name, title, &newaxes, keepTargetAxis);
+
+ if (keepTargetAxis) {
+ // make the whole axes visible, i.e. unset the range
+ if (wantNDim) {
+ for (Int_t d = 0; d < ndim; ++d) {
+ hn->GetAxis(d)->SetRange(0, 0);
+ }
+ } else {
+ hist->GetXaxis()->SetRange(0, 0);
+ hist->GetYaxis()->SetRange(0, 0);
+ hist->GetZaxis()->SetRange(0, 0);
+ }
+ }
+
+ Bool_t haveErrors = GetCalculateErrors();
+ Bool_t wantErrors = haveErrors || (option && (strchr(option, 'E') || strchr(option, 'e')));
+
+ Int_t* bins = new Int_t[ndim];
+ Long64_t myLinBin = 0;
+
+ THnIter iter(this, kTRUE /*use axis range*/);
+
+ while ((myLinBin = iter.Next()) >= 0) {
+ Double_t v = GetBinContent(myLinBin);
+
+ for (Int_t d = 0; d < ndim; ++d) {
+ bins[d] = iter.GetCoord(dim[d]);
+ if (!keepTargetAxis && GetAxis(dim[d])->TestBit(TAxis::kAxisRange)) {
+ bins[d] -= GetAxis(dim[d])->GetFirst() - 1;
+ }
+ }
+
+ Long64_t targetLinBin = -1;
+ if (!wantNDim) {
+ if (ndim == 1) targetLinBin = bins[0];
+ else if (ndim == 2) targetLinBin = hist->GetBin(bins[0], bins[1]);
+ else if (ndim == 3) targetLinBin = hist->GetBin(bins[0], bins[1], bins[2]);
+ } else {
+ targetLinBin = hn->GetBin(bins, kTRUE /*allocate*/);
+ }
+
+ if (wantErrors) {
+ Double_t err2 = 0.;
+ if (haveErrors) {
+ err2 = GetBinError2(myLinBin);
+ } else {
+ err2 = v;
+ }
+ if (wantNDim) {
+ hn->AddBinError2(targetLinBin, err2);
+ } else {
+ Double_t preverr = hist->GetBinError(targetLinBin);
+ hist->SetBinError(targetLinBin, TMath::Sqrt(preverr * preverr + err2));
+ }
+ }
+
+ // only _after_ error calculation, or sqrt(v) is taken into account!
+ if (wantNDim)
+ hn->AddBinContent(targetLinBin, v);
+ else
+ hist->AddBinContent(targetLinBin, v);
+ }
+
+ delete [] bins;
+
+ if (wantNDim) {
+ hn->SetEntries(fEntries);
+ } else {
+ if (!iter.HaveSkippedBin()) {
+ hist->SetEntries(fEntries);
+ } else {
+ // re-compute the entries
+ // in case of error calculation (i.e. when Sumw2() is set)
+ // use the effective entries for the entries
+ // since this is the only way to estimate them
+ hist->ResetStats();
+ Double_t entries = hist->GetEffectiveEntries();
+ if (!wantErrors) {
+ // to avoid numerical rounding
+ entries = TMath::Floor(entries + 0.5);
+ }
+ hist->SetEntries(entries);
+ }
+ }
+
+ if (hadRange) {
+ // reset kAxisRange bit:
+ for (Int_t d = 0; d < ndim; ++d)
+ GetAxis(dim[d])->SetBit(TAxis::kAxisRange, hadRange[d]);
+
+ delete [] hadRange;
+ }
+
+ return ret;
+}
+
+//______________________________________________________________________________
+void THnBase::Scale(Double_t c)
+{
+ // Scale contents and errors of this histogram by c:
+ // this = this * c
+ // It does not modify the histogram's number of entries.
+
+
+ Double_t nEntries = GetEntries();
+ // Scale the contents & errors
+ Bool_t haveErrors = GetCalculateErrors();
+ Long64_t i = 0;
+ THnIter iter(this);
+ while ((i = iter.Next()) >= 0) {
+ // Get the content of the bin from the current histogram
+ Double_t v = GetBinContent(i);
+ SetBinContent(i, c * v);
+ if (haveErrors) {
+ Double_t err2 = GetBinError2(i);
+ SetBinError2(i, c * c * err2);
+ }
+ }
+ SetEntries(nEntries);
+}
+
+//______________________________________________________________________________
+void THnBase::AddInternal(const THnBase* h, Double_t c, Bool_t rebinned)
+{
+ // Add() implementation for both rebinned histograms and those with identical
+ // binning. See THnBase::Add().
+
+ if (fNdimensions != h->GetNdimensions()) {
+ Warning("RebinnedAdd", "Different number of dimensions, cannot carry out operation on the histograms");
+ return;
+ }
+
+ // Trigger error calculation if h has it
+ if (!GetCalculateErrors() && h->GetCalculateErrors())
+ Sumw2();
+ Bool_t haveErrors = GetCalculateErrors();
+
+ Double_t* x = 0;
+ if (rebinned) {
+ x = new Double_t[fNdimensions];
+ }
+ Int_t* coord = new Int_t[fNdimensions];
+
+ // Expand the exmap if needed, to reduce collisions
+ Long64_t numTargetBins = GetNbins() + h->GetNbins();
+ Reserve(numTargetBins);
+
+ Long64_t i = 0;
+ THnIter iter(h);
+ // Add to this whatever is found inside the other histogram
+ while ((i = iter.Next(coord)) >= 0) {
+ // Get the content of the bin from the second histogram
+ Double_t v = h->GetBinContent(i);
+
+ Long64_t mybinidx = -1;
+ if (rebinned) {
+ // Get the bin center given a coord
+ for (Int_t j = 0; j < fNdimensions; ++j)
+ x[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
+
+ mybinidx = GetBin(x, kTRUE /* allocate*/);
+ } else {
+ mybinidx = GetBin(coord, kTRUE /*allocate*/);
+ }
+
+ if (haveErrors) {
+ Double_t err2 = h->GetBinError2(i) * c * c;
+ AddBinError2(mybinidx, err2);
+ }
+ // only _after_ error calculation, or sqrt(v) is taken into account!
+ AddBinContent(mybinidx, c * v);
+ }
+
+ delete [] coord;
+ delete [] x;
+
+ Double_t nEntries = GetEntries() + c * h->GetEntries();
+ SetEntries(nEntries);
+}
+
+//______________________________________________________________________________
+void THnBase::Add(const THnBase* h, Double_t c)
+{
+ // Add contents of h scaled by c to this histogram:
+ // this = this + c * h
+ // Note that if h has Sumw2 set, Sumw2 is automatically called for this
+ // if not already set.
+
+ // Check consistency of the input
+ if (!CheckConsistency(h, "Add")) return;
+
+ AddInternal(h, c, kFALSE);
+}
+
+//______________________________________________________________________________
+void THnBase::RebinnedAdd(const THnBase* h, Double_t c)
+{
+ // Add contents of h scaled by c to this histogram:
+ // this = this + c * h
+ // Note that if h has Sumw2 set, Sumw2 is automatically called for this
+ // if not already set.
+ // In contrast to Add(), RebinnedAdd() does not require consistent binning of
+ // this and h; instead, each bin's center is used to determine the target bin.
+
+ AddInternal(h, c, kTRUE);
+}
+
+
+//______________________________________________________________________________
+Long64_t THnBase::Merge(TCollection* list)
+{
+ // Merge this with a list of THnBase's. All THnBase's provided
+ // in the list must have the same bin layout!
+
+ if (!list) return 0;
+ if (list->IsEmpty()) return (Long64_t)GetEntries();
+
+ Long64_t sumNbins = GetNbins();
+ TIter iter(list);
+ const TObject* addMeObj = 0;
+ while ((addMeObj = iter())) {
+ const THnBase* addMe = dynamic_cast<const THnBase*>(addMeObj);
+ if (addMe) {
+ sumNbins += addMe->GetNbins();
+ }
+ }
+ Reserve(sumNbins);
+
+ iter.Reset();
+ while ((addMeObj = iter())) {
+ const THnBase* addMe = dynamic_cast<const THnBase*>(addMeObj);
+ if (!addMe)
+ Error("Merge", "Object named %s is not THnBase! Skipping it.",
+ addMeObj->GetName());
+ else
+ Add(addMe);
+ }
+ return (Long64_t)GetEntries();
+}
+
+
+//______________________________________________________________________________
+void THnBase::Multiply(const THnBase* h)
+{
+ // Multiply this histogram by histogram h
+ // this = this * h
+ // Note that if h has Sumw2 set, Sumw2 is automatically called for this
+ // if not already set.
+
+ // Check consistency of the input
+ if(!CheckConsistency(h, "Multiply"))return;
+
+ // Trigger error calculation if h has it
+ Bool_t wantErrors = kFALSE;
+ if (GetCalculateErrors() || h->GetCalculateErrors())
+ wantErrors = kTRUE;
+
+ if (wantErrors) Sumw2();
+
+ Double_t nEntries = GetEntries();
+ // Now multiply the contents: in this case we have the intersection of the sets of bins
+ Int_t* coord = new Int_t[fNdimensions];
+ Long64_t i = 0;
+ THnIter iter(this);
+ // Add to this whatever is found inside the other histogram
+ while ((i = iter.Next(coord)) >= 0) {
+ // Get the content of the bin from the current histogram
+ Double_t v1 = GetBinContent(i);
+ // Now look at the bin with the same coordinates in h
+ Long64_t idxh = h->GetBin(coord);
+ Double_t v2 = 0.;
+ if (idxh >= 0) v2 = h->GetBinContent(idxh);
+ SetBinContent(i, v1 * v2);
+ if (wantErrors) {
+ Double_t err1 = GetBinError(i) * v2;
+ Double_t err2 = 0.;
+ if (idxh >= 0) err2 = h->GetBinError(idxh) * v1;
+ SetBinError(i, TMath::Sqrt((err2 * err2 + err1 * err1)));
+ }
+ }
+ SetEntries(nEntries);
+
+ delete [] coord;
+}
+
+//______________________________________________________________________________
+void THnBase::Multiply(TF1* f, Double_t c)
+{
+ // Performs the operation: this = this*c*f1
+ // if errors are defined, 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
+ // calling Multiply()
+
+ Int_t* coord = new Int_t[fNdimensions];
+ Double_t* x = new Double_t[fNdimensions];
+
+ Bool_t wantErrors = GetCalculateErrors();
+ if (wantErrors) Sumw2();
+
+ Long64_t i = 0;
+ THnIter iter(this);
+ // Add to this whatever is found inside the other histogram
+ while ((i = iter.Next(coord)) >= 0) {
+ Double_t value = GetBinContent(i);
+
+ // Get the bin coordinates given an index array
+ for (Int_t j = 0; j < fNdimensions; ++j)
+ x[j] = GetAxis(j)->GetBinCenter(coord[j]);
+
+ if (!f->IsInside(x))
+ continue;
+ TF1::RejectPoint(kFALSE);
+
+ // Evaluate function at points
+ Double_t fvalue = f->EvalPar(x, NULL);
+
+ SetBinContent(i, c * fvalue * value);
+ if (wantErrors) {
+ Double_t error = GetBinError(i);
+ SetBinError(i, c * fvalue * error);
+ }
+ }
+
+ delete [] x;
+ delete [] coord;
+}
+
+//______________________________________________________________________________
+void THnBase::Divide(const THnBase *h)
+{
+ // Divide this histogram by h
+ // this = this/(h)
+ // Note that if h has Sumw2 set, Sumw2 is automatically called for
+ // this if not already set.
+ // The resulting errors are calculated assuming uncorrelated content.
+
+ // Check consistency of the input
+ if (!CheckConsistency(h, "Divide"))return;
+
+ // Trigger error calculation if h has it
+ Bool_t wantErrors=GetCalculateErrors();
+ if (!GetCalculateErrors() && h->GetCalculateErrors())
+ wantErrors=kTRUE;
+
+ // Remember original histogram statistics
+ Double_t nEntries = fEntries;
+
+ if (wantErrors) Sumw2();
+ Bool_t didWarn = kFALSE;
+
+ // Now divide the contents: also in this case we have the intersection of the sets of bins
+ Int_t* coord = new Int_t[fNdimensions];
+ Long64_t i = 0;
+ THnIter iter(this);
+ // Add to this whatever is found inside the other histogram
+ while ((i = iter.Next(coord)) >= 0) {
+ // Get the content of the bin from the first histogram
+ Double_t v1 = GetBinContent(i);
+ // Now look at the bin with the same coordinates in h
+ Long64_t hbin = h->GetBin(coord);
+ Double_t v2 = h->GetBinContent(hbin);
+ if (!v2) {
+ v1 = 0.;
+ v2 = 1.;
+ if (!didWarn) {
+ Warning("Divide(h)", "Histogram h has empty bins - division by zero! Setting bin to 0.");
+ didWarn = kTRUE;
+ }
+ }
+ SetBinContent(i, v1 / v2);
+ if (wantErrors) {
+ Double_t err1 = GetBinError(i) * v2;
+ Double_t err2 = h->GetBinError(hbin) * v1;
+ Double_t b22 = v2 * v2;
+ Double_t err = (err1 * err1 + err2 * err2) / (b22 * b22);
+ SetBinError2(i, err);
+ }
+ }
+ delete [] coord;
+ SetEntries(nEntries);
+}
+
+//______________________________________________________________________________
+void THnBase::Divide(const THnBase *h1, const THnBase *h2, Double_t c1, Double_t c2, Option_t *option)
+{
+ // Replace contents of this histogram by multiplication of h1 by h2
+ // this = (c1*h1)/(c2*h2)
+ // 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 content.
+ // 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.
+
+
+ TString opt = option;
+ opt.ToLower();
+ Bool_t binomial = kFALSE;
+ if (opt.Contains("b")) binomial = kTRUE;
+
+ // Check consistency of the input
+ if (!CheckConsistency(h1, "Divide") || !CheckConsistency(h2, "Divide"))return;
+ if (!c2) {
+ Error("Divide","Coefficient of dividing histogram cannot be zero");
+ return;
+ }
+
+ Reset();
+
+ // Trigger error calculation if h1 or h2 have it
+ if (!GetCalculateErrors() && (h1->GetCalculateErrors()|| h2->GetCalculateErrors() != 0))
+ Sumw2();
+
+ // Count filled bins
+ Long64_t nFilledBins=0;
+
+ // Now divide the contents: we have the intersection of the sets of bins
+
+ Int_t* coord = new Int_t[fNdimensions];
+ memset(coord, 0, sizeof(Int_t) * fNdimensions);
+ Bool_t didWarn = kFALSE;
+
+ Long64_t i = 0;
+ THnIter iter(h1);
+ // Add to this whatever is found inside the other histogram
+ while ((i = iter.Next(coord)) >= 0) {
+ // Get the content of the bin from the first histogram
+ Double_t v1 = h1->GetBinContent(i);
+ // Now look at the bin with the same coordinates in h2
+ Long64_t h2bin = h2->GetBin(coord);
+ Double_t v2 = h2->GetBinContent(h2bin);
+ if (!v2) {
+ v1 = 0.;
+ v2 = 1.;
+ if (!didWarn) {
+ Warning("Divide(h1, h2)", "Histogram h2 has empty bins - division by zero! Setting bin to 0.");
+ didWarn = kTRUE;
+ }
+ }
+ nFilledBins++;
+ Long64_t myBin = GetBin(coord);
+ SetBinContent(myBin, c1 * v1 / c2 / v2);
+ if(GetCalculateErrors()){
+ Double_t err1 = h1->GetBinError(i);
+ Double_t err2 = h2->GetBinError(h2bin);
+ Double_t errSq = 0.;
+ if (binomial) {
+ if (v1 != v2) {
+ Double_t w = v1 / v2;
+ err2 *= w;
+ errSq = TMath::Abs( ( (1. - 2.*w) * err1 * err1 + err2 * err2 ) / (v2 * v2) );
+ }
+ } else {
+ c1 *= c1;
+ c2 *= c2;
+ Double_t b22 = v2 * v2 * c2;
+ err1 *= v2;
+ err2 *= v1;
+ errSq = c1 * c2 * (err1 * err1 + err2 * err2) / (b22 * b22);
+ }
+ SetBinError2(myBin, errSq);
+ }
+ }
+
+ delete [] coord;
+ SetFilledBins(nFilledBins);
+
+ // Set as entries in the result histogram the entries in the numerator
+ SetEntries(h1->GetEntries());
+}
+
+//______________________________________________________________________________
+Bool_t THnBase::CheckConsistency(const THnBase *h, const char *tag) const
+{
+ // Consistency check on (some of) the parameters of two histograms (for operations).
+
+ if (fNdimensions != h->GetNdimensions()) {
+ Warning(tag, "Different number of dimensions, cannot carry out operation on the histograms");
+ return kFALSE;
+ }
+ for (Int_t dim = 0; dim < fNdimensions; dim++){
+ if (GetAxis(dim)->GetNbins() != h->GetAxis(dim)->GetNbins()) {
+ Warning(tag, "Different number of bins on axis %i, cannot carry out operation on the histograms", dim);
+ return kFALSE;
+ }
+ }
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+void THnBase::SetBinEdges(Int_t idim, const Double_t* bins)
+{
+ // Set the axis # of bins and bin limits on dimension idim
+
+ TAxis* axis = (TAxis*) fAxes[idim];
+ axis->Set(axis->GetNbins(), bins);
+}
+
+//______________________________________________________________________________
+void THnBase::SetTitle(const char *title)
+{
+ // Change (i.e. set) the title.
+ //
+ // If title is in the form "stringt;string0;string1;string2 ..."
+ // the histogram title is set to stringt, the title of axis0 to string0,
+ // of axis1 to string1, of axis2 to string2, etc, just like it is done
+ // for TH1/TH2/TH3.
+ // To insert the character ";" in one of the titles, one should use "#;"
+ // or "#semicolon".
+
+ fTitle = title;
+ fTitle.ReplaceAll("#;",2,"#semicolon",10);
+
+ Int_t endHistTitle = fTitle.First(';');
+ if (endHistTitle >= 0) {
+ // title contains a ';' so parse the axis titles
+ Int_t posTitle = endHistTitle + 1;
+ Int_t lenTitle = fTitle.Length();
+ Int_t dim = 0;
+ while (posTitle > 0 && posTitle < lenTitle && dim < fNdimensions){
+ Int_t endTitle = fTitle.Index(";", posTitle);
+ TString axisTitle = fTitle(posTitle, endTitle - posTitle);
+ axisTitle.ReplaceAll("#semicolon", 10, ";", 1);
+ GetAxis(dim)->SetTitle(axisTitle);
+ dim++;
+ if (endTitle > 0)
+ posTitle = endTitle + 1;
+ else
+ posTitle = -1;
+ }
+ // Remove axis titles from histogram title
+ fTitle.Remove(endHistTitle, lenTitle - endHistTitle);
+ }
+
+ fTitle.ReplaceAll("#semicolon", 10, ";", 1);
+
+}
+
+//______________________________________________________________________________
+THnBase* THnBase::RebinBase(Int_t group) const
+{
+ // Combine the content of "group" neighboring bins into
+ // a new bin and return the resulting THnBase.
+ // For group=2 and a 3 dimensional histogram, all "blocks"
+ // of 2*2*2 bins will be put into a bin.
+
+ Int_t* ngroup = new Int_t[GetNdimensions()];
+ for (Int_t d = 0; d < GetNdimensions(); ++d)
+ ngroup[d] = group;
+ THnBase* ret = RebinBase(ngroup);
+ delete [] ngroup;
+ return ret;
+}
+
+//______________________________________________________________________________
+THnBase* THnBase::RebinBase(const Int_t* group) const
+{
+ // Combine the content of "group" neighboring bins for each dimension
+ // into a new bin and return the resulting THnBase.
+ // For group={2,1,1} and a 3 dimensional histogram, pairs of x-bins
+ // will be grouped.
+
+ Int_t ndim = GetNdimensions();
+ TString name(GetName());
+ for (Int_t d = 0; d < ndim; ++d)
+ name += Form("_%d", group[d]);
+
+
+ TString title(GetTitle());
+ Ssiz_t posInsert = title.First(';');
+ if (posInsert == kNPOS) {
+ title += " rebin ";
+ for (Int_t d = 0; d < ndim; ++d)
+ title += Form("{%d}", group[d]);
+ } else {
+ for (Int_t d = ndim - 1; d >= 0; --d)
+ title.Insert(posInsert, Form("{%d}", group[d]));
+ title.Insert(posInsert, " rebin ");
+ }
+
+ TObjArray newaxes(ndim);
+ newaxes.SetOwner();
+ for (Int_t d = 0; d < ndim; ++d) {
+ newaxes.AddAt(GetAxis(d)->Clone(),d);
+ if (group[d] > 1) {
+ TAxis* newaxis = (TAxis*) newaxes.At(d);
+ Int_t newbins = (newaxis->GetNbins() + group[d] - 1) / group[d];
+ if (newaxis->GetXbins() && newaxis->GetXbins()->GetSize()) {
+ // variable bins
+ Double_t *edges = new Double_t[newbins + 1];
+ for (Int_t i = 0; i < newbins + 1; ++i)
+ if (group[d] * i <= newaxis->GetNbins())
+ edges[i] = newaxis->GetXbins()->At(group[d] * i);
+ else edges[i] = newaxis->GetXmax();
+ newaxis->Set(newbins, edges);
+ delete [] edges;
+ } else {
+ newaxis->Set(newbins, newaxis->GetXmin(), newaxis->GetXmax());
+ }
+ }
+ }
+
+ THnBase* h = CloneEmpty(name.Data(), title.Data(), &newaxes, kTRUE);
+ Bool_t haveErrors = GetCalculateErrors();
+ Bool_t wantErrors = haveErrors;
+
+ Int_t* bins = new Int_t[ndim];
+ Int_t* coord = new Int_t[fNdimensions];
+
+ Long64_t i = 0;
+ THnIter iter(this);
+ while ((i = iter.Next(coord)) >= 0) {
+ Double_t v = GetBinContent(i);
+ for (Int_t d = 0; d < ndim; ++d) {
+ bins[d] = TMath::CeilNint( (double) coord[d]/group[d] );
+ }
+ Long64_t idxh = h->GetBin(bins, kTRUE /*allocate*/);
+
+ if (wantErrors) {
+ Double_t err2 = 0.;
+ if (haveErrors) {
+ err2 = GetBinError2(i);
+ } else err2 = v;
+ h->AddBinError2(idxh, err2);
+ }
+
+ // only _after_ error calculation, or sqrt(v) is taken into account!
+ h->AddBinContent(idxh, v);
+ }
+
+ delete [] bins;
+ delete [] coord;
+
+ h->SetEntries(fEntries);
+
+ return h;
+
+}
+
+//______________________________________________________________________________
+void THnBase::ResetBase(Option_t * /*option = ""*/)
+{
+ // Clear the histogram
+ fEntries = 0.;
+ fTsumw = 0.;
+ fTsumw2 = -1.;
+ if (fIntegralStatus != kNoInt) {
+ delete [] fIntegral;
+ fIntegralStatus = kNoInt;
+ }
+}
+
+//______________________________________________________________________________
+Double_t THnBase::ComputeIntegral()
+{
+ // Calculate the integral of the histogram
+
+ // delete old integral
+ if (fIntegralStatus != kNoInt) {
+ delete [] fIntegral;
+ fIntegralStatus = kNoInt;
+ }
+
+ // check number of bins
+ if (GetNbins() == 0) {
+ Error("ComputeIntegral", "The histogram must have at least one bin.");
+ return 0.;
+ }
+
+ // allocate integral array
+ fIntegral = new Double_t [GetNbins() + 1];
+ fIntegral[0] = 0.;
+
+ // fill integral array with contents of regular bins (non over/underflow)
+ Int_t* coord = new Int_t[fNdimensions];
+ Long64_t i = 0;
+ THnIter iter(this);
+ while ((i = iter.Next(coord)) >= 0) {
+ Double_t v = GetBinContent(i);
+
+ // check whether the bin is regular
+ bool regularBin = true;
+ for (Int_t dim = 0; dim < fNdimensions; dim++) {
+ if (coord[dim] < 1 || coord[dim] > GetAxis(dim)->GetNbins()) {
+ regularBin = false;
+ break;
+ }
+ }
+
+ // if outlayer, count it with zero weight
+ if (!regularBin) v = 0.;
+
+ fIntegral[i + 1] = fIntegral[i] + v;
+ }
+ delete [] coord;
+
+ // check sum of weights
+ if (fIntegral[GetNbins()] == 0.) {
+ Error("ComputeIntegral", "No hits in regular bins (non over/underflow).");
+ delete [] fIntegral;
+ return 0.;
+ }
+
+ // normalize the integral array
+ for (Long64_t j = 0; j <= GetNbins(); ++j)
+ fIntegral[j] = fIntegral[j] / fIntegral[GetNbins()];
+
+ // set status to valid
+ fIntegralStatus = kValidInt;
+ return fIntegral[GetNbins()];
+}
+
+//______________________________________________________________________________
+void THnBase::PrintBin(Long64_t idx, Option_t* options) const
+{
+ // Print bin with linex index "idx".
+ // For valid options see PrintBin(Long64_t idx, Int_t* bin, Option_t* options).
+ Int_t* coord = new Int_t[fNdimensions];
+ PrintBin(idx, coord, options);
+ delete [] coord;
+}
+
+//______________________________________________________________________________
+Bool_t THnBase::PrintBin(Long64_t idx, Int_t* bin, Option_t* options) const
+{
+ // Print one bin. If "idx" is != -1 use that to determine the bin,
+ // otherwise (if "idx" == -1) use the coordinate in "bin".
+ // If "options" contains:
+ // '0': only print bins with an error or content != 0
+ // Return whether the bin was printed (depends on options)
+
+ Double_t v = -42;
+ if (idx == -1) {
+ idx = GetBin(bin);
+ v = GetBinContent(idx);
+ } else {
+ v = GetBinContent(idx, bin);
+ }
+
+ Double_t err = 0.;
+ if (GetCalculateErrors()) {
+ if (idx != -1) {
+ err = GetBinError(idx);
+ }
+ }
+
+ if (v == 0. && err == 0. && options && strchr(options, '0')) {
+ // suppress zeros, and we have one.
+ return kFALSE;
+ }
+
+ TString coord;
+ for (Int_t dim = 0; dim < fNdimensions; ++dim) {
+ coord += bin[dim];
+ coord += ',';
+ }
+ coord.Remove(coord.Length() - 1);
+
+ if (GetCalculateErrors()) {
+ Printf("Bin at (%s) = %g (+/- %g)", coord.Data(), v, err);
+ } else {
+ Printf("Bin at (%s) = %g", coord.Data(), v);
+ }
+
+ return kTRUE;
+}
+
+//______________________________________________________________________________
+void THnBase::PrintEntries(Long64_t from /*=0*/, Long64_t howmany /*=-1*/,
+ Option_t* options /*=0*/) const
+{
+ // Print "howmany" entries starting at "from". If "howmany" is -1, print all.
+ // If "options" contains:
+ // 'x': print in the order of axis bins, i.e. (0,0,...,0), (0,0,...,1),...
+ // '0': only print bins with content != 0
+
+ if (from < 0) from = 0;
+ if (howmany == -1) howmany = GetNbins();
+
+ Int_t* bin = new Int_t[fNdimensions];
+
+ if (options && (strchr(options, 'x') || strchr(options, 'X'))) {
+ Int_t* nbins = new Int_t[fNdimensions];
+ for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
+ nbins[dim] = GetAxis(dim)->GetNbins();
+ bin[dim] = from % nbins[dim];
+ from /= nbins[dim];
+ }
+
+ for (Long64_t i = 0; i < howmany; ++i) {
+ if (!PrintBin(-1, bin, options))
+ ++howmany;
+ // Advance to next bin:
+ ++bin[fNdimensions - 1];
+ for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
+ if (bin[dim] >= nbins[dim]) {
+ bin[dim] = 0;
+ if (dim > 0) {
+ ++bin[dim - 1];
+ } else {
+ howmany = -1; // aka "global break"
+ }
+ }
+ }
+ }
+ delete [] nbins;
+ } else {
+ for (Long64_t i = from; i < from + howmany; ++i) {
+ if (!PrintBin(i, bin, options))
+ ++howmany;
+ }
+ }
+ delete [] bin;
+}
+
+//______________________________________________________________________________
+void THnBase::Print(Option_t* options) const
+{
+ // Print a THnBase. If "option" contains:
+ // 'a': print axis details
+ // 'm': print memory usage
+ // 's': print statistics
+ // 'c': print its content, too (this can generate a LOT of output!)
+ // Other options are forwarded to PrintEntries().
+
+ Bool_t optAxis = options && (strchr(options, 'A') || (strchr(options, 'a')));
+ Bool_t optMem = options && (strchr(options, 'M') || (strchr(options, 'm')));
+ Bool_t optStat = options && (strchr(options, 'S') || (strchr(options, 's')));
+ Bool_t optContent = options && (strchr(options, 'C') || (strchr(options, 'c')));
+
+ Printf("%s (*0x%lx): \"%s\" \"%s\"", IsA()->GetName(), (unsigned long)this, GetName(), GetTitle());
+ Printf(" %d dimensions, %g entries in %lld filled bins", GetNdimensions(), GetEntries(), GetNbins());
+
+ if (optAxis) {
+ for (Int_t dim = 0; dim < fNdimensions; ++dim) {
+ TAxis* axis = GetAxis(dim);
+ Printf(" axis %d \"%s\": %d bins (%g..%g), %s bin sizes", dim,
+ axis->GetTitle(), axis->GetNbins(), axis->GetXmin(), axis->GetXmax(),
+ (axis->GetXbins() ? "variable" : "fixed"));
+ }
+ }
+
+ if (optStat) {
+ Printf(" %s error calculation", (GetCalculateErrors() ? "with" : "without"));
+ if (GetCalculateErrors()) {
+ Printf(" Sum(w)=%g, Sum(w^2)=%g", GetSumw(), GetSumw2());
+ for (Int_t dim = 0; dim < fNdimensions; ++dim) {
+ Printf(" axis %d: Sum(w*x)=%g, Sum(w*x^2)=%g", dim, GetSumwx(dim), GetSumwx2(dim));
+ }
+ }
+ }
+
+ if (optMem && InheritsFrom(THnSparse::Class())) {
+ const THnSparse* hsparse = dynamic_cast<const THnSparse*>(this);
+ Printf(" coordinates stored in %d chunks of %d entries\n %g of bins filled using %g of memory compared to an array",
+ hsparse->GetNChunks(), hsparse->GetChunkSize(),
+ hsparse->GetSparseFractionBins(), hsparse->GetSparseFractionMem());
+ }
+
+ if (optContent) {
+ Printf(" BIN CONTENT:");
+ PrintEntries(0, -1, options);
+ }
+}
+
+
+//______________________________________________________________________________
+void THnBase::Browse(TBrowser *b)
+{
+ // Browse a THnSparse: create an entry (ROOT::THnSparseBrowsable) for each
+ // dimension.
+ if (fBrowsables.IsEmpty()) {
+ for (Int_t dim = 0; dim < fNdimensions; ++dim) {
+ fBrowsables[dim] = new ROOT::THnBaseBrowsable(this, dim);
+ }
+ fBrowsables.SetOwner();
+ }
+
+ for (Int_t dim = 0; dim < fNdimensions; ++dim) {
+ b->Add(fBrowsables[dim]);
+ }
+}
+
+
+
+//______________________________________________________________________________
+//
+//
+// Iterator over THnBase bins; internal implementation.
+//
+//______________________________________________________________________________
+ROOT::THnBaseBinIter::~THnBaseBinIter() {
+ // Destruct a bin iterator.
+
+ // Not much to do, but pin vtable
+}
+
+
+//______________________________________________________________________________
+//
+//
+// Iterator over THnBase bins
+//
+//______________________________________________________________________________
+ClassImp(THnIter);
+
+THnIter::~THnIter() {
+ // Destruct a bin iterator.
+ delete fIter;
+}
+
+
+
+
+//______________________________________________________________________________
+//
+// TBrowser helper for THnBase.
+//
+//______________________________________________________________________________
+ClassImp(ROOT::THnBaseBrowsable);
+
+//______________________________________________________________________________
+ROOT::THnBaseBrowsable::THnBaseBrowsable(THnBase* hist, Int_t axis):
+ TNamed(TString::Format("axis %d", axis),
+ TString::Format("Projection on axis %d of %s", axis, hist->IsA()->GetName())),
+ fHist(hist), fAxis(axis), fProj(0)
+{
+ // Construct a THnBaseBrowsable.
+}
+
+//______________________________________________________________________________
+ROOT::THnBaseBrowsable::~THnBaseBrowsable()
+{
+ // Destruct a THnBaseBrowsable.
+ delete fProj;
+}
+
+//______________________________________________________________________________
+void ROOT::THnBaseBrowsable::Browse(TBrowser* b)
+{
+ // Browse an axis of a THnBase, i.e. draw its projection.
+ if (!fProj) {
+ fProj = fHist->Projection(fAxis);
+ }
+ fProj->Draw(b ? b->GetDrawOption() : "");
+ gPad->Update();
+}
+
diff --git a/hist/hist/src/THnSparse.cxx b/hist/hist/src/THnSparse.cxx
index a5a2ee9766ac1e3a9fefd5c8d0294ed6077c77cb..a1cee5284f972b887b82ff7bf58beca937bfa48c 100644
--- a/hist/hist/src/THnSparse.cxx
+++ b/hist/hist/src/THnSparse.cxx
@@ -2,7 +2,7 @@
// Author: Axel Naumann (2007-09-11)
/*************************************************************************
- * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers. *
+ * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
@@ -11,28 +11,86 @@
#include "THnSparse.h"
-#include "TArrayI.h"
#include "TAxis.h"
-#include "TBrowser.h"
#include "TClass.h"
-#include "TCollection.h"
#include "TDataMember.h"
#include "TDataType.h"
-#include "TH1D.h"
-#include "TH2D.h"
-#include "TH3D.h"
-#include "TF1.h"
-#include "TInterpreter.h"
-#include "TMath.h"
-#include "TRandom.h"
-#include "TVirtualPad.h"
-
-#include "TError.h"
-
-#include "HFitInterface.h"
-#include "Fit/SparseData.h"
-#include "Math/MinimizerOptions.h"
-#include "Math/WrappedMultiTF1.h"
+
+namespace {
+//______________________________________________________________________________
+//
+// THnSparseBinIter iterates over all filled bins of a THnSparse.
+//______________________________________________________________________________
+
+ class THnSparseBinIter: public ROOT::THnBaseBinIter {
+ public:
+ THnSparseBinIter(Bool_t respectAxisRange, const THnSparse* hist):
+ ROOT::THnBaseBinIter(respectAxisRange), fHist(hist),
+ fNbins(hist->GetNbins()), fIndex(-1) {
+ // Construct a THnSparseBinIter
+ fCoord = new Int_t[hist->GetNdimensions()];
+ fCoord[0] = -1;
+ }
+ virtual ~THnSparseBinIter() { delete [] fCoord; }
+
+ virtual Int_t GetCoord(Int_t dim) const;
+ virtual Long64_t Next(Int_t* coord = 0);
+
+ private:
+ const THnSparse* fHist;
+ Int_t* fCoord; // coord buffer for fIndex; fCoord[0] == -1 if not yet calculated
+ Long64_t fNbins; // number of bins to iterate over
+ Long64_t fIndex; // current bin index
+ };
+}
+
+Int_t THnSparseBinIter::GetCoord(Int_t dim) const
+{
+ if (fCoord[0] == -1) {
+ fHist->GetBinContent(fIndex, fCoord);
+ }
+ return fCoord[dim];
+}
+
+Long64_t THnSparseBinIter::Next(Int_t* coord /*= 0*/)
+{
+ // Get next bin index (in range if RespectsAxisRange()).
+ // If coords != 0, set it to the index's axis coordinates
+ // (i.e. coord must point to an array of Int_t[fNdimension]
+ if (!fHist) return -1;
+
+ fCoord[0] = -1;
+ Int_t* useCoordBuf = fCoord;
+ if (coord) {
+ useCoordBuf = coord;
+ coord[0] = -1;
+ }
+
+ do {
+ ++fIndex;
+ if (fIndex >= fHist->GetNbins()) {
+ fHist = 0;
+ return -1;
+ }
+ if (RespectsAxisRange()) {
+ fHist->GetBinContent(fIndex, useCoordBuf);
+ }
+ } while (RespectsAxisRange()
+ && !fHist->IsInRange(useCoordBuf)
+ && (fHaveSkippedBin = kTRUE /* assignment! */));
+
+ if (coord && coord[0] == -1) {
+ if (fCoord[0] == -1) {
+ fHist->GetBinContent(fIndex, coord);
+ } else {
+ memcpy(coord, fCoord, fHist->GetNdimensions() * sizeof(Int_t));
+ }
+ }
+
+ return fIndex;
+}
+
+
//______________________________________________________________________________
//
@@ -493,8 +551,7 @@ ClassImp(THnSparse);
//______________________________________________________________________________
THnSparse::THnSparse():
- fNdimensions(0), fChunkSize(1024), fFilledBins(0), fEntries(0),
- fTsumw(0), fTsumw2(-1.), fCompactCoord(0), fIntegral(0), fIntegralStatus(kNoInt)
+ fChunkSize(1024), fFilledBins(0), fCompactCoord(0)
{
// Construct an empty THnSparse.
fBinContent.SetOwner();
@@ -504,9 +561,8 @@ THnSparse::THnSparse():
THnSparse::THnSparse(const char* name, const char* title, Int_t dim,
const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
Int_t chunksize):
- TNamed(name, title), fNdimensions(dim), fChunkSize(chunksize), fFilledBins(0),
- fAxes(dim), fEntries(0), fTsumw(0), fTsumw2(-1.), fTsumwx(dim), fTsumwx2(dim),
- fCompactCoord(0), fIntegral(0), fIntegralStatus(kNoInt)
+ THnBase(name, title, dim, nbins, xmin, xmax),
+ fChunkSize(chunksize), fFilledBins(0), fCompactCoord(0)
{
// Construct a THnSparse with "dim" dimensions,
// with chunksize as the size of the chunks.
@@ -515,13 +571,6 @@ THnSparse::THnSparse(const char* name, const char* title, Int_t dim,
// The arrays "xmin" and "xmax" can be NULL; in that case SetBinEdges()
// must be called for each dimension.
- for (Int_t i = 0; i < fNdimensions; ++i) {
- TAxis* axis = new TAxis(nbins[i], xmin ? xmin[i] : 0., xmax ? xmax[i] : 1.);
- fAxes.AddAtAndExpand(axis, i);
- }
- SetTitle(title);
- fAxes.SetOwner();
-
fCompactCoord = new THnSparseCompactBinCoord(dim, nbins);
fBinContent.SetOwner();
}
@@ -531,17 +580,6 @@ THnSparse::~THnSparse() {
// Destruct a THnSparse
delete fCompactCoord;
- if (fIntegralStatus != kNoInt) delete [] fIntegral;
-}
-
-//______________________________________________________________________________
-void THnSparse::AddBinContent(const Int_t* coord, Double_t v)
-{
- // Add "v" to the content of bin with coordinates "coord"
-
- GetCompactCoord()->SetCoord(coord);
- Long64_t bin = GetBinIndexForCurrentBin(kTRUE);
- return AddBinContent(bin, v);
}
//______________________________________________________________________________
@@ -567,9 +605,8 @@ THnSparseArrayChunk* THnSparse::AddChunk()
}
//______________________________________________________________________________
-THnSparse* THnSparse::CloneEmpty(const char* name, const char* title,
- const TObjArray* axes, Int_t chunksize,
- Bool_t keepTargetAxis) const
+THnBase* THnSparse::CloneEmpty(const char* name, const char* title,
+ const TObjArray* axes, Bool_t keepTargetAxis) const
{
// Create a new THnSparse object that is of the same type as *this,
// but with dimensions and bins given by axes.
@@ -579,7 +616,7 @@ THnSparse* THnSparse::CloneEmpty(const char* name, const char* title,
THnSparse* ret = (THnSparse*)IsA()->New();
ret->SetNameTitle(name, title);
ret->fNdimensions = axes->GetEntriesFast();
- ret->fChunkSize = chunksize;
+ ret->fChunkSize = fChunkSize;
TIter iAxis(axes);
const TAxis* axis = 0;
@@ -612,130 +649,6 @@ THnSparse* THnSparse::CloneEmpty(const char* name, const char* title,
return ret;
}
-//______________________________________________________________________________
-TH1* THnSparse::CreateHist(const char* name, const char* title,
- const TObjArray* axes,
- Bool_t keepTargetAxis ) const {
- // Create an empty histogram with name and title with a given
- // set of axes. Create a TH1D/TH2D/TH3D, depending on the number
- // of elements in axes.
-
- const int ndim = axes->GetSize();
-
- TH1* hist = 0;
- // create hist with dummy axes, we fix them later.
- if (ndim == 1)
- hist = new TH1D(name, title, 1, 0., 1.);
- else if (ndim == 2)
- hist = new TH2D(name, title, 1, 0., 1., 1, 0., 1.);
- else if (ndim == 3)
- hist = new TH3D(name, title, 1, 0., 1., 1, 0., 1., 1, 0., 1.);
- else {
- Error("CreateHist", "Cannot create histogram %s with %d dimensions!", name, ndim);
- return 0;
- }
-
- TAxis* hax[3] = {hist->GetXaxis(), hist->GetYaxis(), hist->GetZaxis()};
- for (Int_t d = 0; d < ndim; ++d) {
- TAxis* reqaxis = (TAxis*)(*axes)[d];
- hax[d]->SetTitle(reqaxis->GetTitle());
- if (!keepTargetAxis && reqaxis->TestBit(TAxis::kAxisRange)) {
- Int_t binFirst = reqaxis->GetFirst();
- if (binFirst == 0) binFirst = 1;
- Int_t binLast = reqaxis->GetLast();
- Int_t nBins = binLast - binFirst + 1;
- if (reqaxis->GetXbins()->GetSize()) {
- // non-uniform bins:
- hax[d]->Set(nBins, reqaxis->GetXbins()->GetArray() + binFirst - 1);
- } else {
- // uniform bins:
- hax[d]->Set(nBins, reqaxis->GetBinLowEdge(binFirst), reqaxis->GetBinUpEdge(binLast));
- }
- } else {
- if (reqaxis->GetXbins()->GetSize()) {
- // non-uniform bins:
- hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXbins()->GetArray());
- } else {
- // uniform bins:
- hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXmin(), reqaxis->GetXmax());
- }
- }
- }
-
- hist->Rebuild();
-
- return hist;
-}
-
-//______________________________________________________________________________
-THnSparse* THnSparse::CreateSparse(const char* name, const char* title,
- const TH1* h, Int_t chunkSize)
-{
- // Create a THnSparse object from a histogram deriving from TH1.
-
- // Get the dimension of the TH1
- int ndim = h->GetDimension();
-
- // Axis properties
- int nbins[3] = {0,0,0};
- double minRange[3] = {0.,0.,0.};
- double maxRange[3] = {0.,0.,0.};
- TAxis* axis[3] = { h->GetXaxis(), h->GetYaxis(), h->GetZaxis() };
- for (int i = 0; i < ndim; ++i) {
- nbins[i] = axis[i]->GetNbins();
- minRange[i] = axis[i]->GetXmin();
- maxRange[i] = axis[i]->GetXmax();
- }
-
- // Create the corresponding THnSparse, depending on the storage
- // type of the TH1. The class name will be "TH??\0" where the first
- // ? is 1,2 or 3 and the second ? indicates the stograge as C, S,
- // I, F or D.
- THnSparse* s = 0;
- const char* cname( h->ClassName() );
- if (cname[0] == 'T' && cname[1] == 'H' && cname[2] >= '1' && cname[2] <= '3' && cname[4] == 0) {
- if (cname[3] == 'F') {
- s = new THnSparseF(name, title, ndim, nbins, minRange, maxRange, chunkSize);
- } else if (cname[3] == 'D') {
- s = new THnSparseD(name, title, ndim, nbins, minRange, maxRange, chunkSize);
- } else if (cname[3] == 'I') {
- s = new THnSparseI(name, title, ndim, nbins, minRange, maxRange, chunkSize);
- } else if (cname[3] == 'S') {
- s = new THnSparseS(name, title, ndim, nbins, minRange, maxRange, chunkSize);
- } else if (cname[3] == 'C') {
- s = new THnSparseC(name, title, ndim, nbins, minRange, maxRange, chunkSize);
- }
- }
- if (!s) {
- ::Warning("THnSparse::CreateSparse", "Unknown Type of Histogram");
- return 0;
- }
-
- for (int i = 0; i < ndim; ++i) {
- s->GetAxis(i)->SetTitle(axis[i]->GetTitle());
- }
-
- // Get the array to know the number of entries of the TH1
- const TArray *array = dynamic_cast<const TArray*>(h);
- if (!array) {
- ::Warning("THnSparse::CreateSparse", "Unknown Type of Histogram");
- return 0;
- }
-
- // Fill the THnSparse with the bins that have content.
- for (int i = 0; i < array->GetSize(); ++i) {
- double value = h->GetBinContent(i);
- double error = h->GetBinError(i);
- if (!value && !error) continue;
- int x[3] = {0,0,0};
- h->GetBinXYZ(i, x[0], x[1], x[2]);
- s->SetBinContent(x, value);
- s->SetBinError(x, error);
- }
-
- return s;
-}
-
//______________________________________________________________________________
void THnSparse::FillExMap()
{
@@ -770,46 +683,14 @@ void THnSparse::FillExMap()
}
//______________________________________________________________________________
-TFitResultPtr THnSparse::Fit(TF1 *f ,Option_t *option ,Option_t *goption)
-{
-// Fit a THnSparse with function f
-//
-// since the data is sparse by default a likelihood fit is performed
-// merging all the regions with empty bins for betetr performance efficiency
-//
-// Since the THnSparse is not drawn no graphics options are passed
-// Here is the list of possible options
-//
-// = "I" Use integral of function in bin instead of value at bin center
-// = "X" Use chi2 method (default is log-likelihood method)
-// = "U" Use a User specified fitting algorithm (via SetFCN)
-// = "Q" Quiet mode (minimum printing)
-// = "V" Verbose mode (default is between Q and V)
-// = "E" Perform better Errors estimation using Minos technique
-// = "B" Use this option when you want to fix one or more parameters
-// and the fitting function is like "gaus", "expo", "poln", "landau".
-// = "M" More. Improve fit results
-// = "R" Use the Range specified in the function range
-
-
- Foption_t fitOption;
-
- if (!TH1::FitOptionsMake(option,fitOption)) return 0;
-
- // The function used to fit cannot be stored in a THnSparse. It
- // cannot be drawn either. Perhaps in the future.
- fitOption.Nostore = true;
- // Use likelihood fit if not specified
- if (!fitOption.Chi2) fitOption.Like = true;
- // create range and minimizer options with default values
- ROOT::Fit::DataRange range(GetNdimensions());
- for ( int i = 0; i < GetNdimensions(); ++i ) {
- TAxis *axis = GetAxis(i);
- range.AddRange(i, axis->GetXmin(), axis->GetXmax());
+void THnSparse::Reserve(Long64_t nbins) {
+ // Initialize storage for nbins
+ if (!fBins.GetSize() && fBinContent.GetSize()) {
+ FillExMap();
+ }
+ if (2 * nbins > fBins.Capacity()) {
+ fBins.Expand(3 * nbins);
}
- ROOT::Math::MinimizerOptions minOption;
-
- return ROOT::Fit::FitObject(this, f , fitOption , minOption, goption, range);
}
//______________________________________________________________________________
@@ -852,16 +733,6 @@ Long64_t THnSparse::GetBin(const Int_t* coord, Bool_t allocate /*= kTRUE*/)
return GetBinIndexForCurrentBin(allocate);
}
-//______________________________________________________________________________
-Double_t THnSparse::GetBinContent(const Int_t *coord) const {
- // Get content of bin with coordinates "coord"
- GetCompactCoord()->SetCoord(coord);
- Long64_t idx = const_cast<THnSparse*>(this)->GetBinIndexForCurrentBin(kFALSE);
- if (idx < 0) return 0.;
- THnSparseArrayChunk* chunk = GetChunk(idx / fChunkSize);
- return chunk->fContent->GetAt(idx % fChunkSize);
-}
-
//______________________________________________________________________________
Double_t THnSparse::GetBinContent(Long64_t idx, Int_t* coord /* = 0 */) const
{
@@ -889,34 +760,15 @@ Double_t THnSparse::GetBinContent(Long64_t idx, Int_t* coord /* = 0 */) const
}
//______________________________________________________________________________
-Double_t THnSparse::GetBinError(const Int_t *coord) const {
- // Get error of bin with coordinates "coord" as
- // BEGIN_LATEX #sqrt{#sum weight^{2}}
+Double_t THnSparse::GetBinError2(Long64_t linidx) const {
+ // Get square of the error of bin addressed by linidx as
+ // BEGIN_LATEX #sum weight^{2}
// END_LATEX
// If errors are not enabled (via Sumw2() or CalculateErrors())
- // return sqrt(contents).
+ // return contents.
if (!GetCalculateErrors())
- return TMath::Sqrt(GetBinContent(coord));
-
- GetCompactCoord()->SetCoord(coord);
- Long64_t idx = const_cast<THnSparse*>(this)->GetBinIndexForCurrentBin(kFALSE);
- if (idx < 0) return 0.;
-
- THnSparseArrayChunk* chunk = GetChunk(idx / fChunkSize);
- return TMath::Sqrt(chunk->fSumw2->GetAt(idx % fChunkSize));
-}
-
-//______________________________________________________________________________
-Double_t THnSparse::GetBinError(Long64_t linidx) const {
- // Get error of bin addressed by linidx as
- // BEGIN_LATEX #sqrt{#sum weight^{2}}
- // END_LATEX
- // If errors are not enabled (via Sumw2() or CalculateErrors())
- // return sqrt(contents).
-
- if (!GetCalculateErrors())
- return TMath::Sqrt(GetBinContent(linidx));
+ return GetBinContent(linidx);
if (linidx < 0) return 0.;
THnSparseArrayChunk* chunk = GetChunk(linidx / fChunkSize);
@@ -924,7 +776,7 @@ Double_t THnSparse::GetBinError(Long64_t linidx) const {
if (!chunk || chunk->fContent->GetSize() < linidx)
return 0.;
- return TMath::Sqrt(chunk->fSumw2->GetAt(linidx));
+ return chunk->fSumw2->GetAt(linidx);
}
@@ -994,44 +846,6 @@ THnSparseCompactBinCoord* THnSparse::GetCompactCoord() const
return fCompactCoord;
}
-//______________________________________________________________________________
-void THnSparse::GetRandom(Double_t *rand, Bool_t subBinRandom /* = kTRUE */)
-{
- // Generate an n-dimensional random tuple based on the histogrammed
- // distribution. If subBinRandom, the returned tuple will be additionally
- // randomly distributed within the randomized bin, using a flat
- // distribution.
-
- // check whether the integral array is valid
- if (fIntegralStatus != kValidInt)
- ComputeIntegral();
-
- // generate a random bin
- Double_t p = gRandom->Rndm();
- Long64_t idx = TMath::BinarySearch(GetNbins() + 1, fIntegral, p);
- const Int_t nStaticBins = 40;
- Int_t bin[nStaticBins];
- Int_t* pBin = bin;
- if (GetNdimensions() > nStaticBins) {
- pBin = new Int_t[GetNdimensions()];
- }
- GetBinContent(idx, pBin);
-
- // convert bin coordinates to real values
- for (Int_t i = 0; i < fNdimensions; i++) {
- rand[i] = GetAxis(i)->GetBinCenter(pBin[i]);
-
- // randomize the vector withing a bin
- if (subBinRandom)
- rand[i] += (gRandom->Rndm() - 0.5) * GetAxis(i)->GetBinWidth(pBin[i]);
- }
- if (pBin != bin) {
- delete [] pBin;
- }
-
- return;
-}
-
//______________________________________________________________________________
Double_t THnSparse::GetSparseFractionBins() const {
// Return the amount of filled bins over all bins
@@ -1074,680 +888,11 @@ Double_t THnSparse::GetSparseFractionMem() const {
}
//______________________________________________________________________________
-Bool_t THnSparse::IsInRange(Int_t *coord) const
-{
- // Check whether bin coord is in range, as defined by TAxis::SetRange().
- // Currently, TAxis::SetRange() does not allow to select all but over- and
- // underflow bins (it instead resets the axis to "no range selected").
- // Instead, simply call
- // TAxis* axis12 = hsparse.GetAxis(12);
- // axis12->SetRange(1, axis12->GetNbins());
- // axis12->SetBit(TAxis::kAxisRange);
- // to deselect the under- and overflow bins in the 12th dimension.
-
- Int_t min = 0;
- Int_t max = 0;
- for (Int_t i = 0; i < fNdimensions; ++i) {
- TAxis *axis = GetAxis(i);
- if (!axis->TestBit(TAxis::kAxisRange)) continue;
- min = axis->GetFirst();
- max = axis->GetLast();
- if (min == 0 && max == 0) {
- // special case where TAxis::SetBit(kAxisRange) and
- // over- and underflow bins are de-selected.
- // first and last are == 0 due to axis12->SetRange(1, axis12->GetNbins());
- min = 1;
- max = axis->GetNbins();
- }
- if (coord[i] < min || coord[i] > max)
- return kFALSE;
- }
- return kTRUE;
-}
-
-//______________________________________________________________________________
-TH1D* THnSparse::Projection(Int_t xDim, Option_t* option /*= ""*/) const
-{
- // Project all bins into a 1-dimensional histogram,
- // keeping only axis "xDim".
- // If "option" contains "E" errors will be calculated.
- // "A" ranges of the taget axes will be ignored.
- // "O" original axis range of the taget axes will be
- // kept, but only bins inside the selected range
- // will be filled.
-
- return (TH1D*) ProjectionAny(1, &xDim, false, option);
-}
-
-//______________________________________________________________________________
-TH2D* THnSparse::Projection(Int_t yDim, Int_t xDim, Option_t* option /*= ""*/) const
-{
- // Project all bins into a 2-dimensional histogram,
- // keeping only axes "xDim" and "yDim".
- //
- // WARNING: just like TH3::Project3D("yx") and TTree::Draw("y:x"),
- // Projection(y,x) uses the first argument to define the y-axis and the
- // second for the x-axis!
- //
- // If "option" contains "E" errors will be calculated.
- // "A" ranges of the taget axes will be ignored.
-
- const Int_t dim[2] = {xDim, yDim};
- return (TH2D*) ProjectionAny(2, dim, false, option);
-}
-
-//______________________________________________________________________________
-TH3D* THnSparse::Projection(Int_t xDim, Int_t yDim, Int_t zDim,
- Option_t* option /*= ""*/) const
-{
- // Project all bins into a 3-dimensional histogram,
- // keeping only axes "xDim", "yDim", and "zDim".
- // If "option" contains "E" errors will be calculated.
- // "A" ranges of the taget axes will be ignored.
- // "O" original axis range of the taget axes will be
- // kept, but only bins inside the selected range
- // will be filled.
-
- const Int_t dim[3] = {xDim, yDim, zDim};
- return (TH3D*) ProjectionAny(3, dim, false, option);
-}
-
-//______________________________________________________________________________
-THnSparse* THnSparse::Projection(Int_t ndim, const Int_t* dim,
- Option_t* option /*= ""*/) const
-{
- // Project all bins into a ndim-dimensional THnSparse histogram,
- // keeping only axes in dim (specifying ndim dimensions)
- // If "option" contains "E" errors will be calculated.
- // "A" ranges of the taget axes will be ignored.
- // "O" original axis range of the taget axes will be
- // kept, but only bins inside the selected range
- // will be filled.
-
- return (THnSparse*) ProjectionAny(ndim, dim, true, option);
-}
-
-
-//______________________________________________________________________________
-TObject* THnSparse::ProjectionAny(Int_t ndim, const Int_t* dim,
- Bool_t wantSparse, Option_t* option /*= ""*/) const
-{
- // Project all bins into a ndim-dimensional THnSparse histogram,
- // keeping only axes in dim (specifying ndim dimensions)
- // If "option" contains "E" errors will be calculated.
- // "A" ranges of the taget axes will be ignored.
- // "O" original axis range of the taget axes will be
- // kept, but only bins inside the selected range
- // will be filled.
-
- TString name(GetName());
- name +="_proj";
-
- for (Int_t d = 0; d < ndim; ++d) {
- name += "_";
- name += dim[d];
- }
-
- TString title(GetTitle());
- Ssiz_t posInsert = title.First(';');
- if (posInsert == kNPOS) {
- title += " projection ";
- for (Int_t d = 0; d < ndim; ++d)
- title += GetAxis(dim[d])->GetTitle();
- } else {
- for (Int_t d = ndim - 1; d >= 0; --d) {
- title.Insert(posInsert, GetAxis(d)->GetTitle());
- if (d)
- title.Insert(posInsert, ", ");
- }
- title.Insert(posInsert, " projection ");
- }
-
- TObjArray newaxes(ndim);
- for (Int_t d = 0; d < ndim; ++d) {
- newaxes.AddAt(GetAxis(dim[d]),d);
- }
-
- THnSparse* sparse = 0;
- TH1* hist = 0;
- TObject* ret = 0;
-
- Bool_t* hadRange = 0;
- Bool_t ignoreTargetRange = (option && (strchr(option, 'A') || strchr(option, 'a')));
- Bool_t keepTargetAxis = ignoreTargetRange || (option && (strchr(option, 'O') || strchr(option, 'o')));
- if (ignoreTargetRange) {
- hadRange = new Bool_t[ndim];
- for (Int_t d = 0; d < ndim; ++d){
- TAxis *axis = GetAxis(dim[d]);
- hadRange[d] = axis->TestBit(TAxis::kAxisRange);
- axis->SetBit(TAxis::kAxisRange, kFALSE);
- }
- }
-
- if (wantSparse)
- ret = sparse = CloneEmpty(name, title, &newaxes, fChunkSize, keepTargetAxis);
- else
- ret = hist = CreateHist(name, title, &newaxes, keepTargetAxis);
-
- if (keepTargetAxis) {
- // make the whole axes visible, i.e. unset the range
- if (wantSparse) {
- for (Int_t d = 0; d < ndim; ++d) {
- sparse->GetAxis(d)->SetRange(0, 0);
- }
- } else {
- hist->GetXaxis()->SetRange(0, 0);
- hist->GetYaxis()->SetRange(0, 0);
- hist->GetZaxis()->SetRange(0, 0);
- }
- }
-
- Bool_t haveErrors = GetCalculateErrors();
- Bool_t wantErrors = (option && (strchr(option, 'E') || strchr(option, 'e'))) || haveErrors;
-
- Int_t* bins = new Int_t[ndim];
- Int_t linbin = 0;
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
-
- Double_t err = 0.;
- Double_t preverr = 0.;
- Double_t v = 0.;
-
- Bool_t haveSkippedBin = false;
-
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- v = GetBinContent(i, coord);
-
- if (!IsInRange(coord)) {
- haveSkippedBin = true;
- continue;
- }
-
- for (Int_t d = 0; d < ndim; ++d) {
- bins[d] = coord[dim[d]];
- if (!keepTargetAxis && GetAxis(dim[d])->TestBit(TAxis::kAxisRange)) {
- bins[d] -= GetAxis(dim[d])->GetFirst() - 1;
- }
- }
-
- if (!wantSparse) {
- if (ndim == 1) linbin = bins[0];
- else if (ndim == 2) linbin = hist->GetBin(bins[0], bins[1]);
- else if (ndim == 3) linbin = hist->GetBin(bins[0], bins[1], bins[2]);
- }
-
- if (wantErrors) {
- if (haveErrors) {
- err = GetBinError(i);
- err *= err;
- } else err = v;
- if (wantSparse) {
- preverr = sparse->GetBinError(bins);
- sparse->SetBinError(bins, TMath::Sqrt(preverr * preverr + err));
- } else {
- preverr = hist->GetBinError(linbin);
- hist->SetBinError(linbin, TMath::Sqrt(preverr * preverr + err));
- }
- }
-
- // only _after_ error calculation, or sqrt(v) is taken into account!
- if (wantSparse)
- sparse->AddBinContent(bins, v);
- else
- hist->AddBinContent(linbin, v);
- }
-
- delete [] bins;
- delete [] coord;
-
- if (wantSparse) {
- // need to check also when producing a THnSparse how to reset the number of entries
- sparse->SetEntries(fEntries);
- } else {
- if (!haveSkippedBin) {
- hist->SetEntries(fEntries);
- } else {
- // re-compute the entries
- // in case of error calculation (i.e. when Sumw2() is set)
- // use the effective entries for the entries
- // since this is the only way to estimate them
- hist->ResetStats();
- Double_t entries = hist->GetEffectiveEntries();
- if (!wantErrors) {
- // to avoid numerical rounding
- entries = TMath::Floor(entries + 0.5);
- }
- hist->SetEntries(entries);
- }
- }
-
- if (hadRange) {
- // reset kAxisRange bit:
- for (Int_t d = 0; d < ndim; ++d)
- GetAxis(dim[d])->SetBit(TAxis::kAxisRange, hadRange[d]);
-
- delete [] hadRange;
- }
-
- return ret;
-}
-
-//______________________________________________________________________________
-void THnSparse::Scale(Double_t c)
-{
- // Scale contents and errors of this histogram by c:
- // this = this * c
- // It does not modify the histogram's number of entries.
-
-
- Double_t nEntries = GetEntries();
- // Scale the contents & errors
- Bool_t haveErrors = GetCalculateErrors();
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- // Get the content of the bin from the current histogram
- Double_t v = GetBinContent(i);
- SetBinContent(i, c * v);
- if (haveErrors) {
- Double_t err = GetBinError(i);
- SetBinError(i, c * err);
- }
- }
- SetEntries(nEntries);
-}
-
-//______________________________________________________________________________
-void THnSparse::AddInternal(const THnSparse* h, Double_t c, Bool_t rebinned)
-{
- // Add() implementation for both rebinned histograms and those with identical
- // binning. See THnSparse::Add().
-
- if (fNdimensions != h->GetNdimensions()) {
- Warning("RebinnedAdd", "Different number of dimensions, cannot carry out operation on the histograms");
- return;
- }
-
- // Trigger error calculation if h has it
- if (!GetCalculateErrors() && h->GetCalculateErrors())
- Sumw2();
- Bool_t haveErrors = GetCalculateErrors();
-
- // Now add the contents: in this case we have the union of the sets of bins
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
- Double_t* x = 0;
- if (rebinned) {
- x = new Double_t[fNdimensions];
- }
-
- if (!fBins.GetSize() && fBinContent.GetSize()) {
- FillExMap();
- }
-
- // Expand the exmap if needed, to reduce collisions
- Long64_t numTargetBins = GetNbins() + h->GetNbins();
- if (2 * numTargetBins > fBins.Capacity()) {
- fBins.Expand(3 * numTargetBins);
- }
-
- // Add to this whatever is found inside the other histogram
- for (Long64_t i = 0; i < h->GetNbins(); ++i) {
- // Get the content of the bin from the second histogram
- Double_t v = h->GetBinContent(i, coord);
-
- Long64_t binidx = -1;
- if (rebinned) {
- // Get the bin center given a coord
- for (Int_t j = 0; j < fNdimensions; ++j)
- x[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
-
- binidx = GetBin(x);
- } else {
- binidx = GetBin(coord);
- }
-
- if (haveErrors) {
- Double_t err1 = GetBinError(binidx);
- Double_t err2 = h->GetBinError(i) * c;
- SetBinError(binidx, TMath::Sqrt(err1 * err1 + err2 * err2));
- }
- // only _after_ error calculation, or sqrt(v) is taken into account!
- AddBinContent(binidx, c * v);
- }
-
-
- delete [] coord;
- delete [] x;
-
- Double_t nEntries = GetEntries() + c * h->GetEntries();
- SetEntries(nEntries);
-}
-
-//______________________________________________________________________________
-void THnSparse::Add(const THnSparse* h, Double_t c)
-{
- // Add contents of h scaled by c to this histogram:
- // this = this + c * h
- // Note that if h has Sumw2 set, Sumw2 is automatically called for this
- // if not already set.
-
- // Check consistency of the input
- if (!CheckConsistency(h, "Add")) return;
-
- AddInternal(h, c, kFALSE);
-}
-
-//______________________________________________________________________________
-void THnSparse::RebinnedAdd(const THnSparse* h, Double_t c)
-{
- // Add contents of h scaled by c to this histogram:
- // this = this + c * h
- // Note that if h has Sumw2 set, Sumw2 is automatically called for this
- // if not already set.
- // In contrast to Add(), RebinnedAdd() does not require consist binning of
- // this and h; instead, each bin's center is used to determine the target bin.
-
- AddInternal(h, c, kTRUE);
-}
-
-
-//______________________________________________________________________________
-Long64_t THnSparse::Merge(TCollection* list)
-{
- // Merge this with a list of THnSparses. All THnSparses provided
- // in the list must have the same bin layout!
-
- if (!list) return 0;
- if (list->IsEmpty()) return (Long64_t)GetEntries();
-
- Long64_t sumNbins = GetNbins();
- TIter iter(list);
- const TObject* addMeObj = 0;
- while ((addMeObj = iter())) {
- const THnSparse* addMe = dynamic_cast<const THnSparse*>(addMeObj);
- if (addMe) {
- sumNbins += addMe->GetNbins();
- }
- }
- if (!fBins.GetSize() && fBinContent.GetSize()) {
- FillExMap();
- }
- if (2 * sumNbins > fBins.Capacity()) {
- fBins.Expand(3 * sumNbins);
- }
-
- iter.Reset();
- while ((addMeObj = iter())) {
- const THnSparse* addMe = dynamic_cast<const THnSparse*>(addMeObj);
- if (!addMe)
- Error("Merge", "Object named %s is not THnSpase! Skipping it.",
- addMeObj->GetName());
- else
- Add(addMe);
- }
- return (Long64_t)GetEntries();
-}
-
-
-//______________________________________________________________________________
-void THnSparse::Multiply(const THnSparse* h)
-{
- // Multiply this histogram by histogram h
- // this = this * h
- // Note that if h has Sumw2 set, Sumw2 is automatically called for this
- // if not already set.
-
- // Check consistency of the input
- if(!CheckConsistency(h, "Multiply"))return;
-
- // Trigger error calculation if h has it
- Bool_t wantErrors = kFALSE;
- if (GetCalculateErrors() || h->GetCalculateErrors())
- wantErrors = kTRUE;
-
- if (wantErrors) Sumw2();
-
- Double_t nEntries = GetEntries();
- // Now multiply the contents: in this case we have the intersection of the sets of bins
- Int_t* coord = new Int_t[fNdimensions];
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- // Get the content of the bin from the current histogram
- Double_t v1 = GetBinContent(i, coord);
- // Now look at the bin with the same coordinates in h
- Double_t v2 = h->GetBinContent(coord);
- SetBinContent(coord, v1 * v2);;
- if (wantErrors) {
- Double_t err1 = GetBinError(coord) * v2;
- Double_t err2 = h->GetBinError(coord) * v1;
- SetBinError(coord,TMath::Sqrt((err2 * err2 + err1 * err1)));
- }
- }
- SetEntries(nEntries);
-
- delete [] coord;
-}
-
-//______________________________________________________________________________
-void THnSparse::Multiply(TF1* f, Double_t c)
-{
- // Performs the operation: this = this*c*f1
- // if errors are defined, 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 THnSparse::Multiply
-
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
- Double_t* points = new Double_t[fNdimensions];
-
-
- Bool_t wantErrors = GetCalculateErrors();
- if (wantErrors) Sumw2();
-
- ULong64_t nEntries = GetNbins();
- for (ULong64_t i = 0; i < nEntries; ++i) {
- Double_t value = GetBinContent(i, coord);
-
- // Get the bin co-ordinates given an coord
- for (Int_t j = 0; j < fNdimensions; ++j)
- points[j] = GetAxis(j)->GetBinCenter(coord[j]);
-
- if (!f->IsInside(points))
- continue;
- TF1::RejectPoint(kFALSE);
-
- // Evaulate function at points
- Double_t fvalue = f->EvalPar(points, NULL);
-
- SetBinContent(coord, c * fvalue * value);
- if (wantErrors) {
- Double_t error(GetBinError(i));
- SetBinError(coord, c * fvalue * error);
- }
- }
-
- delete [] points;
- delete [] coord;
-}
-
-//______________________________________________________________________________
-void THnSparse::Divide(const THnSparse *h)
-{
- // Divide this histogram by h
- // this = this/(h)
- // Note that if h has Sumw2 set, Sumw2 is automatically called for
- // this if not already set.
- // The resulting errors are calculated assuming uncorrelated content.
-
- // Check consistency of the input
- if (!CheckConsistency(h, "Divide"))return;
-
- // Trigger error calculation if h has it
- Bool_t wantErrors=GetCalculateErrors();
- if (!GetCalculateErrors() && h->GetCalculateErrors())
- wantErrors=kTRUE;
-
- // Remember original histogram statistics
- Double_t nEntries = fEntries;
-
- if (wantErrors) Sumw2();
- Bool_t didWarn = kFALSE;
-
- // Now divide the contents: also in this case we have the intersection of the sets of bins
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
- Double_t err = 0.;
- Double_t b22 = 0.;
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- // Get the content of the bin from the first histogram
- Double_t v1 = GetBinContent(i, coord);
- // Now look at the bin with the same coordinates in h
- Double_t v2 = h->GetBinContent(coord);
- if (!v2) {
- v1 = 0.;
- v2 = 1.;
- if (!didWarn) {
- Warning("Divide(h)", "Histogram h has empty bins - division by zero! Setting bin to 0.");
- didWarn = kTRUE;
- }
- }
- SetBinContent(coord, v1 / v2);
- if (wantErrors) {
- Double_t err1 = GetBinError(coord) * v2;
- Double_t err2 = h->GetBinError(coord) * v1;
- b22 = v2 * v2;
- err = (err1 * err1 + err2 * err2) / (b22 * b22);
- SetBinError(coord, TMath::Sqrt(err));
- }
- }
- delete [] coord;
- SetEntries(nEntries);
-}
-
-//______________________________________________________________________________
-void THnSparse::Divide(const THnSparse *h1, const THnSparse *h2, Double_t c1, Double_t c2, Option_t *option)
-{
- // Replace contents of this histogram by multiplication of h1 by h2
- // this = (c1*h1)/(c2*h2)
- // 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 content.
- // 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.
-
-
- TString opt = option;
- opt.ToLower();
- Bool_t binomial = kFALSE;
- if (opt.Contains("b")) binomial = kTRUE;
-
- // Check consistency of the input
- if (!CheckConsistency(h1, "Divide") || !CheckConsistency(h2, "Divide"))return;
- if (!c2) {
- Error("Divide","Coefficient of dividing histogram cannot be zero");
- return;
- }
-
- Reset();
-
- // Trigger error calculation if h1 or h2 have it
- if (!GetCalculateErrors() && (h1->GetCalculateErrors()|| h2->GetCalculateErrors() != 0))
- Sumw2();
-
- // Count filled bins
- Long64_t nFilledBins=0;
-
- // Now divide the contents: we have the intersection of the sets of bins
-
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
- Float_t w = 0.;
- Float_t err = 0.;
- Float_t b22 = 0.;
- Bool_t didWarn = kFALSE;
-
- for (Long64_t i = 0; i < h1->GetNbins(); ++i) {
- // Get the content of the bin from the first histogram
- Double_t v1 = h1->GetBinContent(i, coord);
- // Now look at the bin with the same coordinates in h2
- Double_t v2 = h2->GetBinContent(coord);
- if (!v2) {
- v1 = 0.;
- v2 = 1.;
- if (!didWarn) {
- Warning("Divide(h1, h2)", "Histogram h2 has empty bins - division by zero! Setting bin to 0.");
- didWarn = kTRUE;
- }
- }
- nFilledBins++;
- SetBinContent(coord, c1 * v1 / c2 / v2);
- if(GetCalculateErrors()){
- Double_t err1=h1->GetBinError(coord);
- Double_t err2=h2->GetBinError(coord);
- if (binomial) {
- if (v1 != v2) {
- w = v1 / v2;
- err2 *= w;
- err = TMath::Abs( ( (1. - 2.*w) * err1 * err1 + err2 * err2 ) / (v2 * v2) );
- } else {
- err = 0;
- }
- } else {
- c1 *= c1;
- c2 *= c2;
- b22 = v2 * v2 * c2;
- err1 *= v2;
- err2 *= v1;
- err = c1 * c2 * (err1 * err1 + err2 * err2) / (b22 * b22);
- }
- SetBinError(coord,TMath::Sqrt(err));
- }
- }
-
- delete [] coord;
- fFilledBins = nFilledBins;
-
- // Set as entries in the result histogram the entries in the numerator
- SetEntries(h1->GetEntries());
-}
-
-//______________________________________________________________________________
-Bool_t THnSparse::CheckConsistency(const THnSparse *h, const char *tag) const
+ROOT::THnBaseBinIter* THnSparse::CreateIter(Bool_t respectAxisRange) const
{
- // Consistency check on (some of) the parameters of two histograms (for operations).
-
- if (fNdimensions!=h->GetNdimensions()) {
- Warning(tag,"Different number of dimensions, cannot carry out operation on the histograms");
- return kFALSE;
- }
- for (Int_t dim = 0; dim < fNdimensions; dim++){
- if (GetAxis(dim)->GetNbins()!=h->GetAxis(dim)->GetNbins()) {
- Warning(tag,"Different number of bins on axis %i, cannot carry out operation on the histograms", dim);
- return kFALSE;
- }
- }
- return kTRUE;
-}
-
-//______________________________________________________________________________
-void THnSparse::SetBinEdges(Int_t idim, const Double_t* bins)
-{
- // Set the axis # of bins and bin limits on dimension idim
-
- TAxis* axis = (TAxis*) fAxes[idim];
- axis->Set(axis->GetNbins(), bins);
-}
-
-//______________________________________________________________________________
-void THnSparse::SetBinContent(const Int_t* coord, Double_t v)
-{
- // Set content of bin with coordinates "coord" to "v"
-
- GetCompactCoord()->SetCoord(coord);
- Long_t bin = GetBinIndexForCurrentBin(kTRUE);
- SetBinContent(bin, v);
+ // Create an iterator over all filled bins of a THnSparse.
+ // Use THnIter instead.
+ return new THnSparseBinIter(respectAxisRange, this);
}
//______________________________________________________________________________
@@ -1761,19 +906,26 @@ void THnSparse::SetBinContent(Long64_t bin, Double_t v)
}
//______________________________________________________________________________
-void THnSparse::SetBinError(const Int_t* coord, Double_t e)
+void THnSparse::SetBinError2(Long64_t bin, Double_t e2)
{
- // Set error of bin with coordinates "coord" to "e", enable errors if needed
+ // Set error of bin with index "bin" to "e", enable errors if needed
- GetCompactCoord()->SetCoord(coord);
- Long_t bin = GetBinIndexForCurrentBin(kTRUE);
- SetBinError(bin, e);
+ THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
+ if (!chunk->fSumw2 ) {
+ // if fSumw2 is zero GetCalculateErrors should return false
+ if (GetCalculateErrors()) {
+ Error("SetBinError", "GetCalculateErrors() logic error!");
+ }
+ Sumw2(); // enable error calculation
+ }
+
+ chunk->fSumw2->SetAt(e2, bin % fChunkSize);
}
//______________________________________________________________________________
-void THnSparse::SetBinError(Long64_t bin, Double_t e)
+void THnSparse::AddBinError2(Long64_t bin, Double_t e2)
{
- // Set error of bin with index "bin" to "e", enable errors if needed
+ // Add "e" to error of bin with index "bin", enable errors if needed
THnSparseArrayChunk* chunk = GetChunk(bin / fChunkSize);
if (!chunk->fSumw2 ) {
@@ -1784,7 +936,7 @@ void THnSparse::SetBinError(Long64_t bin, Double_t e)
Sumw2(); // enable error calculation
}
- chunk->fSumw2->SetAt(e * e, bin % fChunkSize);
+ (*chunk->fSumw2)[bin % fChunkSize] += e2;
}
//______________________________________________________________________________
@@ -1802,422 +954,13 @@ void THnSparse::Sumw2()
}
//______________________________________________________________________________
-THnSparse* THnSparse::Rebin(Int_t group) const
-{
- // Combine the content of "group" neighboring bins into
- // a new bin and return the resulting THnSparse.
- // For group=2 and a 3 dimensional histogram, all "blocks"
- // of 2*2*2 bins will be put into a bin.
-
- Int_t* ngroup = new Int_t[GetNdimensions()];
- for (Int_t d = 0; d < GetNdimensions(); ++d)
- ngroup[d] = group;
- THnSparse* ret = Rebin(ngroup);
- delete [] ngroup;
- return ret;
-}
-
-//______________________________________________________________________________
-void THnSparse::SetTitle(const char *title)
-{
- // Change (i.e. set) the title.
- //
- // If title is in the form "stringt;string0;string1;string2 ..."
- // the histogram title is set to stringt, the title of axis0 to string0,
- // of axis1 to string1, of axis2 to string2, etc, just like it is done
- // for TH1/TH2/TH3.
- // To insert the character ";" in one of the titles, one should use "#;"
- // or "#semicolon".
-
- fTitle = title;
- fTitle.ReplaceAll("#;",2,"#semicolon",10);
-
- Int_t endHistTitle = fTitle.First(';');
- if (endHistTitle >= 0) {
- // title contains a ';' so parse the axis titles
- Int_t posTitle = endHistTitle + 1;
- Int_t lenTitle = fTitle.Length();
- Int_t dim = 0;
- while (posTitle > 0 && posTitle < lenTitle && dim < fNdimensions){
- Int_t endTitle = fTitle.Index(";", posTitle);
- TString axisTitle = fTitle(posTitle, endTitle - posTitle);
- axisTitle.ReplaceAll("#semicolon", 10, ";", 1);
- GetAxis(dim)->SetTitle(axisTitle);
- dim++;
- if (endTitle > 0)
- posTitle = endTitle + 1;
- else
- posTitle = -1;
- }
- // Remove axis titles from histogram title
- fTitle.Remove(endHistTitle, lenTitle - endHistTitle);
- }
-
- fTitle.ReplaceAll("#semicolon", 10, ";", 1);
-
-}
-//______________________________________________________________________________
-THnSparse* THnSparse::Rebin(const Int_t* group) const
-{
- // Combine the content of "group" neighboring bins for each dimension
- // into a new bin and return the resulting THnSparse.
- // For group={2,1,1} and a 3 dimensional histogram, pairs of x-bins
- // will be grouped.
-
- Int_t ndim = GetNdimensions();
- TString name(GetName());
- for (Int_t d = 0; d < ndim; ++d)
- name += Form("_%d", group[d]);
-
-
- TString title(GetTitle());
- Ssiz_t posInsert = title.First(';');
- if (posInsert == kNPOS) {
- title += " rebin ";
- for (Int_t d = 0; d < ndim; ++d)
- title += Form("{%d}", group[d]);
- } else {
- for (Int_t d = ndim - 1; d >= 0; --d)
- title.Insert(posInsert, Form("{%d}", group[d]));
- title.Insert(posInsert, " rebin ");
- }
-
- TObjArray newaxes(ndim);
- newaxes.SetOwner();
- for (Int_t d = 0; d < ndim; ++d) {
- newaxes.AddAt(GetAxis(d)->Clone(),d);
- if (group[d] > 1) {
- TAxis* newaxis = (TAxis*) newaxes.At(d);
- Int_t newbins = (newaxis->GetNbins() + group[d] - 1) / group[d];
- if (newaxis->GetXbins() && newaxis->GetXbins()->GetSize()) {
- // variable bins
- Double_t *edges = new Double_t[newbins + 1];
- for (Int_t i = 0; i < newbins + 1; ++i)
- if (group[d] * i <= newaxis->GetNbins())
- edges[i] = newaxis->GetXbins()->At(group[d] * i);
- else edges[i] = newaxis->GetXmax();
- newaxis->Set(newbins, edges);
- delete [] edges;
- } else {
- newaxis->Set(newbins, newaxis->GetXmin(), newaxis->GetXmax());
- }
- }
- }
-
- THnSparse* h = CloneEmpty(name.Data(), title.Data(), &newaxes, fChunkSize, kTRUE);
- Bool_t haveErrors = GetCalculateErrors();
- Bool_t wantErrors = haveErrors;
-
- Int_t* bins = new Int_t[ndim];
- Int_t* coord = new Int_t[fNdimensions];
- memset(coord, 0, sizeof(Int_t) * fNdimensions);
- Double_t err = 0.;
- Double_t preverr = 0.;
- Double_t v = 0.;
-
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- v = GetBinContent(i, coord);
- for (Int_t d = 0; d < ndim; ++d)
- bins[d] = TMath::CeilNint( (double) coord[d]/group[d] );
-
- if (wantErrors) {
- if (haveErrors) {
- err = GetBinError(i);
- err *= err;
- } else err = v;
- preverr = h->GetBinError(bins);
- h->SetBinError(bins, TMath::Sqrt(preverr * preverr + err));
- }
-
- // only _after_ error calculation, or sqrt(v) is taken into account!
- h->AddBinContent(bins, v);
- }
-
- delete [] bins;
- delete [] coord;
-
- h->SetEntries(fEntries);
-
- return h;
-
-}
-
-//______________________________________________________________________________
-void THnSparse::Reset(Option_t * /*option = ""*/)
+void THnSparse::Reset(Option_t *option /*= ""*/)
{
// Clear the histogram
fFilledBins = 0;
- fEntries = 0.;
- fTsumw = 0.;
- fTsumw2 = -1.;
fBins.Delete();
fBinsContinued.Clear();
fBinContent.Delete();
- if (fIntegralStatus != kNoInt) {
- delete [] fIntegral;
- fIntegralStatus = kNoInt;
- }
-}
-
-//______________________________________________________________________________
-Double_t THnSparse::ComputeIntegral()
-{
- // Calculate the integral of the histogram
-
- // delete old integral
- if (fIntegralStatus != kNoInt) {
- delete [] fIntegral;
- fIntegralStatus = kNoInt;
- }
-
- // check number of bins
- if (GetNbins() == 0) {
- Error("ComputeIntegral", "The histogram must have at least one bin.");
- return 0.;
- }
-
- // allocate integral array
- fIntegral = new Double_t [GetNbins() + 1];
- fIntegral[0] = 0.;
-
- // fill integral array with contents of regular bins (non over/underflow)
- Int_t* coord = new Int_t[fNdimensions];
- for (Long64_t i = 0; i < GetNbins(); ++i) {
- Double_t v = GetBinContent(i, coord);
-
- // check whether the bin is regular
- bool regularBin = true;
- for (Int_t dim = 0; dim < fNdimensions; dim++)
- if (coord[dim] < 1 || coord[dim] > GetAxis(dim)->GetNbins()) {
- regularBin = false;
- break;
- }
-
- // if outlayer, count it with zero weight
- if (!regularBin) v = 0.;
-
- fIntegral[i + 1] = fIntegral[i] + v;
- }
- delete [] coord;
-
- // check sum of weights
- if (fIntegral[GetNbins()] == 0.) {
- Error("ComputeIntegral", "No hits in regular bins (non over/underflow).");
- delete [] fIntegral;
- return 0.;
- }
-
- // normalize the integral array
- for (Long64_t i = 0; i <= GetNbins(); ++i)
- fIntegral[i] = fIntegral[i] / fIntegral[GetNbins()];
-
- // set status to valid
- fIntegralStatus = kValidInt;
- return fIntegral[GetNbins()];
-}
-
-//______________________________________________________________________________
-void THnSparse::PrintBin(Long64_t idx, Option_t* options) const
-{
- // Print bin with linex index "idx".
- // For valid options see PrintBin(Long64_t idx, Int_t* bin, Option_t* options).
- Int_t* coord = new Int_t[fNdimensions];
- PrintBin(idx, coord, options);
- delete [] coord;
-}
-
-//______________________________________________________________________________
-Bool_t THnSparse::PrintBin(Long64_t idx, Int_t* bin, Option_t* options) const
-{
- // Print one bin. If "idx" is != -1 use that to determine the bin,
- // otherwise (if "idx" == -1) use the coordinate in "bin".
- // If "options" contains:
- // '0': only print bins with an error or content != 0
- // Return whether the bin was printed (depends on options)
-
- Double_t v = -42;
- if (idx == -1) {
- idx = const_cast<THnSparse*>(this)->GetBin(bin, kFALSE);
- v = GetBinContent(idx);
- } else {
- v = GetBinContent(idx, bin);
- }
-
- if (options && strchr(options, '0') && v == 0.) {
- if (GetCalculateErrors()) {
- if (GetBinError(idx) == 0.) {
- // suppress zeros, and we have one.
- return kFALSE;
- }
- } else {
- // suppress zeros, and we have one.
- return kFALSE;
- }
- }
-
- TString coord;
- for (Int_t dim = 0; dim < fNdimensions; ++dim) {
- coord += bin[dim];
- coord += ',';
- }
- coord.Remove(coord.Length() - 1);
-
- if (GetCalculateErrors()) {
- Double_t err = 0.;
- if (idx != -1) {
- err = GetBinError(idx);
- }
- Printf("Bin at (%s) = %g (+/- %g)", coord.Data(), v, err);
- } else {
- Printf("Bin at (%s) = %g", coord.Data(), v);
- }
-
- return kTRUE;
-}
-
-//______________________________________________________________________________
-void THnSparse::PrintEntries(Long64_t from /*=0*/, Long64_t howmany /*=-1*/,
- Option_t* options /*=0*/) const
-{
- // Print "howmany" entries starting at "from". If "howmany" is -1, print all.
- // If "options" contains:
- // 'x': print in the order of axis bins, i.e. (0,0,...,0), (0,0,...,1),...
- // '0': only print bins with content != 0
-
- if (from < 0) from = 0;
- if (howmany == -1) howmany = GetNbins();
-
- Int_t* bin = new Int_t[fNdimensions];
-
- if (options && (strchr(options, 'x') || strchr(options, 'X'))) {
- Int_t* nbins = new Int_t[fNdimensions];
- for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
- nbins[dim] = GetAxis(dim)->GetNbins();
- bin[dim] = from % nbins[dim];
- from /= nbins[dim];
- }
-
- for (Long64_t i = 0; i < howmany; ++i) {
- if (!PrintBin(-1, bin, options))
- ++howmany;
- // Advance to next bin:
- ++bin[fNdimensions - 1];
- for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
- if (bin[dim] >= nbins[dim]) {
- bin[dim] = 0;
- if (dim > 0) {
- ++bin[dim - 1];
- } else {
- howmany = -1; // aka "global break"
- }
- }
- }
- }
- delete [] nbins;
- } else {
- for (Long64_t i = from; i < from + howmany; ++i) {
- if (!PrintBin(i, bin, options))
- ++howmany;
- }
- }
- delete [] bin;
-}
-
-//______________________________________________________________________________
-void THnSparse::Print(Option_t* options) const
-{
- // Print a THnSparse. If "option" contains:
- // 'a': print axis details
- // 'm': print memory usage
- // 's': print statistics
- // 'c': print its content, too (this can generate a LOT of output!)
- // Other options are forwarded to PrintEntries().
-
- Bool_t optAxis = options && (strchr(options, 'A') || (strchr(options, 'a')));
- Bool_t optMem = options && (strchr(options, 'M') || (strchr(options, 'm')));
- Bool_t optStat = options && (strchr(options, 'S') || (strchr(options, 's')));
- Bool_t optContent = options && (strchr(options, 'C') || (strchr(options, 'c')));
-
- Printf("%s (*0x%lx): \"%s\" \"%s\"", IsA()->GetName(), (unsigned long)this, GetName(), GetTitle());
- Printf(" %d dimensions, %g entries in %lld filled bins", GetNdimensions(), GetEntries(), GetNbins());
-
- if (optAxis) {
- for (Int_t dim = 0; dim < fNdimensions; ++dim) {
- TAxis* axis = GetAxis(dim);
- Printf(" axis %d \"%s\": %d bins (%g..%g), %s bin sizes", dim,
- axis->GetTitle(), axis->GetNbins(), axis->GetXmin(), axis->GetXmax(),
- (axis->GetXbins() ? "variable" : "fixed"));
- }
- }
-
- if (optStat) {
- Printf(" %s error calculation", (GetCalculateErrors() ? "with" : "without"));
- if (GetCalculateErrors()) {
- Printf(" Sum(w)=%g, Sum(w^2)=%g", GetSumw(), GetSumw2());
- for (Int_t dim = 0; dim < fNdimensions; ++dim) {
- Printf(" axis %d: Sum(w*x)=%g, Sum(w*x^2)=%g", dim, GetSumwx(dim), GetSumwx2(dim));
- }
- }
- }
-
- if (optMem) {
- Printf(" coordinates stored in %d chunks of %d entries\n %g of bins filled using %g of memory compared to an array",
- GetNChunks(), GetChunkSize(), GetSparseFractionBins(), GetSparseFractionMem());
- }
-
- if (optContent) {
- Printf(" BIN CONTENT:");
- PrintEntries(0, -1, options);
- }
-}
-
-
-//______________________________________________________________________________
-void THnSparse::Browse(TBrowser *b)
-{
- // Browse a THnSparse: create an entry (ROOT::THnSparseBrowsable) for each
- // dimension.
- if (fBrowsables.IsEmpty()) {
- for (Int_t dim = 0; dim < fNdimensions; ++dim) {
- fBrowsables[dim] = new ROOT::THnSparseBrowsable(this, dim);
- }
- fBrowsables.SetOwner();
- }
-
- for (Int_t dim = 0; dim < fNdimensions; ++dim) {
- b->Add(fBrowsables[dim]);
- }
-}
-
-
-//______________________________________________________________________________
-//______________________________________________________________________________
-ClassImp(ROOT::THnSparseBrowsable);
-
-//______________________________________________________________________________
-ROOT::THnSparseBrowsable::THnSparseBrowsable(THnSparse* hist, Int_t axis):
- TNamed(TString::Format("axis %d", axis),
- TString::Format("Projection on axis %d of sparse histogram", axis)),
- fHist(hist), fAxis(axis), fProj(0)
-{
- // Construct a THnSparseBrowsable.
-}
-
-//______________________________________________________________________________
-ROOT::THnSparseBrowsable::~THnSparseBrowsable()
-{
- // Destruct a THnSparseBrowsable.
- delete fProj;
-}
-
-//______________________________________________________________________________
-void ROOT::THnSparseBrowsable::Browse(TBrowser* b)
-{
- // Browse an axis of a THnSparse, i.e. draw its projection.
- if (!fProj) {
- fProj = fHist->Projection(fAxis);
- }
- fProj->Draw(b ? b->GetDrawOption() : "");
- gPad->Update();
+ ResetBase(option);
}
diff --git a/test/stressHistogram.cxx b/test/stressHistogram.cxx
index 4855f2ffa5d2110db89981d598f86574acd8877c..63b6c88a32cc0276d070cb4dabbeddd0b8d0d593 100644
--- a/test/stressHistogram.cxx
+++ b/test/stressHistogram.cxx
@@ -45,8 +45,8 @@
// Test 12: Interpolation tests for Histograms...............................OK //
// Test 13: Scale tests for Profiles.........................................OK //
// Test 14: Integral tests for Histograms....................................OK //
-// Test 15: TH1-THnSparse Conversion tests...................................OK //
-// Test 16: Filldata tests for Histograms and Sparses........................OK //
+// Test 15: TH1-THn[Sparse] Conversion tests.................................OK //
+// Test 16: Filldata tests for Histograms and THn[Sparse]....................OK //
// Test 17: Reference File Read for Histograms and Profiles..................OK //
// **************************************************************************** //
// stressHistogram: Real Time = 64.01 seconds Cpu Time = 63.89 seconds //
@@ -62,6 +62,7 @@
#include "TH2.h"
#include "TH3.h"
#include "TH2.h"
+#include "THn.h"
#include "THnSparse.h"
#include "TProfile.h"
@@ -145,8 +146,8 @@ void FillProfiles(TProfile* p1, TProfile* p2, Double_t c1 = 1.0, Double_t c2 = 1
int equals(const char* msg, TH1D* h1, TH1D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
int equals(const char* msg, TH2D* h1, TH2D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
int equals(const char* msg, TH3D* h1, TH3D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
-int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
-int equals(const char* msg, THnSparse* h1, TH1* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
+int equals(const char* msg, THnBase* h1, THnBase* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
+int equals(const char* msg, THnBase* h1, TH1* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
int equals(Double_t n1, Double_t n2, double ERRORLIMIT = defaultErrorLimit);
int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT = defaultErrorLimit);
ostream& operator<<(ostream& out, TH1D* h);
@@ -896,9 +897,10 @@ bool testAdd3DProfile2()
return ret;
}
-bool testAddSparse()
+template<typename HIST>
+bool testAddHn()
{
- // Tests the Add method for Sparse Histograms
+ // Tests the Add method for n-dimensional Histograms
Double_t c = r.Rndm();
@@ -908,9 +910,9 @@ bool testAddSparse()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("tS-s1", "s1", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("tS-s2", "s2", 3, bsize, xmin, xmax);
- THnSparseD* s3 = new THnSparseD("tS-s3", "s3=s1+c*s2", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("tS-s1", "s1", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("tS-s2", "s2", 3, bsize, xmin, xmax);
+ HIST* s3 = new HIST("tS-s3", "s3=s1+c*s2", 3, bsize, xmin, xmax);
s1->Sumw2();s2->Sumw2();s3->Sumw2();
@@ -933,7 +935,7 @@ bool testAddSparse()
}
s1->Add(s2, c);
- bool ret = equals("AddSparse", s3, s1, cmpOptStats , 1E-10);
+ bool ret = equals(TString::Format("AddHn<%s>", HIST::Class()->GetName()), s3, s1, cmpOptStats , 1E-10);
delete s2;
delete s3;
return ret;
@@ -1390,7 +1392,8 @@ bool testMul3D2()
return ret;
}
-bool testMulSparse()
+template<typename HIST>
+bool testMulHn()
{
// Tests the Multiply method for Sparse Histograms
@@ -1400,9 +1403,9 @@ bool testMulSparse()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("m3D2-s1", "s1-Title", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("m3D2-s2", "s2-Title", 3, bsize, xmin, xmax);
- THnSparseD* s3 = new THnSparseD("m3D2-s3", "s3=s1*s2", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("m3D2-s1", "s1-Title", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("m3D2-s2", "s2-Title", 3, bsize, xmin, xmax);
+ HIST* s3 = new HIST("m3D2-s3", "s3=s1*s2", 3, bsize, xmin, xmax);
s1->Sumw2();s2->Sumw2();s3->Sumw2();
@@ -1453,7 +1456,7 @@ bool testMulSparse()
s1->Multiply(s2);
- bool ret = equals("MultSparse", s3, s1, cmpOptNone, 1E-10);
+ bool ret = equals(TString::Format("MultHn<%s>", HIST::Class()->GetName()), s3, s1, cmpOptNone, 1E-10);
delete s2;
delete s3;
return ret;
@@ -1670,6 +1673,7 @@ bool testMulF3D2()
return status;
}
+template <typename HIST>
bool testMulFND()
{
const UInt_t nDims = 3;
@@ -1681,8 +1685,8 @@ bool testMulFND()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("mfND-s2", "s2=f*s2", nDims, bsize, xmin, xmax);
+ HIST* s1 = new HIST("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
+ HIST* s2 = new HIST("mfND-s2", "s2=f*s2", nDims, bsize, xmin, xmax);
TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);
@@ -1701,12 +1705,13 @@ bool testMulFND()
s1->Multiply(f, c1);
- int status = equals("MULF HND", s1, s2);
+ int status = equals(TString::Format("MULF HND<%s>", HIST::Class()->GetName()), s1, s2);
delete s1;
delete f;
return status;
}
+template<typename HIST>
bool testMulFND2()
{
const UInt_t nDims = 3;
@@ -1718,8 +1723,8 @@ bool testMulFND2()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("mfND-s2", "s2=f*s2", nDims, bsize, xmin, xmax);
+ HIST* s1 = new HIST("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
+ HIST* s2 = new HIST("mfND-s2", "s2=f*s2", nDims, bsize, xmin, xmax);
TF2* f = new TF2("sin2", "sin(x)*cos(y)",
minRange - 2, maxRange + 2,
@@ -1742,7 +1747,7 @@ bool testMulFND2()
s1->Multiply(f, c1);
- int status = equals("MULF HND2", s1, s2);
+ int status = equals(TString::Format("MULF HND2<%s>", HIST::Class()->GetName()), s1, s2);
delete s1;
delete f;
return status;
@@ -2217,7 +2222,8 @@ bool testDivide3D2()
return ret;
}
-bool testDivSparse1()
+template <typename HIST>
+bool testDivHn1()
{
// Tests the first Divide method for 3D Histograms
@@ -2231,9 +2237,9 @@ bool testDivSparse1()
const Double_t c1 = 1;
const Double_t c2 = 1;
- THnSparseD* s1 = new THnSparseD("dND1-s1", "s1-Title", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("dND1-s2", "s2-Title", 3, bsize, xmin, xmax);
- THnSparseD* s4 = new THnSparseD("dND1-s4", "s4=s3*s2)", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("dND1-s1", "s1-Title", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("dND1-s2", "s2-Title", 3, bsize, xmin, xmax);
+ HIST* s4 = new HIST("dND1-s4", "s4=s3*s2)", 3, bsize, xmin, xmax);
s1->Sumw2();s2->Sumw2();s4->Sumw2();
@@ -2254,7 +2260,7 @@ bool testDivSparse1()
s4->Fill(points, 1.0);
}
- THnSparseD* s3 = new THnSparseD("dND1-s3", "s3=(c1*s1)/(c2*s2)", 3, bsize, xmin, xmax);
+ HIST* s3 = new HIST("dND1-s3", "s3=(c1*s1)/(c2*s2)", 3, bsize, xmin, xmax);
s3->Divide(s1, s2, c1, c2);
s4->Multiply(s3);
@@ -2271,15 +2277,15 @@ bool testDivSparse1()
s4->SetBinError(coord, sqrt(error));
}
- bool ret = equals("DivideND1", s1, s4, cmpOptNone, 1E-6);
+ bool ret = equals(TString::Format("DivideND1<%s>", HIST::Class()->GetName()), s1, s4, cmpOptNone, 1E-6);
delete s1;
delete s2;
delete s3;
return ret;
}
-
-bool testDivSparse2()
+template <typename HIST>
+bool testDivHn2()
{
// Tests the second Divide method for 3D Histograms
@@ -2293,9 +2299,9 @@ bool testDivSparse2()
const Double_t c1 = 1;
const Double_t c2 = 1;
- THnSparseD* s1 = new THnSparseD("dND2-s1", "s1-Title", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("dND2-s2", "s2-Title", 3, bsize, xmin, xmax);
- THnSparseD* s4 = new THnSparseD("dND2-s4", "s4=s3*s2)", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("dND2-s1", "s1-Title", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("dND2-s2", "s2-Title", 3, bsize, xmin, xmax);
+ HIST* s4 = new HIST("dND2-s4", "s4=s3*s2)", 3, bsize, xmin, xmax);
s1->Sumw2();s2->Sumw2();s4->Sumw2();
@@ -2316,10 +2322,10 @@ bool testDivSparse2()
s4->Fill(points, 1.0);
}
- THnSparseD* s3 = static_cast<THnSparseD*>( s1->Clone() );
+ HIST* s3 = static_cast<HIST*>( s1->Clone() );
s3->Divide(s2);
- THnSparseD* s5 = new THnSparseD("dND2-s5", "s5=(c1*s1)/(c2*s2)", 3, bsize, xmin, xmax);
+ HIST* s5 = new HIST("dND2-s5", "s5=(c1*s1)/(c2*s2)", 3, bsize, xmin, xmax);
s5->Divide(s1,s2);
s4->Multiply(s3);
@@ -2336,7 +2342,7 @@ bool testDivSparse2()
s4->SetBinError(coord, sqrt(error));
}
- bool ret = equals("DivideND2", s1, s4, cmpOptNone, 1E-6);
+ bool ret = equals(TString::Format("DivideND2<%s>", HIST::Class()->GetName()), s1, s4, cmpOptNone, 1E-6);
delete s1;
delete s2;
@@ -2897,7 +2903,8 @@ bool testCloneProfile3D()
return ret;
}
-bool testCloneSparse()
+template <typename HIST>
+bool testCloneHn()
{
// Tests the clone method for Sparse histograms
@@ -2908,7 +2915,7 @@ bool testCloneSparse()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("clS-s1","s1-Title", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("clS-s1","s1-Title", 3, bsize, xmin, xmax);
for ( Int_t i = 0; i < nEvents * nEvents; ++i ) {
Double_t points[3];
@@ -2918,9 +2925,9 @@ bool testCloneSparse()
s1->Fill(points);
}
- THnSparseD* s2 = (THnSparseD*) s1->Clone();
+ HIST* s2 = (HIST*) s1->Clone();
- bool ret = equals("Clone Function THnSparse", s1, s2);
+ bool ret = equals(TString::Format("Clone Function %s", HIST::Class()->GetName()), s1, s2);
delete s1;
return ret;
}
@@ -3147,9 +3154,10 @@ bool testWriteReadProfile3D()
return ret;
}
-bool testWriteReadSparse()
+template <typename HIST>
+bool testWriteReadHn()
{
- // Tests the write and read methods for Sparse Histograms
+ // Tests the write and read methods for n-dim Histograms
Int_t bsize[] = { TMath::Nint( r.Uniform(1, 5) ),
TMath::Nint( r.Uniform(1, 5) ),
@@ -3158,7 +3166,7 @@ bool testWriteReadSparse()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("wrS-s1","s1-Title", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("wrS-s1","s1-Title", 3, bsize, xmin, xmax);
s1->Sumw2();
for ( Int_t i = 0; i < nEvents * nEvents; ++i ) {
@@ -3174,9 +3182,9 @@ bool testWriteReadSparse()
f.Close();
TFile f2("tmpHist.root");
- THnSparseD* s2 = static_cast<THnSparseD*> ( f2.Get("wrS-s1") );
+ HIST* s2 = static_cast<HIST*> ( f2.Get("wrS-s1") );
- bool ret = equals("Read/Write Hist 3D", s1, s2, cmpOptNone);
+ bool ret = equals(TString::Format("Read/Write Hist %s", HIST::Class()->GetName()), s1, s2, cmpOptNone);
delete s1;
return ret;
}
@@ -3525,9 +3533,10 @@ bool testMergeProf3D()
return ret;
}
-bool testMergeSparse()
+template <typename HIST>
+bool testMergeHn()
{
- // Tests the merge method for Sparse Histograms
+ // Tests the merge method for n-dim Histograms
Int_t bsize[] = { TMath::Nint( r.Uniform(1, 5) ),
TMath::Nint( r.Uniform(1, 5) ),
@@ -3536,10 +3545,10 @@ bool testMergeSparse()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("mergeS-s1", "s1-Title", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("mergeS-s2", "s2-Title", 3, bsize, xmin, xmax);
- THnSparseD* s3 = new THnSparseD("mergeS-s3", "s3-Title", 3, bsize, xmin, xmax);
- THnSparseD* s4 = new THnSparseD("mergeS-s4", "s4-Title", 3, bsize, xmin, xmax);
+ HIST* s1 = new HIST("mergeS-s1", "s1-Title", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("mergeS-s2", "s2-Title", 3, bsize, xmin, xmax);
+ HIST* s3 = new HIST("mergeS-s3", "s3-Title", 3, bsize, xmin, xmax);
+ HIST* s4 = new HIST("mergeS-s4", "s4-Title", 3, bsize, xmin, xmax);
s1->Sumw2();s2->Sumw2();s3->Sumw2();
@@ -3576,7 +3585,7 @@ bool testMergeSparse()
s1->Merge(list);
- bool ret = equals("MergeSparse", s1, s4, cmpOptNone, 1E-10);
+ bool ret = equals(TString::Format("MergeHn<%s>", HIST::Class()->GetName()), s1, s4, cmpOptNone, 1E-10);
delete s1;
delete s2;
delete s3;
@@ -6154,6 +6163,12 @@ bool testConversion1D()
THnSparse* s1f = THnSparse::CreateSparse("s1f", "s1fTitle", h1f);
THnSparse* s1d = THnSparse::CreateSparse("s1d", "s1dTitle", h1d);
+ TH1* h1cn = (TH1*) h1c->Clone("h1cn");
+ TH1* h1sn = (TH1*) h1s->Clone("h1sn");
+ TH1* h1in = (TH1*) h1i->Clone("h1in");
+ TH1* h1fn = (TH1*) h1f->Clone("h1fn");
+ TH1* h1dn = (TH1*) h1s->Clone("h1dn");
+
int status = 0;
status += equals("TH1-THnSparseC", s1c, h1c);
status += equals("TH1-THnSparseS", s1s, h1s);
@@ -6167,6 +6182,24 @@ bool testConversion1D()
delete s1f;
delete s1d;
+ THn* n1c = THn::CreateHn("n1c", "n1cTitle", h1cn);
+ THn* n1s = THn::CreateHn("n1s", "n1sTitle", h1sn);
+ THn* n1i = THn::CreateHn("n1i", "n1iTitle", h1in);
+ THn* n1f = THn::CreateHn("n1f", "n1fTitle", h1fn);
+ THn* n1d = THn::CreateHn("n1d", "n1dTitle", h1dn);
+
+ status += equals("TH1-THnC", n1c, h1cn);
+ status += equals("TH1-THnS", n1s, h1sn);
+ status += equals("TH1-THnI", n1i, h1in);
+ status += equals("TH1-THnF", n1f, h1fn);
+ status += equals("TH1-THnD", n1d, h1dn);
+
+ delete n1c;
+ delete n1s;
+ delete n1i;
+ delete n1f;
+ delete n1d;
+
return status;
}
@@ -6213,6 +6246,12 @@ bool testConversion2D()
THnSparse* s2f = THnSparse::CreateSparse("s2f", "s2fTitle", h2f);
THnSparse* s2d = THnSparse::CreateSparse("s2d", "s2dTitle", h2d);
+ TH2* h2cn = (TH2*) h2c->Clone("h2cn");
+ TH2* h2sn = (TH2*) h2s->Clone("h2sn");
+ TH2* h2in = (TH2*) h2i->Clone("h2in");
+ TH2* h2fn = (TH2*) h2f->Clone("h2fn");
+ TH2* h2dn = (TH2*) h2d->Clone("h2dn");
+
int status = 0;
status += equals("TH2-THnSparseC", s2c, h2c);
status += equals("TH2-THnSparseS", s2s, h2s);
@@ -6226,6 +6265,24 @@ bool testConversion2D()
delete s2f;
delete s2d;
+ THn* n2c = THn::CreateHn("n2c", "n2cTitle", h2cn);
+ THn* n2s = THn::CreateHn("n2s", "n2sTitle", h2sn);
+ THn* n2i = THn::CreateHn("n2i", "n2iTitle", h2in);
+ THn* n2f = THn::CreateHn("n2f", "n2fTitle", h2fn);
+ THn* n2d = THn::CreateHn("n2d", "n2dTitle", h2dn);
+
+ status += equals("TH2-THnC", n2c, h2cn);
+ status += equals("TH2-THnS", n2s, h2sn);
+ status += equals("TH2-THnI", n2i, h2in);
+ status += equals("TH2-THnF", n2f, h2fn);
+ status += equals("TH2-THnD", n2d, h2dn);
+
+ delete n2c;
+ delete n2s;
+ delete n2i;
+ delete n2f;
+ delete n2d;
+
return status;
}
@@ -6278,6 +6335,12 @@ bool testConversion3D()
THnSparse* s3f = THnSparse::CreateSparse("s3f", "s3fTitle", h3f);
THnSparse* s3d = THnSparse::CreateSparse("s3d", "s3dTitle", h3d);
+ TH3* h3cn = (TH3*) h3c->Clone("h3cn");
+ TH3* h3sn = (TH3*) h3s->Clone("h3sn");
+ TH3* h3in = (TH3*) h3i->Clone("h3in");
+ TH3* h3fn = (TH3*) h3f->Clone("h3fn");
+ TH3* h3dn = (TH3*) h3d->Clone("h3dn");
+
int status = 0;
status += equals("TH3-THnSparseC", s3c, h3c);
status += equals("TH3-THnSparseS", s3s, h3s);
@@ -6291,6 +6354,24 @@ bool testConversion3D()
delete s3f;
delete s3d;
+ THn* n3c = THn::CreateHn("n3c", "n3cTitle", h3cn);
+ THn* n3s = THn::CreateHn("n3s", "n3sTitle", h3sn);
+ THn* n3i = THn::CreateHn("n3i", "n3iTitle", h3in);
+ THn* n3f = THn::CreateHn("n3f", "n3fTitle", h3fn);
+ THn* n3d = THn::CreateHn("n3d", "n3dTitle", h3dn);
+
+ status += equals("TH3-THnC", n3c, h3cn);
+ status += equals("TH3-THnS", n3s, h3sn);
+ status += equals("TH3-THnI", n3i, h3in);
+ status += equals("TH3-THnF", n3f, h3fn);
+ status += equals("TH3-THnD", n3d, h3dn);
+
+ delete n3c;
+ delete n3s;
+ delete n3i;
+ delete n3f;
+ delete n3d;
+
return status;
}
@@ -7513,9 +7594,10 @@ bool test2DRebinProfile()
return ret;
}
-bool testSparseRebin1()
+template <typename HIST>
+bool testHnRebin1()
{
- // Tests rebin method for Sparse Histogram
+ // Tests rebin method for n-dim Histogram
const int rebin = TMath::Nint( r.Uniform(minRebin, maxRebin) );
@@ -7529,8 +7611,8 @@ bool testSparseRebin1()
Double_t xmin[] = {minRange, minRange, minRange};
Double_t xmax[] = {maxRange, maxRange, maxRange};
- THnSparseD* s1 = new THnSparseD("rebin1-s1","s1-Title", 3, bsize, xmin, xmax);
- THnSparseD* s2 = new THnSparseD("rebin1-s2","s2-Title", 3, bsizeRebin, xmin, xmax);
+ HIST* s1 = new HIST("rebin1-s1","s1-Title", 3, bsize, xmin, xmax);
+ HIST* s2 = new HIST("rebin1-s2","s2-Title", 3, bsizeRebin, xmin, xmax);
for ( Int_t i = 0; i < nEvents; ++i ) {
Double_t points[3];
@@ -7541,9 +7623,9 @@ bool testSparseRebin1()
s2->Fill(points);
}
- THnSparse* s3 = s1->Rebin(rebin);
+ HIST* s3 = (HIST*)s1->Rebin(rebin);
- bool ret = equals("THnSparse Rebin 1", s2, s3);
+ bool ret = equals(TString::Format("%s Rebin 1", HIST::Class()->GetName()), s2, s3);
delete s1;
delete s2;
return ret;
@@ -8246,6 +8328,7 @@ private:
TH1D* hw1ZY;
THnSparseD* s3;
+ THnD* n3;
bool buildWithWeights;
@@ -8345,6 +8428,7 @@ public:
Double_t xmin[] = {lower_limit, lower_limit, lower_limit};
Double_t xmax[] = {upper_limit, upper_limit, upper_limit};
s3 = new THnSparseD("s3","s3", 3, bsize, xmin, xmax);
+ n3 = new THnD("n3","n3", 3, bsize, xmin, xmax);
}
@@ -8396,6 +8480,7 @@ public:
delete hw1ZY;
delete s3;
+ delete n3;
// delete all histogram in gROOT
TList * l = gROOT->GetList();
@@ -8414,7 +8499,7 @@ public:
void buildHistograms()
{
- if (h3->GetSumw2N() ) s3->Sumw2();
+ if (h3->GetSumw2N() ) {s3->Sumw2(); n3->Sumw2();}
for ( int ix = 0; ix <= h3->GetXaxis()->GetNbins() + 1; ++ix ) {
double xc = h3->GetXaxis()->GetBinCenter(ix);
@@ -8431,6 +8516,7 @@ public:
Double_t points[] = {x,y,z};
s3->Fill(points);
+ n3->Fill(points);
h2XY->Fill(x,y);
h2XZ->Fill(x,z);
@@ -8494,6 +8580,7 @@ public:
void buildHistogramsWithWeights()
{
s3->Sumw2();
+ n3->Sumw2();
for ( int ix = 0; ix <= h3->GetXaxis()->GetNbins() + 1; ++ix ) {
double xc = h3->GetXaxis()->GetBinCenter(ix);
@@ -8511,6 +8598,7 @@ public:
Double_t points[] = {x,y,z};
s3->Fill(points,w);
+ n3->Fill(points,w);
h2XY->Fill(x,y,w);
h2XZ->Fill(x,z,w);
@@ -8586,6 +8674,7 @@ public:
Double_t points[] = {x,y,z};
s3->Fill(points);
+ n3->Fill(points);
if ( h3->GetXaxis()->FindBin(x) >= xmin && h3->GetXaxis()->FindBin(x) <= xmax &&
h3->GetYaxis()->FindBin(y) >= ymin && h3->GetYaxis()->FindBin(y) <= ymax &&
@@ -8664,12 +8753,16 @@ public:
h1Y->GetXaxis()->SetRange(ymin, ymax);
h1Z->GetXaxis()->SetRange(zmin, zmax);
- // Neet to set up the rest of the ranges!
+ // Need to set up the rest of the ranges!
s3->GetAxis(1)->SetRange(xmin, xmax);
s3->GetAxis(2)->SetRange(ymin, ymax);
s3->GetAxis(3)->SetRange(zmin, zmax);
+ n3->GetAxis(1)->SetRange(xmin, xmax);
+ n3->GetAxis(2)->SetRange(ymin, ymax);
+ n3->GetAxis(3)->SetRange(zmin, zmax);
+
buildWithWeights = false;
}
@@ -8912,6 +9005,13 @@ public:
status += equals("STH3 -> YZ", h2YZ, (TH2D*) s3->Projection(2,1), options);
status += equals("STH3 -> ZX", h2ZX, (TH2D*) s3->Projection(0,2), options);
status += equals("STH3 -> ZY", h2ZY, (TH2D*) s3->Projection(1,2), options);
+
+ status += equals("THn3 -> XY", h2XY, (TH2D*) n3->Projection(1,0), options);
+ status += equals("THn3 -> XZ", h2XZ, (TH2D*) n3->Projection(2,0), options);
+ status += equals("THn3 -> YX", h2YX, (TH2D*) n3->Projection(0,1), options);
+ status += equals("THn3 -> YZ", h2YZ, (TH2D*) n3->Projection(2,1), options);
+ status += equals("THn3 -> ZX", h2ZX, (TH2D*) n3->Projection(0,2), options);
+ status += equals("THn3 -> ZY", h2ZY, (TH2D*) n3->Projection(1,2), options);
options = 0;
if ( defaultEqualOptions & cmpOptPrint )
cout << "----------------------------------------------" << endl;
@@ -8921,6 +9021,10 @@ public:
status += equals("STH3 -> X", h1X, (TH1D*) s3->Projection(0), options);
status += equals("STH3 -> Y", h1Y, (TH1D*) s3->Projection(1), options);
status += equals("STH3 -> Z", h1Z, (TH1D*) s3->Projection(2), options);
+
+ status += equals("THn3 -> X", h1X, (TH1D*) n3->Projection(0), options);
+ status += equals("THn3 -> Y", h1Y, (TH1D*) n3->Projection(1), options);
+ status += equals("THn3 -> Z", h1Z, (TH1D*) n3->Projection(2), options);
options = 0;
if ( defaultEqualOptions & cmpOptPrint )
cout << "----------------------------------------------" << endl;
@@ -9007,19 +9111,19 @@ int stressHistogram()
rangeTestPointer };
// Test 4
- const unsigned int numberOfRebin = 10;
+ const unsigned int numberOfRebin = 11;
pointer2Test rebinTestPointer[numberOfRebin] = { testIntegerRebin, testIntegerRebinProfile,
testIntegerRebinNoName, testIntegerRebinNoNameProfile,
testArrayRebin, testArrayRebinProfile,
test2DRebin, test3DRebin, test2DRebinProfile,
- testSparseRebin1};
+ testHnRebin1<THnD>, testHnRebin1<THnSparseD>};
struct TTestSuite rebinTestSuite = { numberOfRebin,
"Histogram Rebinning..............................................",
rebinTestPointer };
// Test 5
// Add Tests
- const unsigned int numberOfAdds = 21;
+ const unsigned int numberOfAdds = 22;
pointer2Test addTestPointer[numberOfAdds] = { testAdd1, testAddProfile1,
testAdd2, testAddProfile2,
testAdd3,
@@ -9032,7 +9136,8 @@ int stressHistogram()
testAdd2D2, testAdd2DProfile2,
testAdd3D1, testAdd3DProfile1,
testAdd3D2, testAdd3DProfile2,
- testAddSparse
+ testAddHn<THnSparseD>,
+ testAddHn<THnD>
};
struct TTestSuite addTestSuite = { numberOfAdds,
"Add tests for 1D, 2D and 3D Histograms and Profiles..............",
@@ -9040,16 +9145,20 @@ int stressHistogram()
// Test 6
// Multiply Tests
- const unsigned int numberOfMultiply = 17;
+ const unsigned int numberOfMultiply = 20;
pointer2Test multiplyTestPointer[numberOfMultiply] = { testMul1, testMul2,
testMulVar1, testMulVar2,
testMul2D1, testMul2D2,
testMul3D1, testMul3D2,
- testMulSparse,
+ testMulHn<THnD>,
+ testMulHn<THnSparseD>,
testMulF1D, testMulF1D2,
testMulF2D, testMulF2D2,
testMulF3D, testMulF3D2,
- testMulFND, testMulFND2
+ testMulFND<THnD>,
+ testMulFND<THnSparseD>,
+ testMulFND2<THnD>,
+ testMulFND2<THnSparseD>
};
struct TTestSuite multiplyTestSuite = { numberOfMultiply,
"Multiply tests for 1D, 2D and 3D Histograms......................",
@@ -9057,12 +9166,15 @@ int stressHistogram()
// Test 7
// Divide Tests
- const unsigned int numberOfDivide = 10;
+ const unsigned int numberOfDivide = 12;
pointer2Test divideTestPointer[numberOfDivide] = { testDivide1, testDivide2,
testDivideVar1, testDivideVar2,
testDivide2D1, testDivide2D2,
testDivide3D1, testDivide3D2,
- testDivSparse1, testDivSparse2
+ testDivHn1<THnD>,
+ testDivHn1<THnSparseD>,
+ testDivHn2<THnD>,
+ testDivHn2<THnSparseD>
};
struct TTestSuite divideTestSuite = { numberOfDivide,
"Divide tests for 1D, 2D and 3D Histograms........................",
@@ -9075,7 +9187,7 @@ int stressHistogram()
// Test 8
// Copy Tests
- const unsigned int numberOfCopy = 25;
+ const unsigned int numberOfCopy = 26;
pointer2Test copyTestPointer[numberOfCopy] = { testAssign1D, testAssignProfile1D,
testAssignVar1D, testAssignProfileVar1D,
testCopyConstructor1D, testCopyConstructorProfile1D,
@@ -9088,7 +9200,7 @@ int stressHistogram()
testAssign3D, testAssignProfile3D,
testCopyConstructor3D, testCopyConstructorProfile3D,
testClone3D, testCloneProfile3D,
- testCloneSparse
+ testCloneHn<THnD>, testCloneHn<THnSparseD>
};
struct TTestSuite copyTestSuite = { numberOfCopy,
"Copy tests for 1D, 2D and 3D Histograms and Profiles.............",
@@ -9096,12 +9208,12 @@ int stressHistogram()
// Test 9
// WriteRead Tests
- const unsigned int numberOfReadwrite = 9;
- pointer2Test readwriteTestPointer[numberOfReadwrite] = { testWriteRead1D, testWriteReadProfile1D,
- testWriteReadVar1D, testWriteReadProfileVar1D,
- testWriteRead2D, testWriteReadProfile2D,
- testWriteRead3D, testWriteReadProfile3D,
- testWriteReadSparse
+ const unsigned int numberOfReadwrite = 10;
+ pointer2Test readwriteTestPointer[numberOfReadwrite] = { testWriteRead1D, testWriteReadProfile1D,
+ testWriteReadVar1D, testWriteReadProfileVar1D,
+ testWriteRead2D, testWriteReadProfile2D,
+ testWriteRead3D, testWriteReadProfile3D,
+ testWriteReadHn<THnD>,testWriteReadHn<THnSparseD>
};
struct TTestSuite readwriteTestSuite = { numberOfReadwrite,
"Read/Write tests for 1D, 2D and 3D Histograms and Profiles.......",
@@ -9109,12 +9221,12 @@ int stressHistogram()
// Test 10
// Merge Tests
- const unsigned int numberOfMerge = 42;
+ const unsigned int numberOfMerge = 43;
pointer2Test mergeTestPointer[numberOfMerge] = { testMerge1D, testMergeProf1D,
testMergeVar1D, testMergeProfVar1D,
testMerge2D, testMergeProf2D,
testMerge3D, testMergeProf3D,
- testMergeSparse,
+ testMergeHn<THnD>, testMergeHn<THnSparseD>,
testMerge1DLabelSame, testMergeProf1DLabelSame,
testMerge2DLabelSame, testMergeProf2DLabelSame,
testMerge3DLabelSame, testMergeProf3DLabelSame,
@@ -9182,14 +9294,14 @@ int stressHistogram()
integralTestPointer };
// Test 15
- // TH1-THnSparse Conversions Tests
+ // TH1-THn[Sparse] Conversions Tests
const unsigned int numberOfConversions = 3;
pointer2Test conversionsTestPointer[numberOfConversions] = { testConversion1D,
testConversion2D,
testConversion3D,
};
struct TTestSuite conversionsTestSuite = { numberOfConversions,
- "TH1-THnSparse Conversion tests...................................",
+ "TH1-THn[Sparse] Conversion tests.................................",
conversionsTestPointer };
// Test 16
@@ -9346,7 +9458,7 @@ void FillProfiles(TProfile* p1, TProfile* p2, Double_t c1, Double_t c2)
// Methods for histogram comparisions
-int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options, double ERRORLIMIT)
+int equals(const char* msg, THnBase* h1, THnBase* h2, int options, double ERRORLIMIT)
{
options = options | defaultEqualOptions;
bool print = options & cmpOptPrint;
@@ -9387,7 +9499,7 @@ int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options, double ER
}
// Statistical tests:
- // No statistical tests possible for THnSparse so far...
+ // No statistical tests possible for THnBase so far...
// if ( compareStats )
// differents += compareStatistics( h1, h2, debug, ERRORLIMIT);
@@ -9398,7 +9510,7 @@ int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options, double ER
return differents;
}
-int equals(const char* msg, THnSparse* s, TH1* h2, int options, double ERRORLIMIT)
+int equals(const char* msg, THnBase* s, TH1* h2, int options, double ERRORLIMIT)
{
options = options | defaultEqualOptions;
bool print = options & cmpOptPrint;
@@ -9419,7 +9531,7 @@ int equals(const char* msg, THnSparse* s, TH1* h2, int options, double ERRORLIMI
TArray* array = dynamic_cast<TArray*>(h2);
if ( !array )
- Fatal( "equals(const char* msg, THnSparse* s, TH1* h2, int options, double ERRORLIMIT)" ,"NO ARRAY!");
+ Fatal( "equals(const char* msg, THnBase* s, TH1* h2, int options, double ERRORLIMIT)" ,"NO ARRAY!");
Int_t* coord = new Int_t[3];
for (Long64_t i = 0; i < s->GetNbins(); ++i)