Hist1D division

[aminnj]

Division of two Hist1Ds currently crashes. It works with a small modification

--- a/src/arithmatics.jl
+++ b/src/arithmatics.jl
@@ -36,8 +36,8 @@ for T in (:Hist1D,)
         _f(counts) = any(x -> x<0, counts)
         (_f(h1.hist.weights) || _f(h2.hist.weights)) && error("Can't do / when bin count is negative")
         _hist = /(h1.hist, h2.hist)
-        _sumw2 =  h1.sumw2 / (h2.hist.weights .^ 2) +
-                (sqrt.(h2.sumw2) * h1.hist.weights / (h2.hist.weights .^ 2)) .^ 2
+        _sumw2 =  @. h1.sumw2 / (h2.hist.weights ^ 2) +
+                (sqrt(h2.sumw2) * h1.hist.weights / (h2.hist.weights ^ 2)) ^ 2

         ($T)(_hist, _sumw2)
     end

Though the output _sumw2 is a bunch of zeros if h1 and h2 are filled without weights. Maybe there needs to be a check inside here which "materializes" each of the sumw2s (gaussian/sqrt by default) before doing the error propagation?

[aminnj]

OK, matches ROOT.

julia> using StatsBase
julia> h1 = Hist1D([0.5,1.5,1.5,2.5], Weights(ones(4)), 0:3);
julia> h2 = Hist1D([1.5,2.5,2.5], Weights(ones(3)), 0:3);
julia> (h1/h2).sumw2 |> println
[NaN, 6.0, 0.375]

import ROOT

h1 = ROOT.TH1F("h1", "", 3, 0, 3); h1.Sumw2()
h2 = ROOT.TH1F("h2", "", 3, 0, 3); h2.Sumw2()

for x in [0.5,1.5,1.5,2.5]: h1.Fill(x)
for x in [1.5,2.5,2.5]: h2.Fill(x)

h1.Divide(h2)

print([round(h1.GetBinContent(i), 4) for i in range(1,h1.GetNbinsX()+1)])
# [0.0, 2.0, 0.5]

print([round(h1.GetBinError(i)**2, 4) for i in range(1,h1.GetNbinsX()+1)])
# [0.0, 6.0, 0.375]

[aminnj]

Ignore/close this and let's just move the discussion to #8 .