Extend th covmat plots

validphys2/src/validphys/pdfgrids.py

@@ -374,6 +374,42 @@ def distance_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) =
 
     return newgrids
 
+@check_pdf_normalize_to
+def pull_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = None):
+    """Return an object containing the value of the pull between the two PDFs at the specified values
+    of x and flavour.
+    The parameter ``normalize_to`` identifies the reference PDF set with respect to the
+    pull is computed.
+    This method returns pull grids where the relative pull between both PDF
+    set is computed. At least one grid will be identical to zero.
+    """
+
+    gr2_stats = xplotting_grids[normalize_to].grid_values
+    cv2 = gr2_stats.central_value()
+    sg2 = gr2_stats.std_error()
+
+    newgrids = []
+    for grid, pdf in zip(xplotting_grids, pdfs):
+        if pdf == pdfs[normalize_to]:
+            # Zero the PDF we are normalizing against
+            pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0:1]))
+            newgrid = grid.copy_grid(grid_values=pdf_zero)
+            newgrids.append(newgrid)
+            continue
+
+        g_stats = grid.grid_values
+        cv1 = g_stats.central_value()
+        sg1 = g_stats.std_error()
+
+        # Wrap the pull into a Stats (1, flavours, points)
+        pull = Stats([np.sqrt((cv1 - cv2) ** 2 / (sg1**2 + sg2**2))])
+
+        newgrid = grid.copy_grid(grid_values=pull)
+        newgrids.append(newgrid)
+
+    return newgrids
+
+pull_grids_list = collect(pull_grids, ('pdfs_list',))
 
 @check_pdf_normalize_to
 def variance_distance_grids(pdfs, xplotting_grids, normalize_to: (int, str, type(None)) = None):

validphys2/src/validphys/pdfplots.py

