Skip to content

Commit 91d46f2

Browse files
ferdymercuryferdymercury
authored
[hist] Allow GetRandom to work with density-histograms
Fixes #9530
1 parent 05d4e38 commit 91d46f2

File tree

8 files changed

+34
-19
lines changed

8 files changed

+34
-19
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: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -2505,12 +2505,15 @@ 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 option (optional) Set it to "width" if your non-uniform bin contents
2509+
/// represent a density rather than counts. In that case, we internally multiply
2510+
/// by the bin width/area so that fIntegral is in the "counts" and not on the original "density" domain.
25082511
/// \return 1 if success, 0 if integral is zero, NAN if onlyPositive-test fails
25092512

2510-
Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
2513+
Double_t TH1::ComputeIntegral(Bool_t onlyPositive, Option_t *option)
25112514
{
25122515
if (fBuffer) BufferEmpty();
2513-
2516+
bool useArea = TString(option).ToLower().Contains("width");
25142517
// delete previously computed integral (if any)
25152518
if (fIntegral) delete [] fIntegral;
25162519

@@ -2523,11 +2526,20 @@ Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
25232526
fIntegral = new Double_t[nbins + 2];
25242527
Int_t ibin = 0; fIntegral[ibin] = 0;
25252528

2529+
2530+
2531+
y=y*xWidth*yWidth*zWidth;
25262532
for (Int_t binz=1; binz <= nbinsz; ++binz) {
2533+
Double_t zWidth = (fDimension > 2) ? fZaxis.GetBinWidth(binz) : 1;
25272534
for (Int_t biny=1; biny <= nbinsy; ++biny) {
2535+
Double_t yWidth= (fDimension > 1) ? fYaxis.GetBinWidth(biny) : 1;
25282536
for (Int_t binx=1; binx <= nbinsx; ++binx) {
2537+
Double_t xWidth = fXaxis.GetBinWidth(binx);
25292538
++ibin;
25302539
Double_t y = RetrieveBinContent(GetBin(binx, biny, binz));
2540+
if (useArea)
2541+
y *= xWidth * yWidth * zWidth;
2542+
25312543
if (onlyPositive && y < 0) {
25322544
Error("ComputeIntegral","Bin content is negative - return a NaN value");
25332545
fIntegral[nbins] = TMath::QuietNaN();
@@ -5005,13 +5017,14 @@ void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) cons
50055017
/// is evaluated, normalized to one.
50065018
///
50075019
/// @param rng (optional) Random number generator pointer used (default is gRandom)
5020+
/// @param option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than counts
50085021
///
50095022
/// The integral is automatically recomputed if the number of entries
50105023
/// is not the same then when the integral was computed.
5011-
/// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.
5024+
/// @note Only valid for 1-d histograms. Use GetRandom2 or GetRandom3 otherwise.
50125025
/// If the histogram has a bin with negative content a NaN is returned
50135026

5014-
Double_t TH1::GetRandom(TRandom * rng) const
5027+
Double_t TH1::GetRandom(TRandom * rng, Option_t *option) const
50155028
{
50165029
if (fDimension > 1) {
50175030
Error("GetRandom","Function only valid for 1-d histograms");
@@ -5021,10 +5034,10 @@ Double_t TH1::GetRandom(TRandom * rng) const
50215034
Double_t integral = 0;
50225035
// compute integral checking that all bins have positive content (see ROOT-5894)
50235036
if (fIntegral) {
5024-
if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true);
5037+
if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true, option);
50255038
else integral = fIntegral[nbinsx];
50265039
} else {
5027-
integral = ((TH1*)this)->ComputeIntegral(true);
5040+
integral = ((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: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1156,19 +1156,20 @@ Double_t TH2::GetCovariance(Int_t axis1, Int_t axis2) const
11561156
/// @param[out] x reference to random generated x value
11571157
/// @param[out] y reference to random generated y value
11581158
/// @param[in] rng (optional) Random number generator pointer used (default is gRandom)
1159+
/// @param[in] option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than 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) integral = ComputeIntegral(true, option);
11691170
else integral = fIntegral[nbins];
11701171
} else {
1171-
integral = ComputeIntegral(true);
1172+
integral = ComputeIntegral(true, option);
11721173
}
11731174
if (integral == 0 ) { x = 0; y = 0; return;}
11741175
// 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: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1265,14 +1265,15 @@ Double_t TH3::GetCovariance(Int_t axis1, Int_t axis2) const
12651265

12661266

12671267
////////////////////////////////////////////////////////////////////////////////
1268-
/// Return 3 random numbers along axis x , y and z distributed according
1268+
/// Return 3 random numbers along axis x, y and z distributed according
12691269
/// to the cell-contents of this 3-dim histogram
12701270
/// @param[out] x reference to random generated x value
12711271
/// @param[out] y reference to random generated y value
12721272
/// @param[out] z reference to random generated z value
12731273
/// @param[in] rng (optional) Random number generator pointer used (default is gRandom)
1274+
/// @param[in] option (optional) Set it to "width" if your non-uniform bin contents represent a density rather than 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,10 @@ 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) integral = ComputeIntegral(true, option);
12861287
else integral = fIntegral[nbins];
12871288
} else {
1288-
integral = ComputeIntegral(true);
1289+
integral = ComputeIntegral(true, option);
12891290
}
12901291
if (integral == 0 ) { x = 0; y = 0; z = 0; return;}
12911292
// case histogram has negative bins

0 commit comments

Comments
 (0)