Add script to produce comparisons among RooUnfold and QUnfold

.gitignore

@@ -7,4 +7,5 @@ private
 **/.hypothesis
 **/delphes
 **/*.root
-data
+data
+**/img

requirements.txt

@@ -3,4 +3,6 @@ hypothesis
 # QUnfold
 tqdm
 uproot
-black
+black
+hist
+mplhep

src/PyXSec/core/Spectrum.py

@@ -159,7 +159,7 @@ def _configure(self):
             )
             self.ignore_background = True
         else:
-            log.info("Background path: {}".format(self.bkg_path))
+            log.info("Background path: {}".format(self.input_bkg_path))
 
         # Read histogram paths
         self.histo_reco_path = try_parse_str(self.histo_reco_path, root, "sig", "hpath")

src/PyXSec/core/Unfolder.py

@@ -124,7 +124,7 @@ def do_unfold(self, keep_response=False):
         # Unfolded distribution settings
         if self.error == "kNoError":
             if self.method == "SimNeal":
-                h_unfolded_array = self.m_unfolder.solve_simulated_annealing(num_reads=100)
+                h_unfolded_array = self.m_unfolder.solve_simulated_annealing(num_reads=50)
                 binning = [
                     self.h_data.GetXaxis().GetBinLowEdge(bin)
                     for bin in range(1, self.h_data.GetNbinsX() + 2)

studies/comparisons.py

@@ -0,0 +1,147 @@
+# ---------------------- Metadata ----------------------
+#
+# File name:  comparisons.py
+# Author:     Gianluca Bianco (biancogianluca9@gmail.com)
+# Date:       2023-11-26
+# Copyright:  (c) 2023 Gianluca Bianco under the MIT license.
+
+# TODO: i chi2 vanno calcolati con le covarianze
+# TODO: gira con dati veri
+# TODO: ratio plot
+
+# STD modules
+import argparse as ap
+import os
+
+# Data science modules
+import uproot
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import chisquare
+import math
+
+
+def compute_chi2(unfolded, truth):
+    # Trick for chi2 convergence
+    null_indices = truth == 0
+    truth[null_indices] += 1
+    unfolded[null_indices] += 1
+
+    # Compute chi2
+    chi2, _ = chisquare(
+        unfolded,
+        np.sum(unfolded) / np.sum(truth) * truth,
+    )
+    dof = len(unfolded) - 1
+    chi2_dof = chi2 / dof
+
+    return chi2_dof
+
+
+def plot_hist(qunfold, roounfold, theory):
+    # Plot QUnfold data
+    qunfold_val = qunfold.to_numpy()[0]
+    qunfold_err = qunfold.variances()
+    binning = qunfold.to_numpy()[1]
+    bin_midpoints = 0.5 * (binning[:-1] + binning[1:])
+    chi2_val = compute_chi2(qunfold.to_numpy()[0], theory.to_numpy()[0])
+    expo = math.floor(math.log10(abs(chi2_val))) - 2
+    chi2 = round(chi2_val, -expo)
+    label = rf"SimNeal ($\chi^2 = {chi2}$)"
+
+    plt.errorbar(
+        y=qunfold_val,
+        x=bin_midpoints,
+        yerr=qunfold_err,
+        ms=5,
+        fmt="o",
+        color="g",
+        label=label,
+    )
+
+    # Plot RooUnfold data
+    roounfold_val = roounfold.to_numpy()[0]
+    roounfold_err = roounfold.variances()
+    binning = roounfold.to_numpy()[1]
+
+    chi2_val = compute_chi2(roounfold.to_numpy()[0], theory.to_numpy()[0])
+    expo = math.floor(math.log10(abs(chi2_val))) - 2
+    chi2 = round(chi2_val, -expo)
+    label = rf"IBU ($\chi^2 = {chi2}$)"
+
+    plt.errorbar(
+        y=roounfold_val,
+        x=bin_midpoints,
+        yerr=roounfold_err,
+        ms=5,
+        fmt="v",
+        color="r",
+        label=label,
+    )
+
+    # Plot theory data
+    theory_val = theory.to_numpy()[0]
+    steps = np.append(theory_val, [theory_val[-1]])
+    plt.step(binning, steps, label="Theory", where="post")
+
+    # Plot settings
+    var_name = "{}".format(qunfold.title).replace("(particle)", "").replace("#", "\\")
+    plt.xlabel(
+        r"${}$".format(var_name),
+        loc="right",
+    )
+    if "Relative" in qunfold.name:
+        plt.ylabel(r"$1 / \sigma \cdot d\sigma / d{}$".format(var_name), loc="top")
+    elif "Absolute" in qunfold.name:
+        plt.ylabel(r"$d\sigma / d{}$".format(var_name), loc="top")
+    plt.title(qunfold.name)
+    plt.xlim(binning[0], binning[-1])
+    plt.ylim(0, plt.ylim()[1])
+    plt.legend(loc="upper right")
+
+    # Save
+    plt.tight_layout()
+    plt.savefig("studies/img/AbsoluteDiffXs.png")
+
+
+def main():
+    # Initial message
+    print("RooUnfold file: {}".format(args.roounfold))
+    print("QUnfold file: {}".format(args.qunfold))
+
+    # Create output dirs
+    if not os.path.exists("studies/img"):
+        os.makedirs("studies/img")
+
+    # Read QUnfold information
+    file_QUnfold = uproot.open(args.qunfold)
+    abs_Xs_QUnfold = file_QUnfold["AbsoluteDiffXs"]
+    abs_theory = file_QUnfold["TheoryXs_abs"]
+
+    # Read RooUnfold information
+    file_RooUnfold = uproot.open(args.roounfold)
+    abs_Xs_RooUnfold = file_RooUnfold["AbsoluteDiffXs"]
+
+    # Plot stuff
+    plot_hist(abs_Xs_QUnfold, abs_Xs_RooUnfold, abs_theory)
+
+
+if __name__ == "__main__":
+    # Parser settings
+    parser = ap.ArgumentParser(description="Parsing input arguments.")
+    parser.add_argument(
+        "-q",
+        "--qunfold",
+        default="",
+        help="Input root file with QUnfold results.",
+    )
+    parser.add_argument(
+        "-r",
+        "--roounfold",
+        default="",
+        help="Input root file with RooUnfold results.",
+    )
+    args = parser.parse_args()
+
+    # Main part
+    main()

tests/test_spectrum.py

@@ -14,8 +14,8 @@ def test_histo_comparison():
     Perform comparison among each histogram of the correct output file (produced by TTbarUnfold) and the output file produced by PyXSec.
     """
 
-    file_output = "data/private/RooUnfold/output.root"
-    file_correct = "data/private/RooUnfold/output_correct.root"
+    file_output = "data/private/other/output_RooUnfold.root"
+    file_correct = "data/private/other/output_correct.root"
 
     # Open ROOT files
     root_output = ROOT.TFile(file_output, "READ")
@@ -54,14 +54,14 @@ def test_histo_comparison():
         # Check each bin entry
         if hist_output.GetName() != "Response":
             for bin in range(1, hist_output.GetNbinsX() + 1):
-                assert round(hist_output.GetBinContent(bin), 6) == round(
-                    hist_correct.GetBinContent(bin), 6
+                assert round(hist_output.GetBinContent(bin), 4) == round(
+                    hist_correct.GetBinContent(bin), 4
                 )
         else:
             for bin_x in range(1, hist_output.GetNbinsX() + 1):
                 for bin_y in range(1, hist_output.GetNbinsY() + 1):
-                    assert round(hist_output.GetBinContent(bin_x, bin_y), 6) == round(
-                        hist_correct.GetBinContent(bin_x, bin_y), 6
+                    assert round(hist_output.GetBinContent(bin_x, bin_y), 4) == round(
+                        hist_correct.GetBinContent(bin_x, bin_y), 4
                     )
 
     # Close files

tox.ini

@@ -2,7 +2,7 @@
 [tox]
 requires =
     tox>=4
-env_list = example_RooUnfold, example_QUnfold, tests, data_processing, py{38}
+env_list = example_RooUnfold, example_QUnfold, tests, data_processing, comparisons, py{38}
 
 # Base tox configuration
 [testenv]
@@ -23,8 +23,9 @@ deps =
 commands =
     {[testenv]commands}
     python src/PyXSec/PyXSec.py \
-    --config="data/private/RooUnfold/config.xml" \
-    --output="data/private/RooUnfold/output.root"
+    --config="data/private/pT_lep1_MCstat/config_RooUnfold.xml" \
+    --systematic="MCstat" \
+    --output="data/private/pT_lep1_MCstat/output_RooUnfold.root"
 
 # Run example QUnfold configuration
 [testenv:example_QUnfold]
@@ -34,8 +35,9 @@ deps =
 commands =
     {[testenv]commands}
     python src/PyXSec/PyXSec.py \
-    --config="data/private/QUnfold/config.xml" \
-    --output="data/private/QUnfold/output.root"
+    --config="data/private/pT_lep1_MCstat/config_QUnfold.xml" \
+    --systematic="MCstat" \
+    --output="data/private/pT_lep1_MCstat/output_QUnfold.root"
 
 # Run tests configuration
 [testenv:tests]
@@ -60,4 +62,16 @@ commands =
     --particle="data/simulated/particle.root" \
     --reco-response="data/simulated/reco.root" \
     --particle-response="data/simulated/particle.root" \
-    --output="data/simulated/unfolding_input.root"
+    --output="data/simulated/unfolding_input.root"
+
+# Run comparisons configuration
+[testenv:comparisons]
+description = run comparisons part
+allowlist_externals = find
+deps =
+    {[testenv]deps}
+commands =
+    {[testenv]commands}
+    python studies/comparisons.py \
+    --qunfold="data/private/pT_lep1_MCstat/output_QUnfold.root" \
+    --roounfold="data/private/pT_lep1_MCstat/output_RooUnfold.root"