fix for the qed part

Author: scarlehoff

validphys2/src/validphys/photon/compute.py

@@ -8,6 +8,7 @@
 from scipy.interpolate import interp1d
 import yaml
 
+from eko import basis_rotation
 from eko.io import EKO
 from n3fit.io.writer import XGRID
 from validphys.n3fit_data import replica_luxseed
@@ -165,26 +166,27 @@ def compute_photon_array(self, replica):
         # TODO : the different x points could be even computed in parallel
 
         # Load eko and reshape it
-        with EKO.read(self.path_to_eko_photon) as eko:
+        with EKO.read(self.path_to_eko_photon) as eko_photon:
             # TODO : if the eko has not the correct grid we have to reshape it
             # it has to be done inside vp-setupfit
 
-            # construct PDFs
-            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:
+            # NB: the eko should contain a single operator
+            for _, elem in eko_photon.items():
+                eko_op = elem.operator
+
+                pdfs_init = np.zeros_like(eko_op[0, 0])
+                for j, pid in enumerate(basis_rotation.flavor_basis_pids):
+                    if pid == 22:
+                        pdfs_init[j] = photon_qin
+                        ph_id = j
+                    elif pid not in self.luxpdfset.flavors:
                         continue
-                    pdfs_init[j] = np.array(
-                        [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID]
-                    )
+                    else:
+                        pdfs_init[j] = np.array(
+                            [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID]
+                        )
 
-            # Apply EKO to PDFs
-            for _, elem in eko.items():
-                pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init)
+                pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
 
         photon_Q0 = pdfs_final[ph_id]
 

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

@@ -100,9 +100,7 @@ def test_couplings_exa():
             dict(
                 alphas=theory["alphas"],
                 alphaem=theory["alphaqed"],
-                scale=theory["Qref"],
-                num_flavs_ref=None,
-                max_num_flavs=theory["MaxNfPdf"],
+                ref=(theory["Qref"], theory["MaxNfPdf"]),
                 em_running=True,
             )
         )

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

@@ -4,6 +4,7 @@
 import numpy as np
 import yaml
 
+from eko import basis_rotation
 from eko.io import EKO
 from n3fit.io.writer import XGRID
 from validphys.api import API
@@ -94,8 +95,8 @@ def test_photon():
             photon_fiatlux_qin = np.array([lux.EvaluatePhoton(x, eko.mu20).total for x in XGRID])
             photon_fiatlux_qin /= XGRID
             # construct PDFs
-            pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID)))
-            for j, pid in enumerate(eko.bases.inputpids):
+            pdfs_init = np.zeros((len(basis_rotation.flavor_basis_pids), len(XGRID)))
+            for j, pid in enumerate(basis_rotation.flavor_basis_pids):
                 if pid == 22:
                     pdfs_init[j] = photon_fiatlux_qin
                     ph_id = j