Merge branch 'master' into n3lo_dis_fit

n3fit/src/evolven3fit_new/eko_utils.py

@@ -22,7 +22,7 @@
 }
 
 EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
-    "ev_op_iterations": 10,
+    "ev_op_iterations": 30,
     "ev_op_max_order": (1, 0),
     "evolution_method": "iterate-exact",
     "inversion_method": "exact",
@@ -77,6 +77,11 @@ def construct_eko_cards(
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
+    # if Qedref = Qref alphaem is running, otherwise it's fixed
+    if theory["QED"] > 0:
+        if np.isclose(theory["Qref"], theory["Qedref"]):
+            theory_card.couplings.em_running = True
+
     # construct operator card
     q2_grid = utils.generate_q2grid(
         theory["Q0"],

n3fit/src/n3fit/scripts/evolven3fit_new.py

@@ -112,11 +112,8 @@ def main():
     evolven3fit_parser = construct_evolven3fit_parser(subparsers)
 
     args = parser.parse_args()
-    op_card_info = {
-        "ev_op_max_order": 10,
-        "ev_op_iterations": 1,
-        "n_integration_cores": args.n_cores,
-        "backward_inversion": "expanded",
+    op_card_info = { "configs" : {
+        "n_integration_cores": args.n_cores,}
     }
     theory_card_info = {}
     if args.actions == "evolve":

validphys2/src/validphys/checks.py

@@ -43,6 +43,16 @@ def check_pdf_is_montecarlo(ns, **kwargs):
         raise CheckError(f"Error type of PDF {pdf} must be 'replicas' and not {etype}")
 
 
+@make_argcheck
+def check_pdf_is_montecarlo_or_symmhessian(ns, **kwargs):
+    pdf = ns['pdf']
+    etype = pdf.error_type
+    check(
+        etype in {'replicas', 'symhessian'},
+        f"Error type of PDF {pdf} must be either 'replicas' or 'symmhessian' and not {etype}",
+    )
+
+
 @make_check
 def check_know_errors(ns, **kwargs):
     pdf = ns['pdf']

validphys2/src/validphys/config.py

@@ -1593,7 +1593,7 @@ def produce_group_dataset_inputs_by_metadata(
             if processed_metadata_group == "custom_group":
                 group_name = str(dsinput.custom_group)
             # special case of ALL, grouping everything together
-            if processed_metadata_group == "ALL":
+            elif processed_metadata_group == "ALL":
                 group_name = processed_metadata_group
             # otherwise try and take the key from the metadata.
             else:

validphys2/src/validphys/covmats.py

@@ -15,7 +15,7 @@
     check_data_cuts_match_theorycovmat,
     check_dataset_cuts_match_theorycovmat,
     check_norm_threshold,
-    check_pdf_is_montecarlo,
+    check_pdf_is_montecarlo_or_symmhessian,
     check_speclabels_different,
 )
 from validphys.commondata import loaded_commondata_with_cuts
@@ -677,11 +677,15 @@ def groups_corrmat(groups_covmat):
     return mat
 
 
-@check_pdf_is_montecarlo
+@check_pdf_is_montecarlo_or_symmhessian
 def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
     """For a given `dataset`, returns the sum of the covariance matrix given by
-    `covmat_t0_considered` and the PDF error: a covariance matrix estimated from the
-    replica theory predictions from a given monte carlo `pdf`
+    `covmat_t0_considered` and the PDF error:
+    - If the PDF error_type is 'replicas', a covariance matrix is estimated from
+      the replica theory predictions
+    - If the PDF error_type is 'symmhessian', a covariance matrix is estimated using
+      formulas from (mc2hessian) https://arxiv.org/pdf/1505.06736.pdf
+
 
     Parameters
     ----------
@@ -716,7 +720,22 @@ def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
     True
     """
     th = ThPredictionsResult.from_convolution(pdf, dataset)
-    pdf_cov = np.cov(th.error_members, rowvar=True)
+
+    if pdf.error_type == 'replicas':
+        pdf_cov = np.cov(th.error_members, rowvar=True)
+
+    elif pdf.error_type == 'symmhessian':
+        rescale_fac = pdf._rescale_factor()
+        hessian_eigenvectors = th.error_members
+        central_predictions = th.central_value
+
+        # need to subtract the central set which is not the same as the average of the
+        # Hessian eigenvectors.
+        X = hessian_eigenvectors - central_predictions.reshape((central_predictions.shape[0], 1))
+        # need to rescale the Hessian eigenvectors in case the eigenvector confidence interval is not 68%
+        X = X / rescale_fac
+        pdf_cov = X @ X.T
+
     return pdf_cov + covmat_t0_considered
 
 

validphys2/src/validphys/photon/compute.py

@@ -17,7 +17,6 @@
 
 log = logging.getLogger(__name__)
 
-Q_IN = 100
 FIATLUX_DEFAULT = {
     "apfel": False,
     "eps_base": 1e-5,  # precision on final integration of double integral.
@@ -71,6 +70,9 @@ def __init__(self, theoryid, lux_params, replicas):
         path_to_F2 = theoryid.path / "fastkernel/fiatlux_dis_F2.pineappl.lz4"
         path_to_FL = theoryid.path / "fastkernel/fiatlux_dis_FL.pineappl.lz4"
         self.path_to_eko_photon = theoryid.path / "eko_photon.tar"
+        with EKO.read(self.path_to_eko_photon) as eko:
+            self.q_in = np.sqrt(eko.mu20)
+
 
         # set fiatlux
         self.lux = {}
@@ -129,7 +131,7 @@ def compute_photon_array(self, replica):
         """
         # Compute photon PDF
         log.info(f"Computing photon")
-        photon_qin = np.array([self.lux[replica].EvaluatePhoton(x, Q_IN**2).total for x in XGRID])
+        photon_qin = np.array([self.lux[replica].EvaluatePhoton(x, self.q_in**2).total for x in XGRID])
         photon_qin += self.generate_errors(replica)
         # fiatlux computes x * gamma(x)
         photon_qin /= XGRID
@@ -141,21 +143,20 @@ def compute_photon_array(self, replica):
             # it has to be done inside vp-setupfit
 
             # construct PDFs
-            pdfs_init = np.zeros((len(eko.rotations.inputpids), len(XGRID)))
-            for j, pid in enumerate(eko.rotations.inputpids):
+            pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID)))
+            for j, pid in enumerate(eko.bases.inputpids):
                 if pid == 22:
                     pdfs_init[j] = photon_qin
                     ph_id = j
                 else:
                     if pid not in self.luxpdfset.flavors:
                         continue
                     pdfs_init[j] = np.array(
-                        [self.luxpdfset.xfxQ(x, Q_IN, replica, pid) / x for x in XGRID]
+                        [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID]
                     )
 
             # Apply EKO to PDFs
-            q2 = eko.mu2grid[0]
-            with eko.operator(q2) as elem:
+            for _, elem in eko.items():
                 pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init)
 
         photon_Q0 = pdfs_final[ph_id]
@@ -188,7 +189,7 @@ def error_matrix(self):
         if not self.additional_errors:
             return None
         extra_set = self.additional_errors.load()
-        qs = [Q_IN] * len(XGRID)
+        qs = [self.q_in] * len(XGRID)
         res_central = np.array(extra_set.central_member.xfxQ(22, XGRID, qs))
         res = []
         for idx_member in range(101, 107 + 1):

validphys2/src/validphys/tests/photon/test_compute.py

@@ -8,6 +8,7 @@
 from validphys.core import PDF as PDFset
 
 from ..conftest import PDF
+from eko.io import EKO
 
 
 class FakeTheory:
@@ -68,6 +69,17 @@ def PlugStructureFunctions(self, f2, fl, f2lo):
     def EvaluatePhoton(self, x, q):
         return self.res
 
+class FakeEKO:
+    def __init__(self, path):
+        self.path = path
+        self.mu20 = 100**2
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, exc_type: type, _exc_value, _traceback):
+        pass
+
 
 class FakeStructureFunction:
     def __init__(self, path, pdfs):
@@ -95,6 +107,7 @@ def test_parameters_init(monkeypatch):
     monkeypatch.setattr(structure_functions, "F2LO", FakeF2LO)
     monkeypatch.setattr(fiatlux, "FiatLux", FakeFiatlux)
     monkeypatch.setattr(Photon, "compute_photon_array", lambda *args: np.zeros(196))
+    monkeypatch.setattr(EKO, "read", FakeEKO)
 
     photon = Photon(FakeTheory(), fiatlux_runcard, [1, 2, 3])
     alpha = Alpha(FakeTheory().get_description())