@@ -273,6 +273,104 @@ def plot_pdf_uncertainties(
     PDF's central value is plotted. Otherwise it is the absolute values."""
     yield from UncertaintyPDFPlotter(pdfs, xplotting_grids, xscale, normalize_to, ymin, ymax)
 
+class PullPDFPlotter(metaclass=abc.ABCMeta):
+    """Auxiliary class which groups multiple pulls in one plot."""
+
+    def __init__(self, pdfs_list, pull_grids_list, xscale, normalize_to, ymin, ymax):
+        self.pdfs_list = pdfs_list
+        self.pull_grids_list = pull_grids_list
+        self._xscale = xscale
+        self.normalize_to = normalize_to
+        self.ymin = ymin
+        self.ymax = ymax
+        self.firstgrid = pull_grids_list[0][0]
+
+    def legend(self, flstate):
+        return flstate.ax.legend()
+
+    def get_ylabel(self):
+        return "Pull"
+
+    @property
+    def Q(self):
+        return self.firstgrid.Q
+
+    @property
+    def xscale(self):
+        if self._xscale is None:
+            return scale_from_grid(self.firstgrid)
+        return self._xscale
+
+    def get_title(self, flstate):
+        return '$%s$' % flstate.parton_name + f', Q={self.Q : .1f} GeV'
+
+    def draw(self, pdfs, grid, flstate):
+        ax = flstate.ax
+        flindex = flstate.flindex
+        color = ax._get_lines.get_next_color()
+
+        # The grid for the pull is (1, flavours, points)
+        # take only the flavour we are interested in
+        gv = grid.select_flavour(flindex).grid_values.data.squeeze()
+
+        ax.plot(grid.xgrid, gv, color=color, label = f'{pdfs[0].label}-{pdfs[1].label} pull')
+
+        return gv
+
+    def plot_call(self):
+        basis = self.firstgrid.basis
+        for flindex, fl in enumerate(self.firstgrid.flavours):
+            fig, ax = plotutils.subplots()
+            parton_name = basis.elementlabel(fl)
+            flstate = FlavourState(flindex=flindex, fl=fl, fig=fig, ax=ax, parton_name=parton_name)
+            ax.set_title(self.get_title(flstate))
+
+            all_vals = []
+            for pdf, grids in zip(self.pdfs_list, self.pull_grids_list):
+                limits = self.draw([pdf['pdfs'][0],pdf['pdfs'][1]], grids[1], flstate)
+                if limits is not None:
+                    all_vals.append(np.atleast_2d(limits))
+
+            # Note these two lines do not conmute!
+            ax.set_xscale(self.xscale)
+            plotutils.frame_center(ax, self.firstgrid.xgrid, np.concatenate(all_vals))
+            if self.ymin is not None:
+                ax.set_ylim(ymin=self.ymin)
+            if self.ymax is not None:
+                ax.set_ylim(ymax=self.ymax)
+
+            ax.set_xlabel('$x$')
+            ax.set_xlim(self.firstgrid.xgrid[0])
+
+            ax.set_ylabel(self.get_ylabel())
+
+            ax.set_axisbelow(True)
+
+            self.legend(flstate)
+            yield fig, parton_name
+
+    def __call__(self):
+        for fig, parton_name in self.plot_call():
+            ax = fig.get_axes()[0]
+            ymin, _ = ax.get_ylim()
+            ax.set_ylim(max(0, ymin), None)
+            yield fig, parton_name
+
+
+@figuregen
+@check_scale('xscale', allow_none=True)
+def plot_pdf_pulls(
+    pdfs_list,
+    pull_grids_list,
+    xscale: (str, type(None)) = None,
+    normalize_to: (int, str, type(None)) = None,
+    ymin=None,
+    ymax=None,
+):
+    """Plot the PDF standard deviations as a function of x.
+    If normalize_to is set, the ratio to that
+    PDF's central value is plotted. Otherwise it is the absolute values."""
+    yield from PullPDFPlotter(pdfs_list, pull_grids_list, xscale, normalize_to, ymin, ymax)()
 
 class AllFlavoursPlotter(PDFPlotter):
     """Auxiliary class which groups multiple PDF flavours in one plot."""

validphys2/src/validphys/results.py

@@ -291,10 +291,17 @@ def procs_data_values(proc_result_table):
     return data_central_values
 
 
+def procs_data_values_experiment(proc_result_table_experiment):
+    """Like groups_data_values but grouped by experiment."""
+    data_central_values = proc_result_table_experiment["data_central"]
+    return data_central_values
+
+
 groups_results = collect("dataset_inputs_results", ("group_dataset_inputs_by_metadata",))
 
-procs_results = collect("dataset_inputs_results", ("group_dataset_inputs_by_process",))
+procs_results = collect("dataset_inputs_results_central", ("group_dataset_inputs_by_process",))
 
+procs_results_experiment = collect("dataset_inputs_results_central", ("group_dataset_inputs_by_experiment",),)
 
 def group_result_table_no_table(groups_results, groups_index):
     """Generate a table containing the data central value, the central prediction,
@@ -328,11 +335,14 @@ def group_result_table(group_result_table_no_table):
 def proc_result_table_no_table(procs_results, procs_index):
     return group_result_table_no_table(procs_results, procs_index)
 
-
 @table
 def proc_result_table(proc_result_table_no_table):
     return proc_result_table_no_table
 
+@table
+def proc_result_table_experiment(procs_results_experiment, experiments_index):
+    return group_result_table_no_table(procs_results_experiment, experiments_index)
+
 
 experiment_result_table = collect("group_result_table", ("group_dataset_inputs_by_experiment",))
 
@@ -566,6 +576,12 @@ def results_with_scale_variations(results, theory_covmat_dataset):
     theory_error_result = ThUncertaintiesResult(cv, total_error, label=central_th_result.label)
     return (data_result, theory_error_result)
 
+def dataset_inputs_results_central(
+    data, pdf: PDF, dataset_inputs_covariance_matrix, dataset_inputs_sqrt_covmat
+):
+    """Like `dataset_inputs_results` but for a group of datasets and replica0."""
+    return results_central(data, pdf, dataset_inputs_covariance_matrix, dataset_inputs_sqrt_covmat)
+
 
 def dataset_inputs_results(
     data, pdf: PDF, dataset_inputs_covariance_matrix, dataset_inputs_sqrt_covmat

validphys2/src/validphys/scalevariations/pointprescriptions.yaml

@@ -21,3 +21,7 @@
 'n3lo 3pt hadronic': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)', '(1, 0.5 hadronic)', '(1, 2 hadronic)','(-1, -1)','(1, 1)']
 # N3LO 7 point scale variations
 'n3lo 7 point': ['(0, 0, 0, 0)', '(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)']
+# N3LO 3 point hadronic renormalization scale variations
+'n3lo 3r point hadronic': ['(0, 0, 0, 0)', '(1, 0.5 hadronic)', '(1, 2 hadronic)']
+# N3LO 3 point missing renormalization scale variations
+'n3lo 3r point missing': ['(0, 0, 0, 0)', '(1, 0.5 missing)', '(1, 2 missing)']

validphys2/src/validphys/tableloader.py

@@ -41,6 +41,11 @@ def parse_data_cv(filename):
     df = sane_load(filename, index_col=[0, 1, 2])
     return df
 
+def parse_chi2_data(filename):
+    """Useful for reading whole DataFrames."""
+    df = sane_load(filename)
+    return df
+
 
 def parse_exp_mat(filename):
     """Parse a dump of a matrix like experiments_covmat."""

validphys2/src/validphys/tests/conftest.py

@@ -50,6 +50,14 @@ def tmp(tmpdir):
     {'dataset': 'NMC', 'weight': 100},
 ]
 
+DATA_THCOVMAT = [
+    {'dataset': 'NMC'},
+    {'dataset': 'CHORUSNUPb_dw_ite'},
+    {'dataset': 'CMSZDIFF12', 'cfac':['NRM'], 'sys':10},
+    {'dataset': 'ATLAS_WP_JET_8TEV_PT'},
+    {'dataset': 'LHCBZMU8TEV', 'cfac':['NRM']}
+]
+
 POSITIVITIES = ["POSDYCBD", "POSF2S"]
 
 PDF = "NNPDF40_nnlo_as_01180"
@@ -77,6 +85,17 @@ def tmp(tmpdir):
 def data_config():
     return base_config
 
+@pytest.fixture(scope='module')
+def thcovmat_config(data_config):
+    """Same as data_config but with additional info for the thcovmat production."""
+    new_config = dict(data_config)
+    new_config["point_prescription"] = "3 point"
+    new_config["use_theorycovmat"] = "true"
+    new_config["use_cuts"] = "internal"
+    new_config.update(theoryid=708)
+    new_config["theoryids"] = {"from_": "scale_variation_theories"}
+    new_config.update(dataset_inputs=DATA_THCOVMAT)
+    return new_config
 
 @pytest.fixture(scope='module')
 def data_internal_cuts_config(data_config):

validphys2/src/validphys/tests/regressions/test_thcovmat_chi2.csv

@@ -0,0 +1,6 @@
+	dataset	chi2_value
+0	NMC	278.5292049671733
+1	CHORUSNUPb_dw_ite	400.4162663033812
+2	CMSZDIFF12	38.82979671473409
+3	ATLAS_WP_JET_8TEV_PT	11.501661474293808
+4	LHCBZMU8TEV	22.882044849141472

validphys2/src/validphys/tests/test_regressions.py

@@ -18,7 +18,7 @@
 from validphys.api import API
 from validphys.tests.test_covmats import CORR_DATA
 from validphys.tableloader import (parse_data_cv, parse_exp_mat, load_perreplica_chi2_table,
-                                   sane_load, load_fits_chi2_table)
+                                   sane_load, load_fits_chi2_table, parse_chi2_data)
 from validphys.tests.test_covmats import CORR_DATA
 
 
@@ -67,6 +67,14 @@ def test_mcreplica(data_config):
 def test_expcovmat(data_config):
     return API.groups_covmat_no_table(**data_config)
 
+@make_table_comp(parse_chi2_data)
+def test_thcovmat_chi2(thcovmat_config):
+    chi2_by_dataset = API.each_dataset_chi2_sv(**thcovmat_config)
+    records = []
+    for dataset, chi2 in zip(thcovmat_config['dataset_inputs'], chi2_by_dataset):
+        records.append(dict(dataset=dataset['dataset'], chi2_value=chi2.central_result))
+    return pd.DataFrame.from_records(records)
+
 
 @make_table_comp(parse_exp_mat)
 def test_t0covmat(data_witht0_internal_cuts_config):