From 4b58c061e9813d681a3947d1add14337b7edf045 Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Sun, 28 Dec 2025 20:20:07 +0000
Subject: [PATCH 01/14] Edited loading fk tables

---
 validphys2/src/validphys/config.py  |  8 +++++++-
 validphys2/src/validphys/core.py    |  7 ++++++-
 validphys2/src/validphys/covmats.py | 17 ++++++++++++++++-
 validphys2/src/validphys/results.py |  9 ++++++++-
 4 files changed, 37 insertions(+), 4 deletions(-)

diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py
index 4b2d8cad5b..ece19dc466 100644
--- a/validphys2/src/validphys/config.py
+++ b/validphys2/src/validphys/config.py
@@ -347,7 +347,13 @@ def produce_fitunderlyinglaw(self, fit):
         """
         with self.set_context(ns=self._curr_ns.new_child({"fit": fit})):
             _, datacuts = self.parse_from_("fit", "closuretest", write=False)
-        underlyinglaw = self.parse_pdf(datacuts["fakepdf"])
+        # underlyinglaw = self.parse_pdf(datacuts["fakepdf"])
+        fakepdf = datacuts["fakepdf"]
+
+        if isinstance(fakepdf, str):
+            underlyinglaw = self.parse_pdf(fakepdf)
+        else:
+            underlyinglaw = fakepdf
         return {"pdf": underlyinglaw}
 
     @element_of("hyperscans")
diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index 09e7516fd6..61bde31dbc 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -42,6 +42,7 @@
 from validphys.hyperoptplot import HyperoptTrial
 from validphys.utils import experiments_to_dataset_inputs
 from validphys.lhapdfset import LHAPDFSet
+# from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
@@ -513,7 +514,11 @@ def load(self):
 
         fktables = []
         for p in self.fkspecs:
-            fktable = p.load()
+            try:
+                fktable = p.load()
+            except Exception as e:
+                from validphys.pineparser import pineappl_reader
+                fktable = pineappl_reader(p)
             #IMPORTANT: We need to tell the python garbage collector to NOT free the
             #memory owned by the FKTable on garbage collection.
             #TODO: Do this automatically
diff --git a/validphys2/src/validphys/covmats.py b/validphys2/src/validphys/covmats.py
index f710f7f9af..ecf8561b79 100644
--- a/validphys2/src/validphys/covmats.py
+++ b/validphys2/src/validphys/covmats.py
@@ -410,8 +410,23 @@ def sqrt_covmat(covariance_matrix):
                          f"{dimensions[1]}")
 
     sqrt_diags = np.sqrt(np.diag(covariance_matrix))
-    correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
+    # There are some zero entries in sqrt_diags, this gives errors when finding the cholesky decomposition.
+
+    if np.any(sqrt_diags == 0):
+        # Use safe method if there are zeros
+        outer = np.outer(sqrt_diags, sqrt_diags)
+        correlation_matrix = np.divide(
+            covariance_matrix, outer, out=np.zeros_like(covariance_matrix), where=outer!=0
+        )
+    else:
+        # Normal division if all variances are nonzero
+        correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
+
+    # Always fix the diagonal
+    np.fill_diagonal(correlation_matrix, 1.0)
+    # correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
     decomp = la.cholesky(correlation_matrix)
+
     sqrt_matrix = (decomp * sqrt_diags).T
     return sqrt_matrix
 
diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index 48a9f52220..e7b9a7eb66 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -42,6 +42,7 @@
 
 from validphys.n3fit_data_utils import parse_simu_parameters_names_CF
 
+from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
@@ -385,7 +386,13 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
     if parsed_bsm_facs is None:
         # We want an array of ones that ndata x nrep
         # where ndata is the number of post cut datapoints
-        ndata = len(dataset.load().get_cv())
+        try:
+            ndata = len(dataset.load().get_cv())
+        except Exception:
+            fkspec = dataset.fkspecs[0]
+            fkdata = pineappl_reader(fkspec)
+            ndata = fkdata.ndata
+
         nrep = len(pdf)
         return np.ones((ndata, nrep))
 

From 0af4ae50fc26521a485feccac5cd6d65b53cc37a Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Wed, 31 Dec 2025 16:10:29 +0000
Subject: [PATCH 02/14] Added level0 bsm chi2 function

---
 validphys2/src/validphys/core.py             |   6 +-
 validphys2/src/validphys/simunet_analysis.py | 111 +++++++++++++++++++
 2 files changed, 112 insertions(+), 5 deletions(-)

diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index 61bde31dbc..d0924a577d 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -514,11 +514,7 @@ def load(self):
 
         fktables = []
         for p in self.fkspecs:
-            try:
-                fktable = p.load()
-            except Exception as e:
-                from validphys.pineparser import pineappl_reader
-                fktable = pineappl_reader(p)
+            fktable = p.load()
             #IMPORTANT: We need to tell the python garbage collector to NOT free the
             #memory owned by the FKTable on garbage collection.
             #TODO: Do this automatically
diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 601a153286..8257160978 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2107,6 +2107,117 @@ def compute_datasets_chi2_dist(
 
     return chi2_dict
 
+def compute_datasets_chi2(
+    level0_commondata_wc,
+    sm_predictions,
+    groups_covmat,
+    load_datasets_contamination,
+    read_bsm_facs,
+    dataset_inputs,
+    theoryid,
+):
+    """
+    Parameters
+    ----------
+
+    level0_commondata_wc: level0 data using 'fakepdf'
+
+    sm_predictions: SM predictions using 'pdf'
+
+    groups_covmat
+
+    load_contamination
+
+    read_bsm_facs: BSM factors from the fit
+
+    dataset_inputs
+
+    theoryid
+
+    Returns
+    -------
+
+    dict
+        dictionary of lists of chi2 per dataset
+
+    """
+
+    bsm_facs_df = read_bsm_facs
+
+    means = bsm_facs_df.mean()
+    central_pred = {}
+    if dataset_inputs is not None:
+        for dataset in dataset_inputs:
+            central_sm = sm_predictions[dataset.name]
+            bsm_factors = np.zeros(len(central_sm))
+            # dataset.simu_parameters_linear_combinations includes the contamination linear combinations
+            # This loads the whole dataset but we only need simu_facs - is there a way to clean this?
+            ds = l.check_dataset(
+                name=dataset.name,
+                theoryid=theoryid,
+                cfac=dataset.cfac,
+                simu_parameters_names=dataset.simu_parameters_names,
+                simu_parameters_linear_combinations=dataset.simu_parameters_linear_combinations,
+                use_fixed_predictions=dataset.use_fixed_predictions,
+                new_commondata=dataset.new_commondata,
+            )
+            bsm_fac = parse_simu_parameters_names_CF(
+                ds.simu_parameters_names_CF,
+                ds.simu_parameters_linear_combinations,
+                cuts=ds.cuts,
+            )
+            # bsm_fac = (contamination_value*k-factors)/SM pred in Simu_fac file
+
+            if bsm_fac != None:
+                coefficients = central_sm.to_numpy().T * np.array(
+                    [i.central_value for i in bsm_fac.values()]
+                )
+                for i, key in enumerate(bsm_fac.keys()):
+                    label = key.split("_")[-1]
+
+                    scaled_row = coefficients[i] * means[label]
+
+                    bsm_factors += scaled_row
+
+            central_pred[dataset.name] = central_sm.values.squeeze() * (1 + bsm_factors)
+
+    covmat = groups_covmat
+    data = level0_commondata_wc
+    contamination_factors = load_datasets_contamination
+
+    chi2_dict = {dataset.setname: [] for dataset in data}
+
+    for dataset in data:
+        data_name = dataset.setname
+        cont_fac = contamination_factors[data_name]
+
+        if cont_fac.shape[0] == 1:
+            data_values = dataset.central_values * cont_fac
+        else:
+            indices = dataset.commondata_table_indices
+            data_values = dataset.central_values * cont_fac[indices]
+
+        num_data = dataset.ndata
+
+        covmat_dataset = (
+            covmat.xs(data_name, level=1, drop_level=False)
+            .T.xs(data_name, level=1, drop_level=False)
+            .values
+        )
+
+        theory = central_pred[data_name]
+
+        diff = (data_values - theory).squeeze()
+
+        if diff.size == 1:
+            chi2 = diff**2 / covmat_dataset[0, 0] / num_data
+        else:
+            chi2 = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) / num_data
+
+        chi2_dict[data_name].append(chi2)
+
+    return chi2_dict
+
 
 def write_datasets_chi2_dist_csv(
         pdf,

From 1c366be02ec657273d22eae45f32066811f2a364 Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Wed, 31 Dec 2025 16:39:11 +0000
Subject: [PATCH 03/14] Edited contamination function to include multiple
 parameters

---
 validphys2/src/validphys/simunet_analysis.py | 58 +++++++++++++++++---
 1 file changed, 50 insertions(+), 8 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 8257160978..413e48d126 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2001,10 +2001,6 @@ def load_datasets_contamination(
 
     cont_path = l.datapath / f"theory_{theoryid.id}" / "simu_factors"
 
-    cont_name = contamination_parameters["name"]
-    cont_value = contamination_parameters["value"]
-    cont_lin_comb = contamination_parameters["linear_combination"]
-
     bsm_dict = {}
 
     for dataset in dataset_inputs:
@@ -2032,11 +2028,26 @@ def load_datasets_contamination(
             stream.close()
 
             k_factors = np.zeros(len(simu_card["SM_fixed"]))
-            for op in cont_lin_comb:
-                k_factors += cont_lin_comb[op] * np.array(simu_card[cont_order][op])
-            k_factors = 1. + k_factors * cont_value / np.array(simu_card[cont_order]["SM"])
 
-            bsm_dict[dataset.name] = k_factors
+            for cont_params in contamination_parameters:
+                cont_name = cont_params["name"]
+                cont_value = cont_params["value"]
+                cont_lin_comb = cont_params["linear_combination"]
+
+                if cont_value == 0.:
+                    continue
+
+                for op in cont_lin_comb:
+                    if op in simu_card[cont_order]:
+                        k_factors += cont_lin_comb[op] * np.array(simu_card[cont_order][op])
+                        
+
+            sm_array = np.array(simu_card[cont_order]["SM"])
+            total_factor = 1. + (k_factors*cont_value) / sm_array
+            print(f"Dataset {dataset.name} contaminated with factor {cont_value} on order {cont_order}")
+            print(f"Total BSM factor min: {total_factor.min()}, max: {total_factor.max()}")
+
+            bsm_dict[dataset.name] = total_factor
 
     return bsm_dict
 
@@ -2248,3 +2259,34 @@ def write_datasets_chi2_dist_csv(
     chi2 = pd.concat([df_ndat, df_chi2], ignore_index=True)
 
     chi2.to_csv(f"{pdf}_chi2_dist.csv", index=False)
+
+
+def write_datasets_chi2_csv(
+        pdf,
+        compute_datasets_chi2,
+        level0_commondata_wc
+    ):
+
+    """
+    Parameters
+    ----------
+
+    pdf: core.PDF
+
+    compute_chi2
+
+    level0_commondata_wc
+
+    Returns
+    -------
+    
+    """
+
+    ndata = {dataset.setname: [dataset.ndata] for dataset in level0_commondata_wc}
+
+    df_ndat = pd.DataFrame(ndata)
+    df_chi2 = pd.DataFrame(compute_datasets_chi2)
+
+    chi2 = pd.concat([df_ndat, df_chi2], ignore_index=True)
+
+    chi2.to_csv(f"{pdf}_chi2_dist.csv", index=False)
\ No newline at end of file

From 626c7f342f491d4d2b7585a2b4e0dd5d22e68baa Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Wed, 31 Dec 2025 16:50:00 +0000
Subject: [PATCH 04/14] revert files

---
 validphys2/src/validphys/core.py    |  1 -
 validphys2/src/validphys/covmats.py | 17 +----------------
 validphys2/src/validphys/results.py |  9 +--------
 3 files changed, 2 insertions(+), 25 deletions(-)

diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index d0924a577d..09e7516fd6 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -42,7 +42,6 @@
 from validphys.hyperoptplot import HyperoptTrial
 from validphys.utils import experiments_to_dataset_inputs
 from validphys.lhapdfset import LHAPDFSet
-# from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
diff --git a/validphys2/src/validphys/covmats.py b/validphys2/src/validphys/covmats.py
index ecf8561b79..f710f7f9af 100644
--- a/validphys2/src/validphys/covmats.py
+++ b/validphys2/src/validphys/covmats.py
@@ -410,23 +410,8 @@ def sqrt_covmat(covariance_matrix):
                          f"{dimensions[1]}")
 
     sqrt_diags = np.sqrt(np.diag(covariance_matrix))
-    # There are some zero entries in sqrt_diags, this gives errors when finding the cholesky decomposition.
-
-    if np.any(sqrt_diags == 0):
-        # Use safe method if there are zeros
-        outer = np.outer(sqrt_diags, sqrt_diags)
-        correlation_matrix = np.divide(
-            covariance_matrix, outer, out=np.zeros_like(covariance_matrix), where=outer!=0
-        )
-    else:
-        # Normal division if all variances are nonzero
-        correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
-
-    # Always fix the diagonal
-    np.fill_diagonal(correlation_matrix, 1.0)
-    # correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
+    correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
     decomp = la.cholesky(correlation_matrix)
-
     sqrt_matrix = (decomp * sqrt_diags).T
     return sqrt_matrix
 
diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index e7b9a7eb66..48a9f52220 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -42,7 +42,6 @@
 
 from validphys.n3fit_data_utils import parse_simu_parameters_names_CF
 
-from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
@@ -386,13 +385,7 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
     if parsed_bsm_facs is None:
         # We want an array of ones that ndata x nrep
         # where ndata is the number of post cut datapoints
-        try:
-            ndata = len(dataset.load().get_cv())
-        except Exception:
-            fkspec = dataset.fkspecs[0]
-            fkdata = pineappl_reader(fkspec)
-            ndata = fkdata.ndata
-
+        ndata = len(dataset.load().get_cv())
         nrep = len(pdf)
         return np.ones((ndata, nrep))
 

From 03912be7bb391e3f559597b7ad98d63b06df54fb Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Mon, 5 Jan 2026 17:17:07 +0000
Subject: [PATCH 05/14] Added PDF fit covmat

---
 validphys2/src/validphys/results.py          |  8 ++++++--
 validphys2/src/validphys/simunet_analysis.py | 21 ++++++++++++++++----
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index 48a9f52220..a847fb09c9 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -36,6 +36,7 @@
 from validphys.convolution import (
     predictions,
     PredictionsRequireCutsError,
+    central_predictions
 )
 from validphys.plotoptions.core import get_info
 
@@ -385,7 +386,10 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
     if parsed_bsm_facs is None:
         # We want an array of ones that ndata x nrep
         # where ndata is the number of post cut datapoints
-        ndata = len(dataset.load().get_cv())
+        try:
+            ndata = len(dataset.load().get_cv())
+        except Exception as e:
+            ndata = (len(central_predictions(dataset, pdf)))
         nrep = len(pdf)
         return np.ones((ndata, nrep))
 
@@ -408,7 +412,7 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
     if not read_bsm_facs.empty:
         replica_result = 1 + np.sum(scaled_replicas, axis=2)
     else:
-        replica_result = np.ones((len(dataset.load().get_cv()), len(pdf)-1))
+        replica_result = np.ones((ndata, len(pdf)-1))
     
     average_result = np.mean(replica_result, axis=1, keepdims=True)
     result = np.concatenate((average_result, replica_result), axis=1)
diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 413e48d126..2069fc4013 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -42,7 +42,8 @@
 from validphys.n3fit_data_utils import parse_simu_parameters_names_CF
 from validphys.loader import _get_nnpdf_profile
 
-from validphys.convolution import central_predictions
+from validphys.convolution import central_predictions, predictions
+from validphys.results import ThPredictionsResult, dataset_bsm_factor
 
 log = logging.getLogger(__name__)
 
@@ -2119,6 +2120,7 @@ def compute_datasets_chi2_dist(
     return chi2_dict
 
 def compute_datasets_chi2(
+    pdf,
     level0_commondata_wc,
     sm_predictions,
     groups_covmat,
@@ -2156,7 +2158,9 @@ def compute_datasets_chi2(
     bsm_facs_df = read_bsm_facs
 
     means = bsm_facs_df.mean()
+    stds = bsm_facs_df.std()
     central_pred = {}
+    pdf_covmats = {}
     if dataset_inputs is not None:
         for dataset in dataset_inputs:
             central_sm = sm_predictions[dataset.name]
@@ -2190,9 +2194,18 @@ def compute_datasets_chi2(
 
                     bsm_factors += scaled_row
 
-            central_pred[dataset.name] = central_sm.values.squeeze() * (1 + bsm_factors)
+            central_pred[dataset.name] = central_sm.values.squeeze() * (1 + bsm_factors) 
+            # Do we want to use this or rep_bsm_predictions which is the mean over all replicas?
+
+            # Finding PDF/SMEFT fit covmat from replicas
+            rep_sm_predictions = predictions(ds, pdf)
+            bsm_factor=dataset_bsm_factor(ds,pdf,read_bsm_facs)
+            rep_bsm_predictions = rep_sm_predictions * bsm_factor
+            replicas = rep_bsm_predictions.iloc[:, 1:]
+            pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True)
+            
+    covmat = groups_covmat # This is the experimental covmat
 
-    covmat = groups_covmat
     data = level0_commondata_wc
     contamination_factors = load_datasets_contamination
 
@@ -2214,7 +2227,7 @@ def compute_datasets_chi2(
             covmat.xs(data_name, level=1, drop_level=False)
             .T.xs(data_name, level=1, drop_level=False)
             .values
-        )
+        ) + pdf_covmats[data_name]
 
         theory = central_pred[data_name]
 

From 6280484aa45928f654c0bda1e72aa31e64c33f1c Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Tue, 6 Jan 2026 14:01:12 +0000
Subject: [PATCH 06/14] Full covmat

---
 validphys2/src/validphys/simunet_analysis.py | 27 +-------------------
 1 file changed, 1 insertion(+), 26 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 2069fc4013..db857318eb 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2045,8 +2045,6 @@ def load_datasets_contamination(
 
             sm_array = np.array(simu_card[cont_order]["SM"])
             total_factor = 1. + (k_factors*cont_value) / sm_array
-            print(f"Dataset {dataset.name} contaminated with factor {cont_value} on order {cont_order}")
-            print(f"Total BSM factor min: {total_factor.min()}, max: {total_factor.max()}")
 
             bsm_dict[dataset.name] = total_factor
 
@@ -2163,10 +2161,6 @@ def compute_datasets_chi2(
     pdf_covmats = {}
     if dataset_inputs is not None:
         for dataset in dataset_inputs:
-            central_sm = sm_predictions[dataset.name]
-            bsm_factors = np.zeros(len(central_sm))
-            # dataset.simu_parameters_linear_combinations includes the contamination linear combinations
-            # This loads the whole dataset but we only need simu_facs - is there a way to clean this?
             ds = l.check_dataset(
                 name=dataset.name,
                 theoryid=theoryid,
@@ -2176,26 +2170,6 @@ def compute_datasets_chi2(
                 use_fixed_predictions=dataset.use_fixed_predictions,
                 new_commondata=dataset.new_commondata,
             )
-            bsm_fac = parse_simu_parameters_names_CF(
-                ds.simu_parameters_names_CF,
-                ds.simu_parameters_linear_combinations,
-                cuts=ds.cuts,
-            )
-            # bsm_fac = (contamination_value*k-factors)/SM pred in Simu_fac file
-
-            if bsm_fac != None:
-                coefficients = central_sm.to_numpy().T * np.array(
-                    [i.central_value for i in bsm_fac.values()]
-                )
-                for i, key in enumerate(bsm_fac.keys()):
-                    label = key.split("_")[-1]
-
-                    scaled_row = coefficients[i] * means[label]
-
-                    bsm_factors += scaled_row
-
-            central_pred[dataset.name] = central_sm.values.squeeze() * (1 + bsm_factors) 
-            # Do we want to use this or rep_bsm_predictions which is the mean over all replicas?
 
             # Finding PDF/SMEFT fit covmat from replicas
             rep_sm_predictions = predictions(ds, pdf)
@@ -2203,6 +2177,7 @@ def compute_datasets_chi2(
             rep_bsm_predictions = rep_sm_predictions * bsm_factor
             replicas = rep_bsm_predictions.iloc[:, 1:]
             pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True)
+            central_pred[dataset.name] = rep_bsm_predictions.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
             
     covmat = groups_covmat # This is the experimental covmat
 

From f67c9c717a126727b33ebb22d905116e45f4ead2 Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Thu, 8 Jan 2026 11:34:16 +0000
Subject: [PATCH 07/14] T0 and exp comparison

---
 validphys2/src/validphys/simunet_analysis.py | 66 +++++++++++++-------
 1 file changed, 43 insertions(+), 23 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index db857318eb..29bcef423e 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2126,6 +2126,8 @@ def compute_datasets_chi2(
     read_bsm_facs,
     dataset_inputs,
     theoryid,
+    dataset_inputs_covmat_t0_considered,
+    simu_parameters_scales
 ):
     """
     Parameters
