Skip to content

Commit 97b81ce

Browse files
[hist] Allow GetRandom to work with density-histograms
To compute the CDF of histograms, one has to take into account if the original histogram contains event counts or a density. Support for the latter case is added in this commit. Fixes #9530 Decided not to clang-format the headers, because this would break the existing indentation. Co-authored-by: Stephan Hageboeck <stephan.hageboeck@cern.ch>
1 parent 7f4454a commit 97b81ce

File tree

8 files changed

+39
-22
lines changed

8 files changed

+39
-22
lines changed

hist/hist/inc/TH1.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -444,7 +444,7 @@ class TH1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
444444
virtual Double_t Chisquare(TF1 * f1, Option_t *option = "") const;
445445
static Int_t CheckConsistency(const TH1* h1, const TH1* h2);
446446
virtual void ClearUnderflowAndOverflow();
447-
virtual Double_t ComputeIntegral(Bool_t onlyPositive = false);
447+
virtual Double_t ComputeIntegral(Bool_t onlyPositive = false, Option_t *option = "");
448448
TObject* Clone(const char *newname = "") const override;
449449
void Copy(TObject &hnew) const override;
450450
virtual void DirectoryAutoAdd(TDirectory *);
@@ -550,7 +550,7 @@ class TH1 : public TNamed, public TAttLine, public TAttFill, public TAttMarker {
550550
TVirtualHistPainter *GetPainter(Option_t *option="");
551551

552552
virtual Int_t GetQuantiles(Int_t n, Double_t *xp, const Double_t *p = nullptr);
553-
virtual Double_t GetRandom(TRandom * rng = nullptr) const;
553+
virtual Double_t GetRandom(TRandom *rng = nullptr, Option_t *option = "") const;
554554
virtual void GetStats(Double_t *stats) const;
555555
virtual Double_t GetStdDev(Int_t axis=1) const;
556556
virtual Double_t GetStdDevError(Int_t axis=1) const;

hist/hist/inc/TH2.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,7 +102,7 @@ class TH2 : public TH1 {
102102
virtual Double_t GetBinErrorUp(Int_t binx, Int_t biny) { return TH1::GetBinErrorUp( GetBin(binx, biny) ); }
103103
virtual Double_t GetCorrelationFactor(Int_t axis1=1,Int_t axis2=2) const;
104104
virtual Double_t GetCovariance(Int_t axis1=1,Int_t axis2=2) const;
105-
virtual void GetRandom2(Double_t &x, Double_t &y, TRandom * rng = nullptr);
105+
virtual void GetRandom2(Double_t &x, Double_t &y, TRandom * rng = nullptr, Option_t *option = "");
106106
void GetStats(Double_t *stats) const override;
107107
Double_t Integral(Option_t *option="") const override;
108108
//virtual Double_t Integral(Int_t, Int_t, Option_t * ="") const {return 0;}

hist/hist/inc/TH2Poly.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ class TH2Poly : public TH2 {
8888
using TH2::Interpolate;
8989
Bool_t Divide(TF1 *, Double_t) override;
9090
Bool_t Multiply(TF1 *, Double_t) override;
91-
Double_t ComputeIntegral(Bool_t) override;
91+
Double_t ComputeIntegral(Bool_t, Option_t *) override;
9292
TH1 * FFT(TH1*, Option_t * ) override;
9393
virtual TH1 * GetAsymmetry(TH1* , Double_t, Double_t);
9494
virtual Double_t Interpolate(Double_t, Double_t);

hist/hist/inc/TH3.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ class TH3 : public TH1, public TAtt3D {
119119
virtual Double_t GetBinWithContent3(Double_t c, Int_t &binx, Int_t &biny, Int_t &binz, Int_t firstx=0, Int_t lastx=0,Int_t firsty=0, Int_t lasty=0, Int_t firstz=0, Int_t lastz=0, Double_t maxdiff=0) const;
120120
virtual Double_t GetCorrelationFactor(Int_t axis1=1,Int_t axis2=2) const;
121121
virtual Double_t GetCovariance(Int_t axis1=1,Int_t axis2=2) const;
122-
virtual void GetRandom3(Double_t &x, Double_t &y, Double_t &, TRandom * rng = nullptr);
122+
virtual void GetRandom3(Double_t &x, Double_t &y, Double_t &, TRandom * rng = nullptr, Option_t *option = "");
123123
void GetStats(Double_t *stats) const override;
124124
Double_t Integral(Option_t *option="") const override;
125125
virtual Double_t Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Option_t *option="") const;

hist/hist/src/TH1.cxx

Lines changed: 20 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -2505,12 +2505,17 @@ void TH1::ClearUnderflowAndOverflow()
25052505
/// The resulting integral is normalized to 1.
25062506
/// If the routine is called with the onlyPositive flag set an error will
25072507
/// be produced in case of negative bin content and a NaN value returned
2508+
/// \param onlyPositive If set to true, an error will be produced and NaN will be returned
2509+
/// when a bin with negative number of entries is encountered.
2510+
/// \param option
2511+
/// - `""` (default) Compute the cumulative density function assuming current bin contents represent counts.
2512+
/// - `"width"` Computes the cumulative density function assuming current bin contents represent densities.
25082513
/// \return 1 if success, 0 if integral is zero, NAN if onlyPositive-test fails
25092514

2510-
Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
2515+
Double_t TH1::ComputeIntegral(Bool_t onlyPositive, Option_t *option)
25112516
{
25122517
if (fBuffer) BufferEmpty();
2513-
2518+
bool useArea = TString(option).Contains("width", TString::kIgnoreCase);
25142519
// delete previously computed integral (if any)
25152520
if (fIntegral) delete [] fIntegral;
25162521

@@ -2524,10 +2529,16 @@ Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
25242529
Int_t ibin = 0; fIntegral[ibin] = 0;
25252530

25262531
for (Int_t binz=1; binz <= nbinsz; ++binz) {
2532+
Double_t zWidth = (fDimension > 2) ? fZaxis.GetBinWidth(binz) : 1;
25272533
for (Int_t biny=1; biny <= nbinsy; ++biny) {
2534+
Double_t yWidth = (fDimension > 1) ? fYaxis.GetBinWidth(biny) : 1;
25282535
for (Int_t binx=1; binx <= nbinsx; ++binx) {
2536+
Double_t xWidth = fXaxis.GetBinWidth(binx);
25292537
++ibin;
25302538
Double_t y = RetrieveBinContent(GetBin(binx, biny, binz));
2539+
if (useArea)
2540+
y *= xWidth * yWidth * zWidth;
2541+
25312542
if (onlyPositive && y < 0) {
25322543
Error("ComputeIntegral","Bin content is negative - return a NaN value");
25332544
fIntegral[nbins] = TMath::QuietNaN();
@@ -5005,13 +5016,14 @@ void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) cons
50055016
/// is evaluated, normalized to one.
50065017
///
50075018
/// @param rng (optional) Random number generator pointer used (default is gRandom)
5019+
/// @param option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than counts
50085020
///
50095021
/// The integral is automatically recomputed if the number of entries
50105022
/// is not the same then when the integral was computed.
5011-
/// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.
5012-
/// If the histogram has a bin with negative content a NaN is returned
5023+
/// @note Only valid for 1-d histograms. Use GetRandom2 or GetRandom3 otherwise.
5024+
/// If the histogram has a bin with negative content, a NaN is returned.
50135025

5014-
Double_t TH1::GetRandom(TRandom * rng) const
5026+
Double_t TH1::GetRandom(TRandom *rng, Option_t *option) const
50155027
{
50165028
if (fDimension > 1) {
50175029
Error("GetRandom","Function only valid for 1-d histograms");
@@ -5021,10 +5033,11 @@ Double_t TH1::GetRandom(TRandom * rng) const
50215033
Double_t integral = 0;
50225034
// compute integral checking that all bins have positive content (see ROOT-5894)
50235035
if (fIntegral) {
5024-
if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true);
5036+
if (fIntegral[nbinsx + 1] != fEntries)
5037+
integral = const_cast<TH1 *>(this)->ComputeIntegral(true, option);
50255038
else integral = fIntegral[nbinsx];
50265039
} else {
5027-
integral = ((TH1*)this)->ComputeIntegral(true);
5040+
integral = const_cast<TH1 *>(this)->ComputeIntegral(true, option);
50285041
}
50295042
if (integral == 0) return 0;
50305043
// return a NaN in case some bins have negative content

hist/hist/src/TH2.cxx

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1146,7 +1146,6 @@ Double_t TH2::GetCovariance(Int_t axis1, Int_t axis2) const
11461146
return sumwxy/sumw - sumwx/sumw*sumwy/sumw;
11471147
}
11481148

1149-
11501149
////////////////////////////////////////////////////////////////////////////////
11511150
/// Return 2 random numbers along axis x and y distributed according
11521151
/// to the cell-contents of this 2-D histogram.
@@ -1156,19 +1155,22 @@ Double_t TH2::GetCovariance(Int_t axis1, Int_t axis2) const
11561155
/// @param[out] x reference to random generated x value
11571156
/// @param[out] y reference to random generated y value
11581157
/// @param[in] rng (optional) Random number generator pointer used (default is gRandom)
1158+
/// @param[in] option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than
1159+
/// counts
11591160

1160-
void TH2::GetRandom2(Double_t &x, Double_t &y, TRandom * rng)
1161+
void TH2::GetRandom2(Double_t &x, Double_t &y, TRandom *rng, Option_t *option)
11611162
{
11621163
Int_t nbinsx = GetNbinsX();
11631164
Int_t nbinsy = GetNbinsY();
11641165
Int_t nbins = nbinsx*nbinsy;
11651166
Double_t integral;
11661167
// compute integral checking that all bins have positive content (see ROOT-5894)
11671168
if (fIntegral) {
1168-
if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral(true);
1169+
if (fIntegral[nbins + 1] != fEntries)
1170+
integral = ComputeIntegral(true, option);
11691171
else integral = fIntegral[nbins];
11701172
} else {
1171-
integral = ComputeIntegral(true);
1173+
integral = ComputeIntegral(true, option);
11721174
}
11731175
if (integral == 0 ) { x = 0; y = 0; return;}
11741176
// case histogram has negative bins

hist/hist/src/TH2Poly.cxx

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1695,7 +1695,7 @@ Bool_t TH2Poly::Multiply(TF1 *, Double_t)
16951695
}
16961696
////////////////////////////////////////////////////////////////////////////////
16971697
/// NOT IMPLEMENTED for TH2Poly
1698-
Double_t TH2Poly::ComputeIntegral(Bool_t )
1698+
Double_t TH2Poly::ComputeIntegral(Bool_t, Option_t *)
16991699
{
17001700
Error("ComputeIntegral", "Not implemented for TH2Poly");
17011701
return TMath::QuietNaN();

hist/hist/src/TH3.cxx

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1263,16 +1263,17 @@ Double_t TH3::GetCovariance(Int_t axis1, Int_t axis2) const
12631263
return 0;
12641264
}
12651265

1266-
12671266
////////////////////////////////////////////////////////////////////////////////
1268-
/// Return 3 random numbers along axis x , y and z distributed according
1267+
/// Return 3 random numbers along axis x, y and z distributed according
12691268
/// to the cell-contents of this 3-dim histogram
12701269
/// @param[out] x reference to random generated x value
12711270
/// @param[out] y reference to random generated y value
12721271
/// @param[out] z reference to random generated z value
12731272
/// @param[in] rng (optional) Random number generator pointer used (default is gRandom)
1273+
/// @param[in] option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than
1274+
/// counts
12741275

1275-
void TH3::GetRandom3(Double_t &x, Double_t &y, Double_t &z, TRandom * rng)
1276+
void TH3::GetRandom3(Double_t &x, Double_t &y, Double_t &z, TRandom *rng, Option_t *option)
12761277
{
12771278
Int_t nbinsx = GetNbinsX();
12781279
Int_t nbinsy = GetNbinsY();
@@ -1282,10 +1283,11 @@ void TH3::GetRandom3(Double_t &x, Double_t &y, Double_t &z, TRandom * rng)
12821283
Double_t integral;
12831284
// compute integral checking that all bins have positive content (see ROOT-5894)
12841285
if (fIntegral) {
1285-
if (fIntegral[nbins+1] != fEntries) integral = ComputeIntegral(true);
1286+
if (fIntegral[nbins + 1] != fEntries)
1287+
integral = ComputeIntegral(true, option);
12861288
else integral = fIntegral[nbins];
12871289
} else {
1288-
integral = ComputeIntegral(true);
1290+
integral = ComputeIntegral(true, option);
12891291
}
12901292
if (integral == 0 ) { x = 0; y = 0; z = 0; return;}
12911293
// case histogram has negative bins

0 commit comments

Comments
 (0)