|
| 1 | +module WVZReportExt |
| 2 | + |
| 3 | +export significance_table, print_sigtable |
| 4 | + |
| 5 | +include(joinpath(@__DIR__, "../../src/alltags.jl")) |
| 6 | + |
| 7 | +using PrettyTables, Serialization, Measurements, FHist |
| 8 | + |
| 9 | +""" |
| 10 | + rebinscan(S, B; atleast=1, from=:right, by = (s,b) -> s/sqrt(b)) |
| 11 | +
|
| 12 | +Take two `Hist1D`, try to find the best bin-edges for rebinning, given these two conditions: |
| 13 | +
|
| 14 | +- a metric to maximize, `s` and `b` are signal counts and background counts in each new bin |
| 15 | +- `atleast` specifcy the minimal number of `s+b` in each newly formed bin |
| 16 | +
|
| 17 | +`from` keyword argument can be used to specifcy the direction of the scan, defautls to `:right` assuming |
| 18 | +the high-er end of the histogram is signal-like and one of the histograms is monotonic. |
| 19 | +
|
| 20 | +The function returns the values of new bin-edges including both ends (of course, the ends are identical to before). |
| 21 | +
|
| 22 | +The metric can also be replaced by the more accurate `sqrt(2*((s + b) * log(1 + s/b) - s ))`. |
| 23 | +""" |
| 24 | +function rebinscan(S, B; atleast=1, from=:right, by = (s,b) -> s/sqrt(b)) |
| 25 | + _binedges = binedges(S) |
| 26 | + _binedges == binedges(B) || error("Bin edges aren't compatible") |
| 27 | + Scounts = bincounts(S) |
| 28 | + Bcounts = bincounts(B) |
| 29 | + if from == :right |
| 30 | + Scounts = reverse(Scounts) |
| 31 | + BCounts = reverse(Bcounts) |
| 32 | + end |
| 33 | + |
| 34 | + tempS = tempB = 0.0 |
| 35 | + score = -Inf |
| 36 | + newEdges = [1] # first edge at the beginning |
| 37 | + for i = 2:lastindex(_binedges) # one more edge than bins |
| 38 | + s = Scounts[i-1] |
| 39 | + b = Bcounts[i-1] |
| 40 | + tempS += s |
| 41 | + tempB += b |
| 42 | + |
| 43 | + tempS + tempB < atleast && continue # not enough statistics |
| 44 | + newScore = by(tempS, tempB) |
| 45 | + if newScore < score # score starts to drop |
| 46 | + push!(newEdges, i-1) |
| 47 | + tempS = s |
| 48 | + tempB = b |
| 49 | + score = -Inf |
| 50 | + end |
| 51 | + score = newScore |
| 52 | + end |
| 53 | + push!(newEdges, lastindex(_binedges)) |
| 54 | + return _binedges[newEdges] |
| 55 | +end |
| 56 | + |
| 57 | +function _significance(signal, bkg) |
| 58 | + S = integral(signal) |
| 59 | + B = integral(bkg) |
| 60 | + Sig = sqrt(2*((S + B) * log(1 + S/B) - S)) |
| 61 | + dS = sqrt(sum(signal.sumw2)) |
| 62 | + dB = sqrt(sum(bkg.sumw2)) |
| 63 | + dSigdS = log(1 + S/B) / Sig |
| 64 | + dSigdB = (log(1 + S/B) - S/B) / Sig |
| 65 | + err = sqrt((dSigdS * dS)^2 + (dSigdB * dB)^2) |
| 66 | + return Sig ± err |
| 67 | +end |
| 68 | + |
| 69 | + |
| 70 | +function significance_matrix(input_dir::AbstractString) |
| 71 | + Ms = map(ALL_TAGS) do tag |
| 72 | + deserialize(joinpath(input_dir,"$(tag).jlserialize")) |
| 73 | + end |
| 74 | + significance_matrix(Ms) |
| 75 | +end |
| 76 | + |
| 77 | +function significance_matrix(Ms) |
| 78 | + res = mapreduce(vcat, Ms) do M |
| 79 | + N = nbins(M[:DF__BDT__NOMINAL]) |
| 80 | + hists = rebin.([M[:SFinZ__BDT__NOMINAL], M[:SFnoZ__BDT__NOMINAL], M[:DF__BDT__NOMINAL]], N) |
| 81 | + N = nbins(M[:ZZCR__Njet__NOMINAL]) |
| 82 | + push!(hists, rebin(M[:ZZCR__Njet__NOMINAL], N)) |
| 83 | + N = nbins(M[:ttZCR__Njet__NOMINAL]) |
| 84 | + push!(hists, rebin(M[:ttZCR__Njet__NOMINAL], N)) |
| 85 | + |
| 86 | + permutedims(hists) |
| 87 | + end |
| 88 | + res |
| 89 | +end |
| 90 | + |
| 91 | +""" |
| 92 | + significance_table() |
| 93 | + significance_table(significance_matrix(); ) |
| 94 | +
|
| 95 | +# Examples |
| 96 | +
|
| 97 | +```julia |
| 98 | +julia> M = significance_table() |
| 99 | +14×6 Matrix{Any}: |
| 100 | + "Signal" 10.66±0.067 … 1.072±0.017 2.151±0.037 |
| 101 | + "ZZ" 1219.9±5.4 555.7±2.4 69.19±0.74 |
| 102 | + "Zjets" -0.019±0.13 -0.0004±0.0004 0.63±0.39 |
| 103 | + "Zgamma" 0.0±0.0 0.0±0.0 0.0±0.0 |
| 104 | + "ttbar" 0.0±0.0 0.0±0.0 0.43±0.13 |
| 105 | + "WZ" 0.358±0.096 … 0.0044±0.0032 0.295±0.061 |
| 106 | + "tZ" 0.012±0.012 0.0±0.0 0.179±0.043 |
| 107 | + "ttZ" 1.231±0.08 0.0111±0.0074 73.62±0.61 |
| 108 | + "tWZ" 0.56±0.11 0.0095±0.0095 14.71±0.56 |
| 109 | + "VBS" 10.231±0.085 1.478±0.033 0.833±0.024 |
| 110 | + "VH" 1.29±0.71 … 0.0±0.0 0.0±0.0 |
| 111 | + "Others" 0.0568±0.0088 0.0021±0.0017 5.02±0.088 |
| 112 | + "Bkg Tot." 1233.6±5.4 557.2±2.4 164.9±1.2 |
| 113 | + "Significance" 0.3031±0.002 0.0454±0.00072 0.1672±0.0029 |
| 114 | +``` |
| 115 | +""" |
| 116 | +function significance_table(input_dir::AbstractString) |
| 117 | + body = significance_matrix(input_dir) |
| 118 | + significance_table(body) |
| 119 | +end |
| 120 | + |
| 121 | +function significance_table(body::Matrix; recreate=false) |
| 122 | + total_sig = body[1:1, :] #first row |
| 123 | + total_bkg = mapreduce(sum, hcat, eachcol(body[2:end, :])) #2:end row |
| 124 | + sig_errors = _significance.(total_sig, total_bkg) |
| 125 | + combined_sig = sqrt(sum(x[1]^2 for x in sig_errors)) |
| 126 | + |
| 127 | + full_nums = [ |
| 128 | + @. integral(body) ± only(binerrors(body)) |
| 129 | + @. integral(total_bkg) ± only(binerrors(total_bkg)) |
| 130 | + sig_errors |
| 131 | + ] |
| 132 | + |
| 133 | + full_body = [collect(ALL_TAGS) |
| 134 | + "Bkg Tot." |
| 135 | + "Significance" ;; full_nums |
| 136 | + ]; |
| 137 | +end |
| 138 | + |
| 139 | +const sigtable_fmt = (v, i, j) -> v isa Number ? "$(round(Measurements.value(v); digits=2)) ± $(round(Measurements.uncertainty(v); digits=2))" : v |
| 140 | + |
| 141 | +""" |
| 142 | + print_sigtable(full_table; io=stdout) |
| 143 | +
|
| 144 | +Takes the output of [`significance_table`](@ref) and pretty print it: |
| 145 | +
|
| 146 | +# Example |
| 147 | +
|
| 148 | +```julia |
| 149 | +julia> M = significance_table(<path>); |
| 150 | +
|
| 151 | +julia> print_sigtable(M) |
| 152 | +┌──────────────┬────────────────┬───────────────┬──────────────┬──────────────┬───────────────┐ |
| 153 | +│ │ SF-inZ │ SF-noZ │ DF │ CR-ZZ │ CR-ttZ │ |
| 154 | +├──────────────┼────────────────┼───────────────┼──────────────┼──────────────┼───────────────┤ |
| 155 | +│ Signal │ 10.66 ± 0.07 │ 9.31 ± 0.1 │ 10.73 ± 0.14 │ 1.07 ± 0.02 │ 2.15 ± 0.04 │ |
| 156 | +│ ZZ │ 1219.92 ± 5.38 │ 469.06 ± 2.44 │ 19.78 ± 0.45 │ 555.68 ± 2.4 │ 69.19 ± 0.74 │ |
| 157 | +│ Zjets │ -0.02 ± 0.13 │ 2.6 ± 2.22 │ 6.47 ± 5.51 │ -0.0 ± 0.0 │ 0.63 ± 0.39 │ |
| 158 | +│ Zgamma │ 0.0 ± 0.0 │ 0.0 ± 0.0 │ 0.3 ± 0.29 │ 0.0 ± 0.0 │ 0.0 ± 0.0 │ |
| 159 | +│ ttbar │ 0.0 ± 0.0 │ 0.63 ± 0.18 │ 0.28 ± 0.1 │ 0.0 ± 0.0 │ 0.43 ± 0.13 │ |
| 160 | +│ WZ │ 0.36 ± 0.1 │ 1.79 ± 0.23 │ 2.24 ± 0.29 │ 0.0 ± 0.0 │ 0.29 ± 0.06 │ |
| 161 | +│ tZ │ 0.01 ± 0.01 │ 0.07 ± 0.03 │ 0.06 ± 0.02 │ 0.0 ± 0.0 │ 0.18 ± 0.04 │ |
| 162 | +│ ttZ │ 1.23 ± 0.08 │ 4.71 ± 0.16 │ 5.74 ± 0.18 │ 0.01 ± 0.01 │ 73.62 ± 0.61 │ |
| 163 | +│ tWZ │ 0.56 ± 0.11 │ 2.16 ± 0.23 │ 2.5 ± 0.24 │ 0.01 ± 0.01 │ 14.71 ± 0.56 │ |
| 164 | +│ VBS │ 10.23 ± 0.09 │ 6.4 ± 0.08 │ 0.18 ± 0.01 │ 1.48 ± 0.03 │ 0.83 ± 0.02 │ |
| 165 | +│ VH │ 1.29 ± 0.71 │ 5.76 ± 1.4 │ 5.77 ± 1.29 │ 0.0 ± 0.0 │ 0.0 ± 0.0 │ |
| 166 | +├──────────────┼────────────────┼───────────────┼──────────────┼──────────────┼───────────────┤ |
| 167 | +│ Others │ 0.06 ± 0.01 │ 0.4 ± 0.13 │ 0.56 ± 0.08 │ 0.0 ± 0.0 │ 5.02 ± 0.09 │ |
| 168 | +├──────────────┼────────────────┼───────────────┼──────────────┼──────────────┼───────────────┤ |
| 169 | +│ Bkg Tot. │ 1233.64 ± 5.43 │ 493.58 ± 3.61 │ 43.89 ± 5.7 │ 557.19 ± 2.4 │ 164.91 ± 1.19 │ |
| 170 | +│ Significance │ 0.3 ± 0.0 │ 0.42 ± 0.0 │ 1.56 ± 0.1 │ 0.05 ± 0.0 │ 0.17 ± 0.0 │ |
| 171 | +└──────────────┴────────────────┴───────────────┴──────────────┴──────────────┴───────────────┘ |
| 172 | +``` |
| 173 | +""" |
| 174 | +print_sigtable(full_table; io=stdout) = pretty_table(io, |
| 175 | + full_table; |
| 176 | + header = ["", "SF-inZ", "SF-noZ", "DF", "CR-ZZ", "CR-ttZ"], |
| 177 | + formatters = sigtable_fmt, |
| 178 | + body_hlines = [size(full_table, 1) - 3, size(full_table, 1) - 2], |
| 179 | + highlighters = hl_row([1], crayon"bold"), |
| 180 | + alignment=:c, |
| 181 | +) |
| 182 | + |
| 183 | +end # module WVZReportExt |
0 commit comments