@@ -2152,13 +2154,25 @@ def compute_datasets_chi2(
         dictionary of lists of chi2 per dataset
 
     """
+    # import IPython ; IPython.embed()
+    central_pred_scaled = {}
+    pdf_covmats_scaled = {}
+    t0_covmat = dataset_inputs_covmat_t0_considered
+    if simu_parameters_scales:
+        scaled_bsm = read_bsm_facs / simu_parameters_scales
+    else:
+        scaled_bsm = read_bsm_facs
+    t0_covmats = {}
+    start = 0
+    for dataset in level0_commondata_wc:
+        name = dataset.setname
+        ndata = dataset.ndata
 
-    bsm_facs_df = read_bsm_facs
+        t0_covmats[name] = t0_covmat[start:start+ndata,
+                                    start:start+ndata]
+
+        start += ndata
 
-    means = bsm_facs_df.mean()
-    stds = bsm_facs_df.std()
-    central_pred = {}
-    pdf_covmats = {}
     if dataset_inputs is not None:
         for dataset in dataset_inputs:
             ds = l.check_dataset(
@@ -2173,18 +2187,19 @@ def compute_datasets_chi2(
 
             # Finding PDF/SMEFT fit covmat from replicas
             rep_sm_predictions = predictions(ds, pdf)
-            bsm_factor=dataset_bsm_factor(ds,pdf,read_bsm_facs)
-            rep_bsm_predictions = rep_sm_predictions * bsm_factor
-            replicas = rep_bsm_predictions.iloc[:, 1:]
-            pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True)
-            central_pred[dataset.name] = rep_bsm_predictions.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
-            
+            bsm_factor_scaled = dataset_bsm_factor(ds,pdf,scaled_bsm)
+            rep_bsm_predictions_scaled = rep_sm_predictions * bsm_factor_scaled
+            replicas_scaled = rep_bsm_predictions_scaled.iloc[:, 1:]
+            pdf_covmats_scaled[dataset.name] = np.cov(replicas_scaled, rowvar=True)
+            central_pred_scaled[dataset.name] = rep_bsm_predictions_scaled.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
+
     covmat = groups_covmat # This is the experimental covmat
 
     data = level0_commondata_wc
     contamination_factors = load_datasets_contamination
 
-    chi2_dict = {dataset.setname: [] for dataset in data}
+    chi2_dict_exp = {dataset.setname: [] for dataset in data}
+    chi2_dict_t0 = {dataset.setname: [] for dataset in data}
 
     for dataset in data:
         data_name = dataset.setname
@@ -2201,21 +2216,25 @@ def compute_datasets_chi2(
         covmat_dataset = (
             covmat.xs(data_name, level=1, drop_level=False)
             .T.xs(data_name, level=1, drop_level=False)
-            .values
-        ) + pdf_covmats[data_name]
+            .values) + pdf_covmats_scaled[data_name]
+
+        covmat_dataset_t0 = t0_covmats[data_name] + pdf_covmats_scaled[data_name]
 
-        theory = central_pred[data_name]
+        theory_scaled = central_pred_scaled[data_name]
 
-        diff = (data_values - theory).squeeze()
+        diff_scaled = (data_values - theory_scaled).squeeze()
 
-        if diff.size == 1:
-            chi2 = diff**2 / covmat_dataset[0, 0] / num_data
+        if diff_scaled.size == 1:
+            chi2_exp = diff_scaled**2 / covmat_dataset[0, 0] / num_data
+            chi2_t0 = diff_scaled**2 / covmat_dataset_t0[0, 0] / num_data
         else:
-            chi2 = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) / num_data
+            chi2_exp = (diff_scaled.T @ np.linalg.inv(covmat_dataset) @ diff_scaled) / num_data
+            chi2_t0 = (diff_scaled.T @ np.linalg.inv(covmat_dataset_t0) @ diff_scaled) / num_data
 
-        chi2_dict[data_name].append(chi2)
+        chi2_dict_exp[data_name].append(chi2_exp)
+        chi2_dict_t0[data_name].append(chi2_t0)
 
-    return chi2_dict
+    return chi2_dict_exp, chi2_dict_t0
 
 
 def write_datasets_chi2_dist_csv(
@@ -2273,8 +2292,9 @@ def write_datasets_chi2_csv(
     ndata = {dataset.setname: [dataset.ndata] for dataset in level0_commondata_wc}
 
     df_ndat = pd.DataFrame(ndata)
-    df_chi2 = pd.DataFrame(compute_datasets_chi2)
+    df_chi2_dict_exp = pd.DataFrame(compute_datasets_chi2[0])
+    df_chi2_dict_t0 = pd.DataFrame(compute_datasets_chi2[1])
 
-    chi2 = pd.concat([df_ndat, df_chi2], ignore_index=True)
+    chi2 = pd.concat([df_ndat, df_chi2_dict_exp, df_chi2_dict_t0], ignore_index=True)
 
     chi2.to_csv(f"{pdf}_chi2_dist.csv", index=False)
\ No newline at end of file

From 44979181ca8f4ed870c68aa99832931d2c1635c8 Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Thu, 8 Jan 2026 13:27:00 +0000
Subject: [PATCH 08/14] Removed scaling

---
 validphys2/src/validphys/simunet_analysis.py | 41 +++++++++-----------
 1 file changed, 18 insertions(+), 23 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 29bcef423e..523412e113 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2126,8 +2126,7 @@ def compute_datasets_chi2(
     read_bsm_facs,
     dataset_inputs,
     theoryid,
-    dataset_inputs_covmat_t0_considered,
-    simu_parameters_scales
+    dataset_inputs_covmat_t0_considered
 ):
     """
     Parameters
@@ -2154,14 +2153,10 @@ def compute_datasets_chi2(
         dictionary of lists of chi2 per dataset
 
     """
-    # import IPython ; IPython.embed()
-    central_pred_scaled = {}
-    pdf_covmats_scaled = {}
+    central_pred = {}
+    pdf_covmats = {}
     t0_covmat = dataset_inputs_covmat_t0_considered
-    if simu_parameters_scales:
-        scaled_bsm = read_bsm_facs / simu_parameters_scales
-    else:
-        scaled_bsm = read_bsm_facs
+
     t0_covmats = {}
     start = 0
     for dataset in level0_commondata_wc:
@@ -2187,11 +2182,11 @@ def compute_datasets_chi2(
 
             # Finding PDF/SMEFT fit covmat from replicas
             rep_sm_predictions = predictions(ds, pdf)
-            bsm_factor_scaled = dataset_bsm_factor(ds,pdf,scaled_bsm)
-            rep_bsm_predictions_scaled = rep_sm_predictions * bsm_factor_scaled
-            replicas_scaled = rep_bsm_predictions_scaled.iloc[:, 1:]
-            pdf_covmats_scaled[dataset.name] = np.cov(replicas_scaled, rowvar=True)
-            central_pred_scaled[dataset.name] = rep_bsm_predictions_scaled.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
+            bsm_factor = dataset_bsm_factor(ds,pdf,read_bsm_facs)
+            rep_bsm_predictions = rep_sm_predictions * bsm_factor
+            replicas = rep_bsm_predictions.iloc[:, 1:]
+            pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True)
+            central_pred[dataset.name] = rep_bsm_predictions.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
 
     covmat = groups_covmat # This is the experimental covmat
 
@@ -2216,20 +2211,20 @@ def compute_datasets_chi2(
         covmat_dataset = (
             covmat.xs(data_name, level=1, drop_level=False)
             .T.xs(data_name, level=1, drop_level=False)
-            .values) + pdf_covmats_scaled[data_name]
+            .values) + pdf_covmats[data_name]
 
-        covmat_dataset_t0 = t0_covmats[data_name] + pdf_covmats_scaled[data_name]
+        covmat_dataset_t0 = t0_covmats[data_name] + pdf_covmats[data_name]
 
-        theory_scaled = central_pred_scaled[data_name]
+        theory = central_pred[data_name]
 
-        diff_scaled = (data_values - theory_scaled).squeeze()
+        diff = (data_values - theory).squeeze()
 
-        if diff_scaled.size == 1:
-            chi2_exp = diff_scaled**2 / covmat_dataset[0, 0] / num_data
-            chi2_t0 = diff_scaled**2 / covmat_dataset_t0[0, 0] / num_data
+        if diff.size == 1:
+            chi2_exp = diff**2 / covmat_dataset[0, 0] / num_data
+            chi2_t0 = diff**2 / covmat_dataset_t0[0, 0] / num_data
         else:
-            chi2_exp = (diff_scaled.T @ np.linalg.inv(covmat_dataset) @ diff_scaled) / num_data
-            chi2_t0 = (diff_scaled.T @ np.linalg.inv(covmat_dataset_t0) @ diff_scaled) / num_data
+            chi2_exp = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) / num_data
+            chi2_t0 = (diff.T @ np.linalg.inv(covmat_dataset_t0) @ diff) / num_data
 
         chi2_dict_exp[data_name].append(chi2_exp)
         chi2_dict_t0[data_name].append(chi2_t0)

From c463cdcd3c20424b9d7155407ea80fdf2c3d321f Mon Sep 17 00:00:00 2001
From: Ella Cole <ellagcole02@gmail.com>
Date: Thu, 8 Jan 2026 14:13:30 +0000
Subject: [PATCH 09/14] Added global chi2

---
 validphys2/src/validphys/simunet_analysis.py | 55 ++++++++++++++------
 1 file changed, 40 insertions(+), 15 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 523412e113..f57abcaec8 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2195,7 +2195,8 @@ def compute_datasets_chi2(
 
     chi2_dict_exp = {dataset.setname: [] for dataset in data}
     chi2_dict_t0 = {dataset.setname: [] for dataset in data}
-
+    tot_chi2_exp = 0
+    tot_chi2_t0 = 0
     for dataset in data:
         data_name = dataset.setname
         cont_fac = contamination_factors[data_name]
@@ -2223,13 +2224,20 @@ def compute_datasets_chi2(
             chi2_exp = diff**2 / covmat_dataset[0, 0] / num_data
             chi2_t0 = diff**2 / covmat_dataset_t0[0, 0] / num_data
         else:
-            chi2_exp = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) / num_data
-            chi2_t0 = (diff.T @ np.linalg.inv(covmat_dataset_t0) @ diff) / num_data
+            chi2_exp = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) 
+            chi2_t0 = (diff.T @ np.linalg.inv(covmat_dataset_t0) @ diff)
+            chi2_exp_red = chi2_exp / num_data
+            chi2_t0_red = chi2_t0 / num_data
 
-        chi2_dict_exp[data_name].append(chi2_exp)
-        chi2_dict_t0[data_name].append(chi2_t0)
+            tot_chi2_exp += chi2_exp
+            tot_chi2_t0 += chi2_t0
 
-    return chi2_dict_exp, chi2_dict_t0
+        chi2_dict_exp[data_name].append(chi2_exp_red)
+        chi2_dict_t0[data_name].append(chi2_t0_red)
+    total_ndata = sum([dataset.ndata for dataset in data])
+    tot_chi2_exp_red = tot_chi2_exp / total_ndata
+    tot_chi2_t0_red = tot_chi2_t0 / total_ndata
+    return chi2_dict_exp, chi2_dict_t0, tot_chi2_exp_red, tot_chi2_t0_red
 
 
 def write_datasets_chi2_dist_csv(
@@ -2284,12 +2292,29 @@ def write_datasets_chi2_csv(
     
     """
 
-    ndata = {dataset.setname: [dataset.ndata] for dataset in level0_commondata_wc}
-
-    df_ndat = pd.DataFrame(ndata)
-    df_chi2_dict_exp = pd.DataFrame(compute_datasets_chi2[0])
-    df_chi2_dict_t0 = pd.DataFrame(compute_datasets_chi2[1])
-
-    chi2 = pd.concat([df_ndat, df_chi2_dict_exp, df_chi2_dict_t0], ignore_index=True)
-
-    chi2.to_csv(f"{pdf}_chi2_dist.csv", index=False)
\ No newline at end of file
+    rows = []
+    global_chi2_exp = compute_datasets_chi2[2]
+    global_chi2_t0 = compute_datasets_chi2[3]
+    for dataset in level0_commondata_wc:
+        name = dataset.setname
+        rows.append({
+            "dataset": name,
+            "ndata": dataset.ndata,
+            "chi2_exp": compute_datasets_chi2[0][name],
+            "chi2_t0": compute_datasets_chi2[1][name],
+        })
+
+    df = pd.DataFrame(rows)
+
+    df = pd.concat(
+    [pd.DataFrame([{
+         "dataset": "GLOBAL",
+         "ndata": sum(d.ndata for d in level0_commondata_wc),
+         "chi2_exp": global_chi2_exp,
+         "chi2_t0": global_chi2_t0,
+     }]),
+     df
+    ],
+    ignore_index=True)
+
+    df.to_csv(f"{pdf}_chi2_dist.csv", index=False)

From 7d58b8364f182d881bcb04e8ce5e5cc076bfa2a0 Mon Sep 17 00:00:00 2001
From: Ella Cole <egc41@zprime.damtp.cam.ac.uk>
Date: Fri, 9 Jan 2026 11:01:48 +0000
Subject: [PATCH 10/14] bias=True in covmat

---
 validphys2/src/validphys/simunet_analysis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index f57abcaec8..0d8915da41 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2185,7 +2185,7 @@ def compute_datasets_chi2(
             bsm_factor = dataset_bsm_factor(ds,pdf,read_bsm_facs)
             rep_bsm_predictions = rep_sm_predictions * bsm_factor
             replicas = rep_bsm_predictions.iloc[:, 1:]
-            pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True)
+            pdf_covmats[dataset.name] = np.cov(replicas, rowvar=True,bias=True)
             central_pred[dataset.name] = rep_bsm_predictions.iloc[:, 0].values.squeeze() #Prediction with mean PDF and mean BSM factor
 
     covmat = groups_covmat # This is the experimental covmat

From 08dc2d1a96543723b98ece31418e38a36cea530a Mon Sep 17 00:00:00 2001
From: Ella Cole <egc41@zprime.damtp.cam.ac.uk>
Date: Fri, 9 Jan 2026 13:47:20 +0000
Subject: [PATCH 11/14] fixed bug

---
 validphys2/src/validphys/simunet_analysis.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 0d8915da41..3b4a84aa9f 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2221,8 +2221,10 @@ def compute_datasets_chi2(
         diff = (data_values - theory).squeeze()
 
         if diff.size == 1:
-            chi2_exp = diff**2 / covmat_dataset[0, 0] / num_data
-            chi2_t0 = diff**2 / covmat_dataset_t0[0, 0] / num_data
+            chi2_exp = diff**2 / covmat_dataset[0, 0]
+            chi2_t0 = diff**2 / covmat_dataset_t0[0, 0]
+            chi2_exp_red = chi2_exp / num_data
+            chi2_t0_red = chi2_t0 / num_data
         else:
             chi2_exp = (diff.T @ np.linalg.inv(covmat_dataset) @ diff) 
             chi2_t0 = (diff.T @ np.linalg.inv(covmat_dataset_t0) @ diff)

From c2ddd58587ce0558849da52fbf1d6c0ec2426f8b Mon Sep 17 00:00:00 2001
From: Ella Cole <egc41@zprime.damtp.cam.ac.uk>
Date: Sun, 11 Jan 2026 11:06:54 +0000
Subject: [PATCH 12/14] Fixed load contamination bug

---
 validphys2/src/validphys/simunet_analysis.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 3b4a84aa9f..9d3adc68df 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2040,12 +2040,10 @@ def load_datasets_contamination(
 
                 for op in cont_lin_comb:
                     if op in simu_card[cont_order]:
-                        k_factors += cont_lin_comb[op] * np.array(simu_card[cont_order][op])
-                        
-
+                        k_factors += cont_lin_comb[op] * np.array(simu_card[cont_order][op]) * cont_value
+                                    
             sm_array = np.array(simu_card[cont_order]["SM"])
-            total_factor = 1. + (k_factors*cont_value) / sm_array
-
+            total_factor = 1. + (k_factors / sm_array)
             bsm_dict[dataset.name] = total_factor
 
     return bsm_dict

From dd75d9500463ec7068ed5c03f537f735b443d016 Mon Sep 17 00:00:00 2001
From: Ella Cole <egc41@zprime.damtp.cam.ac.uk>
Date: Sun, 11 Jan 2026 11:39:20 +0000
Subject: [PATCH 13/14] Removed unused variables

---
 validphys2/src/validphys/results.py          | 2 +-
 validphys2/src/validphys/simunet_analysis.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index a847fb09c9..47e4f2e4cb 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -388,7 +388,7 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
         # where ndata is the number of post cut datapoints
         try:
             ndata = len(dataset.load().get_cv())
-        except Exception as e:
+        except Exception:
             ndata = (len(central_predictions(dataset, pdf)))
         nrep = len(pdf)
         return np.ones((ndata, nrep))
diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 9d3adc68df..33e2f3abeb 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2031,7 +2031,6 @@ def load_datasets_contamination(
             k_factors = np.zeros(len(simu_card["SM_fixed"]))
 
             for cont_params in contamination_parameters:
-                cont_name = cont_params["name"]
                 cont_value = cont_params["value"]
                 cont_lin_comb = cont_params["linear_combination"]
 

From 3d5e4c70fb5afa6a029b404f04c1f32ec1d81255 Mon Sep 17 00:00:00 2001
From: Ella Cole <egc41@zprime.damtp.cam.ac.uk>
Date: Thu, 30 Apr 2026 12:25:29 +0100
Subject: [PATCH 14/14] PDF uncertainty option

---
 validphys2/src/validphys/simunet_analysis.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/validphys2/src/validphys/simunet_analysis.py b/validphys2/src/validphys/simunet_analysis.py
index 33e2f3abeb..892417a469 100644
--- a/validphys2/src/validphys/simunet_analysis.py
+++ b/validphys2/src/validphys/simunet_analysis.py
@@ -2123,7 +2123,8 @@ def compute_datasets_chi2(
     read_bsm_facs,
     dataset_inputs,
     theoryid,
-    dataset_inputs_covmat_t0_considered
+    dataset_inputs_covmat_t0_considered,
+    pdf_uncertainty=True
 ):
     """
     Parameters
@@ -2150,6 +2151,7 @@ def compute_datasets_chi2(
         dictionary of lists of chi2 per dataset
 
     """
+    print('PDF_UNCERTAINTY:', pdf_uncertainty)
     central_pred = {}
     pdf_covmats = {}
     t0_covmat = dataset_inputs_covmat_t0_considered
@@ -2209,7 +2211,9 @@ def compute_datasets_chi2(
         covmat_dataset = (
             covmat.xs(data_name, level=1, drop_level=False)
             .T.xs(data_name, level=1, drop_level=False)
-            .values) + pdf_covmats[data_name]
+            .values) 
+        if pdf_uncertainty == True:
+            covmat_dataset += pdf_covmats[data_name]
 
         covmat_dataset_t0 = t0_covmats[data_name] + pdf_covmats[data_name]