From e135fd62fbb83734a6e8723d9a3625310cee76f8 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 Oct 2025 15:16:42 +0000
Subject: [PATCH 01/30] Making the photon a downloadable resource

---
 n3fit/src/n3fit/checks.py                  |  10 ++
 n3fit/src/n3fit/performfit.py              |   2 +-
 validphys2/src/validphys/loader.py         |  57 ++++++
 validphys2/src/validphys/photon/compute.py | 191 +++++++++++----------
 validphys2/src/validphys/scripts/vp_get.py |   7 +-
 5 files changed, 174 insertions(+), 93 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index ab74c14b15..fdc2ec0041 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -459,6 +459,16 @@ def check_multireplica_qed(replicas, fiatlux):
     if fiatlux is not None:
         if len(replicas) > 1:
             raise CheckError("At the moment, running a multireplica QED fits is not allowed.")
+        
+@make_argcheck
+def check_photonQED_exists(theoryid, fiatlux):
+    """Check that the Photon QED set for this theoryid and luxset exists"""
+    if fiatlux is not None:
+      luxset = fiatlux['luxset']
+      try:
+          _ = FallbackLoader().check_photonQED(theoryid, luxset)
+      except FileNotFoundError:
+          raise CheckError(f"Photon QED set for TheoryID {theoryid.id} and luxset {luxset} not found.")
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py
index e0fd1af9b5..312885de4c 100644
--- a/n3fit/src/n3fit/performfit.py
+++ b/n3fit/src/n3fit/performfit.py
@@ -13,7 +13,7 @@
 
 # Action to be called by validphys
 # All information defining the NN should come here in the "parameters" dict
-@n3fit.checks.check_multireplica_qed
+@n3fit.checks.check_photonQED_exists
 @n3fit.checks.check_polarized_configs
 def performfit(
     *,
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index ed10582319..2462e2be40 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -80,6 +80,10 @@ class EkoNotFound(LoadFailedError):
     pass
 
 
+class PhotonQEDNotFound(LoadFailedError):
+    pass
+
+
 class TheoryMetadataNotFound(LoadFailedError):
     pass
 
@@ -143,16 +147,19 @@ def __init__(self, profile=None):
         theories_path = pathlib.Path(profile["theories_path"])
         resultspath = pathlib.Path(profile["results_path"])
         ekos_path = pathlib.Path(profile["ekos_path"])
+        photons_qed = pathlib.Path(profile["photons_qed_path"])
 
         # Create the theories and results paths if they don't exist already
         theories_path.mkdir(exist_ok=True, parents=True)
         ekos_path.mkdir(exist_ok=True, parents=True)
         resultspath.mkdir(exist_ok=True, parents=True)
+        photons_qed.mkdir(exist_ok=True, parents=True)
 
         # And save them up
         self.commondata_folders = tuple(datapaths)
         self._theories_path = theories_path
         self._ekos_path = ekos_path
+        self._photons_qed_path = photons_qed
         self.resultspath = resultspath
         self._extremely_old_fits = set()
         self.nnprofile = profile
@@ -209,6 +216,14 @@ def available_ekos(self):
         return {
             eko_path.parent.name.split("_")[1] for eko_path in self._theories_path.glob("*/eko.tar")
         }
+    
+    @property
+    @functools.lru_cache
+    def available_photons_qed(self):
+        """Return a string token for each of the available theories"""
+        return {
+            eko_path.parent.name.split("_")[1] for eko_path in self._theories_path.glob("*/eko.npz")
+        }
 
     @property
     @functools.lru_cache
@@ -305,6 +320,14 @@ def check_eko(self, theoryID):
         if not eko_path.exists():
             raise EkoNotFound(f"Could not find eko {eko_path} in theory: {theoryID}")
         return eko_path
+    
+    @functools.lru_cache
+    def check_photonQED(self, theoryID, luxset):
+        """Check the Photon QED set exists and return the path to it"""
+        photon_qed_path = self._photons_qed_path / f"photon_qed_{theoryID.id}_{luxset}.tar"
+        if not photon_qed_path.exists():
+            raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {theoryID}")
+        return photon_qed_path
 
     @property
     def theorydb_folder(self):
@@ -815,6 +838,16 @@ def eko_index(self):
     @_key_or_loader_error
     def eko_urls(self):
         return self.nnprofile['eko_urls']
+    
+    @property
+    @_key_or_loader_error
+    def photon_qed_index(self):
+        return self.nnprofile['photon_qed_index']
+    
+    @property
+    @_key_or_loader_error
+    def photon_qed_urls(self):
+        return self.nnprofile['photon_qed_urls']
 
     @property
     @_key_or_loader_error
@@ -838,6 +871,7 @@ def lhapdf_urls(self):
 
     def _remote_files_from_url(self, url, index, thing='files'):
         index_url = url + index
+        import ipdb; ipdb.set_trace()
         try:
             resp = requests.get(index_url)
             resp.raise_for_status()
@@ -885,6 +919,13 @@ def remote_ekos(self):
         token = 'eko_'
         rt = self.remote_files(self.eko_urls, self.eko_index, thing="ekos")
         return {k[len(token) :]: v for k, v in rt.items()}
+    
+    @property
+    @functools.lru_cache
+    def remote_photons_qed(self):
+        token = 'photon_qed_'
+        rt = self.remote_files(self.photon_qed_urls, self.photon_qed_index, thing="photons_qed")
+        return {k[len(token) :]: v for k, v in rt.items()}
 
     @property
     @functools.lru_cache
@@ -919,6 +960,10 @@ def downloadable_theories(self):
     @property
     def downloadable_ekos(self):
         return list(self.remote_ekos)
+    
+    @property
+    def downloadable_photonsQED(self):
+        return list(self.remote_photons_qed)
 
     @property
     def lhapdf_pdfs(self):
@@ -1102,6 +1147,18 @@ def download_eko(self, thid):
         target_path = self._ekos_path / f"eko_{int(thid)}.tar"
         download_file(remote[thid], target_path)
 
+    def download_photonQED(self, thid, luxset: str):
+        """Download the Photon set for a given theory ID"""
+        # thid = str(thid)
+        # remote = self.remote_photons_qed
+        # key = f"{thid}_{luxset}"
+        # if key not in remote:
+        #     raise PhotonQEDNotFound(f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server")
+        # # Check that we have the theory we need
+        # target_path = self._photon_qed_path / f"photon_qed_{int(thid)}_{luxset}.tar"
+        # download_file(remote[key], target_path)
+        log.warning("Downloading Photon QED sets is not implemented yet.")
+    
     def download_vp_output_file(self, filename, **kwargs):
         try:
             root_url = self.nnprofile['reports_root_url']
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index a4a3038d8b..5388b7829e 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -2,6 +2,7 @@
 
 import logging
 import tempfile
+from concurrent.futures import ThreadPoolExecutor
 
 import numpy as np
 from scipy.integrate import trapezoid
@@ -12,11 +13,13 @@
 from eko.io import EKO
 from n3fit.io.writer import XGRID
 from validphys.n3fit_data import replica_luxseed
+from validphys.loader import Loader, PhotonQEDNotFound
 
 from . import structure_functions as sf
 from .alpha import Alpha
 
 log = logging.getLogger(__name__)
+loader = Loader()
 
 # not the complete fiatlux runcard since some parameters are set in the code
 FIATLUX_DEFAULT = {
@@ -49,18 +52,30 @@
 
 
 class Photon:
-    """Photon class computing the photon array with the LuxQED approach."""
-
+    """Photon class computing the photon array with the LuxQED approach.
+
+    Parameters
+    ----------
+    theoryid : validphys.core.TheoryIDSpec
+        TheoryIDSpec object describing the theory to be used as
+        specified in the runcard.
+    lux_params : dict
+        Dictionary containing the LuxQED parameters as specified
+        in the runcard.
+    replica_list: list[int], optional
+        List of replica ids to be computed. If None, all replicas
+        will be computed based on the luxqed pdf set.
+    """
     def __init__(self, theoryid, lux_params, replicas):
         self.theoryid = theoryid
         self.lux_params = lux_params
+        self.replicas = replicas
 
-        theory = theoryid.get_description()
         fiatlux_runcard = FIATLUX_DEFAULT
         # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
         # This may be changed in the future in favor of a bool em_running in the runcard
         fiatlux_runcard["qed_running"] = True
-        fiatlux_runcard["mproton"] = float(theory["MP"])
+        fiatlux_runcard["mproton"] = float(theoryid.get_description()["MP"])
 
         # precision on final integration of double integral
         if "eps_base" in lux_params:
@@ -70,37 +85,46 @@ def __init__(self, theoryid, lux_params, replicas):
             fiatlux_runcard["eps_base"] = 1e-5
             log.info(f"Using default value for fiatlux parameter eps_base")
 
-        self.replicas = replicas
-
-        # structure functions
+        self.fiatlux_runcard = fiatlux_runcard
+        # Metadata for the photon ste
         self.luxpdfset = lux_params["luxset"].load()
         self.additional_errors = lux_params["additional_errors"]
         self.luxseed = lux_params["luxseed"]
+        self.luxpdfset_members = self.luxpdfset.n_members - 1 # Remove replica 0
+
+        try:
+          self.load_photon()
+        except PhotonQEDNotFound:
+          log.info(f"Photon set for theory ID {self.theoryid.id} and luxset {self.luxpdfset._name} not found. Computing it now...")
+          self.compute_photon_set()
+
+
+    def compute_photon_set(self):
+        """Compute the photon set for the desired replicas."""
+
+
+        # load fiatlux
+        try:
+          import fiatlux
+        except ModuleNotFoundError as e:
+          log.error("fiatlux not found, please install fiatlux")
+          raise ModuleNotFoundError("Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`") from e
+
+        theory = self.theoryid.get_description()
 
         if theory["PTO"] > 0:
-            path_to_F2 = theoryid.path / "fastkernel/FIATLUX_DIS_F2.pineappl.lz4"
-            path_to_FL = theoryid.path / "fastkernel/FIATLUX_DIS_FL.pineappl.lz4"
+            path_to_F2 = self.theoryid.path / "fastkernel/FIATLUX_DIS_F2.pineappl.lz4"
+            path_to_FL = self.theoryid.path / "fastkernel/FIATLUX_DIS_FL.pineappl.lz4"
 
-        self.path_to_eko_photon = theoryid.path / "eko_photon.tar"
+        self.path_to_eko_photon = self.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 = {}
-
         mb_thr = theory["kbThr"] * theory["mb"]
         mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
-
-        self.interpolator = []
-        self.integral = []
-
-        try:
-            import fiatlux
-        except ModuleNotFoundError as e:
-            log.error("fiatlux not found, please install fiatlux")
-            raise ModuleNotFoundError(
-                "Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`"
-            ) from e
+        interpolator = []
+        integral = []
 
         for replica in self.replicas:
             # As input replica for the photon computation we take the MOD of the luxset_members to
@@ -117,85 +141,69 @@ def __init__(self, theoryid, lux_params, replicas):
                     log.error(
                         "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
                     )
-                fiatlux_runcard["q2_max"] = float(f2.q2_max)
+                self.fiatlux_runcard["q2_max"] = float(f2.q2_max)
             else:
                 f2 = f2lo
                 fl = sf.FLLO()
                 # using a default value for q2_max
-                fiatlux_runcard["q2_max"] = 1e8
-
-            alpha = Alpha(theory, fiatlux_runcard["q2_max"])
+                self.fiatlux_runcard["q2_max"] = 1e8
 
+            alpha = Alpha(theory, self.fiatlux_runcard["q2_max"])
             with tempfile.NamedTemporaryFile(mode="w") as tmp:
-                yaml.dump(fiatlux_runcard, tmp)
-                self.lux[replica] = fiatlux.FiatLux(tmp.name)
+                yaml.dump(self.fiatlux_runcard, tmp)
+                lux = fiatlux.FiatLux(tmp.name)
+
             # we have a dict but fiatlux wants a yaml file
             # TODO : once that fiatlux will allow dictionaries
             # pass directly fiatlux_runcard
+            lux.PlugAlphaQED(alpha.alpha_em, alpha.qref)
+            lux.InsertInelasticSplitQ([mb_thr, mt_thr])
+            lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
 
-            self.lux[replica].PlugAlphaQED(alpha.alpha_em, alpha.qref)
-            self.lux[replica].InsertInelasticSplitQ([mb_thr, mt_thr])
-            self.lux[replica].PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
+            # Evaluate photon for every point in the grid xgrid
+            def evaluate_at_x(x):
+              return lux.EvaluatePhoton(x, self.q_in**2).total
+            with ThreadPoolExecutor() as executor:
+                photon_qin = np.array(list(executor.map(evaluate_at_x, XGRID)))
 
-            photon_array = self.compute_photon_array(replica, photonreplica)
-            self.interpolator.append(
-                interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic")
-            )
-            self.integral.append(trapezoid(photon_array, XGRID))
+            photon_qin += self.generate_errors(replica)
 
-        self.integral = np.stack(self.integral, axis=-1)
+            # fiatlux computes x * gamma(x)
+            photon_qin /= XGRID
 
-    def compute_photon_array(self, replica, photonreplica):
-        r"""
-        Compute the photon PDF for every point in the grid xgrid.
+            # Load eko and reshape it
+            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
 
-        Parameters
-        ----------
-        replica: int
-            replica id
+                # NB: the eko should contain a single operator
+                for _, elem in eko_photon.items():
+                    eko_op = elem.operator
 
-        Returns
-        -------
-        compute_photon_array: numpy.array
-            photon PDF at the fitting scale Q0
-        """
-        # Compute photon PDF
-        log.info(f"Computing photon")
-        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
-        # 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_photon:
-            # TODO : if the eko has not the correct grid we have to reshape it
-            # it has to be done inside vp-setupfit
-
-            # 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
-                    else:
-                        pdfs_init[j] = np.array(
-                            [self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x for x in XGRID]
-                        )
-
-                pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
-
-        photon_Q0 = pdfs_final[ph_id]
-
-        # we want x * gamma(x)
-        return XGRID * photon_Q0
+                    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
+                        else:
+                            pdfs_init[j] = np.array(
+                                [self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x for x in XGRID]
+                            )
+
+                    pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
+
+            photon_Q0 = pdfs_final[ph_id]
+            photon_array = XGRID * photon_Q0
+            interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
+            integral.append(trapezoid(photon_array, XGRID))
+
+        integral = np.stack(self.integral, axis=-1)
+
+
+        self.integral = integral
+        self.interpolator = interpolator
 
     def __call__(self, xgrid):
         """
@@ -247,3 +255,12 @@ def generate_errors(self, replica):
         u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False)
         errors = u @ (s * rng.normal(size=7))
         return errors
+
+    def load_photon(self):
+      """Load the photon resource using the Loader class."""
+      path_to_photon = loader.check_photonQED(self.theoryid, self.luxpdfset._name)
+      log.info(f"Loading photon QED set from {path_to_photon}")
+
+      log.warning("Loading photon QED set is not yet implemented.")
+      exit(1)
+      return
\ No newline at end of file
diff --git a/validphys2/src/validphys/scripts/vp_get.py b/validphys2/src/validphys/scripts/vp_get.py
index 51f88ed4c9..70f06a87b0 100644
--- a/validphys2/src/validphys/scripts/vp_get.py
+++ b/validphys2/src/validphys/scripts/vp_get.py
@@ -21,8 +21,6 @@
 from reportengine import colors
 from validphys.loader import FallbackLoader as Loader, LoadFailedError
 
-
-
 log = logging.getLogger()
 log.setLevel(logging.INFO)
 log.addHandler(colors.ColorHandler())
@@ -51,12 +49,11 @@ def main():
         sys.exit(1)
     p.add_argument('resource_type', help="Type of the resource to be obtained. "
                    "See --list for a list of resource types.")
-    p.add_argument('resource_name', help="Identifier of the resource.")
+    p.add_argument('resource_name', help="Identifier of the resource.", nargs='+')
     p.add_argument('--list', action=ListAction, loader=l, nargs=0,
                    help="List available resources and exit.")
     args = p.parse_args()
 
-
     tp = args.resource_type
     name = args.resource_name
 
@@ -66,7 +63,7 @@ def main():
         sys.exit(f"No such resource {tp}")
 
     try:
-        res = f(name)
+        res = f(*name)
     except LoadFailedError as e:
         print(ErrorWithAlternatives(f"Could not find resource ({tp}): '{name}'.", name))
         sys.exit("Failed to download resource.")

From 83d66620cb6c7770b184c5650354ee3818c90fed Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 Oct 2025 15:17:43 +0000
Subject: [PATCH 02/30] Add photon resource to nnprofile

---
 validphys2/src/validphys/nnprofile_default.yaml | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/validphys2/src/validphys/nnprofile_default.yaml b/validphys2/src/validphys/nnprofile_default.yaml
index 85e4a351f4..c2d7e77965 100644
--- a/validphys2/src/validphys/nnprofile_default.yaml
+++ b/validphys2/src/validphys/nnprofile_default.yaml
@@ -30,6 +30,7 @@ theories_path: theories
 hyperscan_path: hyperscan
 validphys_cache_path: vp-cache
 ekos_path: ekos
+photons_qed_path: photons_qed
 
 # It is possible to add extra folders in which to find data by filling up the data_path variable
 # the default path within the nnpdf_data package will always be added albeit with lower priority.
@@ -61,6 +62,11 @@ eko_urls:
 
 eko_index: 'ekodata.json'
 
+photon_qed_urls:
+    - 'https://nnpdf.nikhef.nl/nnpdf/photons_qed/'
+
+photon_qed_index: 'photondata.json'
+
 lhapdf_urls:
     - 'http://lhapdfsets.web.cern.ch/lhapdfsets/current/'
 nnpdf_pdfs_urls:

From dbf92e304c7e8183105077839818f7755450082e Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 Oct 2025 17:51:15 +0000
Subject: [PATCH 03/30] Utility script to run FiatLux

---
 pyproject.toml                             |  2 +
 validphys2/src/validphys/loader.py         |  2 +-
 validphys2/src/validphys/photon/compute.py | 46 +++++++++++++++-------
 3 files changed, 35 insertions(+), 15 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 5bdc162905..acb0f039da 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -55,6 +55,8 @@ vp-list = "validphys.scripts.vp_list:main"
 vp-nextfitruncard = "validphys.scripts.vp_nextfitruncard:main"
 vp-hyperoptplot = "validphys.scripts.vp_hyperoptplot:main"
 vp-deltachi2 = "validphys.scripts.vp_deltachi2:main"
+# FiatLux scripts
+fiatlux-run = "n3fit.scripts.fiatlux_exec:main"
 
 [tool.poetry.dependencies]
 # Generic dependencies (i.e., validphys)
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index 2462e2be40..6b30daf4c4 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -324,7 +324,7 @@ def check_eko(self, theoryID):
     @functools.lru_cache
     def check_photonQED(self, theoryID, luxset):
         """Check the Photon QED set exists and return the path to it"""
-        photon_qed_path = self._photons_qed_path / f"photon_qed_{theoryID.id}_{luxset}.tar"
+        photon_qed_path = self._photons_qed_path / f"photon_qed_{theoryID.id}_{luxset}"
         if not photon_qed_path.exists():
             raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {theoryID}")
         return photon_qed_path
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 5388b7829e..f7da3a4482 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -53,7 +53,7 @@
 
 class Photon:
     """Photon class computing the photon array with the LuxQED approach.
-
+    
     Parameters
     ----------
     theoryid : validphys.core.TheoryIDSpec
@@ -66,10 +66,11 @@ class Photon:
         List of replica ids to be computed. If None, all replicas
         will be computed based on the luxqed pdf set.
     """
-    def __init__(self, theoryid, lux_params, replicas):
+    def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False):
         self.theoryid = theoryid
         self.lux_params = lux_params
         self.replicas = replicas
+        self.save_to_disk = save_to_disk
 
         fiatlux_runcard = FIATLUX_DEFAULT
         # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
@@ -98,18 +99,15 @@ def __init__(self, theoryid, lux_params, replicas):
           log.info(f"Photon set for theory ID {self.theoryid.id} and luxset {self.luxpdfset._name} not found. Computing it now...")
           self.compute_photon_set()
 
-
     def compute_photon_set(self):
         """Compute the photon set for the desired replicas."""
-
-
         # load fiatlux
         try:
           import fiatlux
         except ModuleNotFoundError as e:
           log.error("fiatlux not found, please install fiatlux")
           raise ModuleNotFoundError("Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`") from e
-
+        
         theory = self.theoryid.get_description()
 
         if theory["PTO"] > 0:
@@ -129,8 +127,7 @@ def compute_photon_set(self):
         for replica in self.replicas:
             # As input replica for the photon computation we take the MOD of the luxset_members to
             # avoid failing due to limited number of replicas in the luxset
-            luxset_members = self.luxpdfset.n_members - 1  # - 1 because rep0 is included
-            photonreplica = (replica % luxset_members) or luxset_members
+            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
 
             f2lo = sf.F2LO(self.luxpdfset.members[photonreplica], theory)
 
@@ -163,10 +160,11 @@ def compute_photon_set(self):
             # Evaluate photon for every point in the grid xgrid
             def evaluate_at_x(x):
               return lux.EvaluatePhoton(x, self.q_in**2).total
+            
             with ThreadPoolExecutor() as executor:
                 photon_qin = np.array(list(executor.map(evaluate_at_x, XGRID)))
 
-            photon_qin += self.generate_errors(replica)
+            # photon_qin += self.generate_errors(replica)
 
             # fiatlux computes x * gamma(x)
             photon_qin /= XGRID
@@ -175,7 +173,7 @@ def evaluate_at_x(x):
             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
-
+                
                 # NB: the eko should contain a single operator
                 for _, elem in eko_photon.items():
                     eko_op = elem.operator
@@ -196,10 +194,17 @@ def evaluate_at_x(x):
 
             photon_Q0 = pdfs_final[ph_id]
             photon_array = XGRID * photon_Q0
+
+            if self.save_to_disk:
+                path_to_photon = loader._photons_qed_path / f"photon_qed_{self.theoryid.id}_{self.luxpdfset._name}"
+                path_to_photon.mkdir(parents=True, exist_ok=True)
+                np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
+                log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
+
             interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
             integral.append(trapezoid(photon_array, XGRID))
 
-        integral = np.stack(self.integral, axis=-1)
+        integral = np.stack(integral, axis=-1)
 
 
         self.integral = integral
@@ -255,12 +260,25 @@ def generate_errors(self, replica):
         u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False)
         errors = u @ (s * rng.normal(size=7))
         return errors
-
+    
     def load_photon(self):
       """Load the photon resource using the Loader class."""
       path_to_photon = loader.check_photonQED(self.theoryid, self.luxpdfset._name)
       log.info(f"Loading photon QED set from {path_to_photon}")
 
-      log.warning("Loading photon QED set is not yet implemented.")
-      exit(1)
+      interpolator = []
+      integral = []
+
+      # Load the needed replicas
+      for replica in self.replicas:
+            # As input replica for the photon computation we take the MOD of the luxset_members to
+            # avoid failing due to limited number of replicas in the luxset
+            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+
+            photon_array = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_array"]
+            interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
+            integral.append(trapezoid(photon_array, XGRID))
+      
+      self.interpolator = interpolator
+      self.integral = np.stack(integral, axis=-1)
       return
\ No newline at end of file

From 5fe121302d7a35a29783389404ebd92973c64b0d Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 Oct 2025 18:00:00 +0000
Subject: [PATCH 04/30] Remove debugger

---
 n3fit/src/n3fit/scripts/fiatlux_exec.py | 46 +++++++++++++++++++++++++
 validphys2/src/validphys/loader.py      |  1 -
 2 files changed, 46 insertions(+), 1 deletion(-)
 create mode 100755 n3fit/src/n3fit/scripts/fiatlux_exec.py

diff --git a/n3fit/src/n3fit/scripts/fiatlux_exec.py b/n3fit/src/n3fit/scripts/fiatlux_exec.py
new file mode 100755
index 0000000000..fbedb03e09
--- /dev/null
+++ b/n3fit/src/n3fit/scripts/fiatlux_exec.py
@@ -0,0 +1,46 @@
+import argparse
+import pathlib
+import logging
+
+from validphys.loader import Loader
+from validphys.utils import yaml_safe
+from validphys.photon.compute import Photon
+
+def check_positive(value):
+    ivalue = int(value)
+    if ivalue <= 0:
+        raise argparse.ArgumentTypeError("%s is an invalid positive int value." % value)
+    return ivalue
+
+log = logging.getLogger(__name__)
+loader = Loader()
+
+
+
+def main():
+    parser = argparse.ArgumentParser(description="FitLux-exec - compute the photon PDF using FiatLux")
+    parser.add_argument("config_yml", help="Path to the configuration file", type=pathlib.Path)
+    parser.add_argument("replica", help="MC replica number", type=check_positive)
+    parser.add_argument(
+        "-r",
+        "--replica_range",
+        help="End of the range of replicas to compute",
+        type=check_positive,
+    )
+    args = parser.parse_args()
+    config_path = args.config_yml.absolute()
+    replica = args.replica
+    if args.replica_range:
+        replicas = list(range(replica, args.replica_range + 1))
+    else:
+        replicas = [replica]
+
+    runcard = yaml_safe.load(config_path.read_text(encoding="UTF-8"))
+    theoryID = loader.check_theoryID(runcard["theory"]["theoryid"])
+    fiatlux_params = runcard.get("fiatlux", None)
+    fiatlux_params['luxset'] = loader.check_pdf(fiatlux_params['luxset'])
+
+    Photon(theoryID, fiatlux_params, replicas=replicas, save_to_disk=True)
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index 6b30daf4c4..a7bd1bff6b 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -871,7 +871,6 @@ def lhapdf_urls(self):
 
     def _remote_files_from_url(self, url, index, thing='files'):
         index_url = url + index
-        import ipdb; ipdb.set_trace()
         try:
             resp = requests.get(index_url)
             resp.raise_for_status()

From 2fd90aa331842ccd0f3a59956988dc33b3ebe4ff Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 Oct 2025 20:18:09 +0000
Subject: [PATCH 05/30] Allow force computation

---
 n3fit/src/n3fit/scripts/fiatlux_exec.py    |  4 ++--
 validphys2/src/validphys/photon/compute.py | 19 ++++++++++++-------
 2 files changed, 14 insertions(+), 9 deletions(-)

diff --git a/n3fit/src/n3fit/scripts/fiatlux_exec.py b/n3fit/src/n3fit/scripts/fiatlux_exec.py
index fbedb03e09..735bcf3957 100755
--- a/n3fit/src/n3fit/scripts/fiatlux_exec.py
+++ b/n3fit/src/n3fit/scripts/fiatlux_exec.py
@@ -2,7 +2,7 @@
 import pathlib
 import logging
 
-from validphys.loader import Loader
+from validphys.loader import FallbackLoader as Loader
 from validphys.utils import yaml_safe
 from validphys.photon.compute import Photon
 
@@ -40,7 +40,7 @@ def main():
     fiatlux_params = runcard.get("fiatlux", None)
     fiatlux_params['luxset'] = loader.check_pdf(fiatlux_params['luxset'])
 
-    Photon(theoryID, fiatlux_params, replicas=replicas, save_to_disk=True)
+    Photon(theoryID, fiatlux_params, replicas=replicas, save_to_disk=True, force_computation=True)
 
 if __name__ == "__main__":
     main()
\ No newline at end of file
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index f7da3a4482..467dc259e9 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -66,7 +66,7 @@ class Photon:
         List of replica ids to be computed. If None, all replicas
         will be computed based on the luxqed pdf set.
     """
-    def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False):
+    def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False, force_computation=False):
         self.theoryid = theoryid
         self.lux_params = lux_params
         self.replicas = replicas
@@ -93,6 +93,11 @@ def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False):
         self.luxseed = lux_params["luxseed"]
         self.luxpdfset_members = self.luxpdfset.n_members - 1 # Remove replica 0
 
+        if force_computation:
+            self.compute_photon_set()
+            return
+            
+
         try:
           self.load_photon()
         except PhotonQEDNotFound:
@@ -271,13 +276,13 @@ def load_photon(self):
 
       # Load the needed replicas
       for replica in self.replicas:
-            # As input replica for the photon computation we take the MOD of the luxset_members to
-            # avoid failing due to limited number of replicas in the luxset
-            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+          # As input replica for the photon computation we take the MOD of the luxset_members to
+          # avoid failing due to limited number of replicas in the luxset
+          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
 
-            photon_array = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_array"]
-            interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
-            integral.append(trapezoid(photon_array, XGRID))
+          photon_array = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_array"]
+          interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
+          integral.append(trapezoid(photon_array, XGRID))
       
       self.interpolator = interpolator
       self.integral = np.stack(integral, axis=-1)

From d191408b523be647194bce6da2be0c7b2d2b2130 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 4 Nov 2025 17:50:07 +0000
Subject: [PATCH 06/30] Apply suggestions + download

---
 n3fit/src/n3fit/checks.py                     |  4 +--
 n3fit/src/n3fit/performfit.py                 |  1 -
 validphys2/src/validphys/loader.py            | 27 +++++++++----------
 .../src/validphys/nnprofile_default.yaml      |  2 +-
 4 files changed, 16 insertions(+), 18 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index fdc2ec0041..18727ce8d2 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -467,8 +467,8 @@ def check_photonQED_exists(theoryid, fiatlux):
       luxset = fiatlux['luxset']
       try:
           _ = FallbackLoader().check_photonQED(theoryid, luxset)
-      except FileNotFoundError:
-          raise CheckError(f"Photon QED set for TheoryID {theoryid.id} and luxset {luxset} not found.")
+      except FileNotFoundError as e:
+          raise CheckError(f"Photon QED set for TheoryID {theoryid.id} and luxset {luxset} not found.") from e
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py
index 312885de4c..d155b7431a 100644
--- a/n3fit/src/n3fit/performfit.py
+++ b/n3fit/src/n3fit/performfit.py
@@ -13,7 +13,6 @@
 
 # Action to be called by validphys
 # All information defining the NN should come here in the "parameters" dict
-@n3fit.checks.check_photonQED_exists
 @n3fit.checks.check_polarized_configs
 def performfit(
     *,
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index a7bd1bff6b..3111706527 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -217,12 +217,11 @@ def available_ekos(self):
             eko_path.parent.name.split("_")[1] for eko_path in self._theories_path.glob("*/eko.tar")
         }
     
-    @property
-    @functools.lru_cache
+    @functools.cached_property
     def available_photons_qed(self):
         """Return a string token for each of the available theories"""
         return {
-            eko_path.parent.name.split("_")[1] for eko_path in self._theories_path.glob("*/eko.npz")
+            photon_path.name.split("photon_")[1] for photon_path in self._photons_qed_path.glob("photon_*")
         }
 
     @property
@@ -324,7 +323,7 @@ def check_eko(self, theoryID):
     @functools.lru_cache
     def check_photonQED(self, theoryID, luxset):
         """Check the Photon QED set exists and return the path to it"""
-        photon_qed_path = self._photons_qed_path / f"photon_qed_{theoryID.id}_{luxset}"
+        photon_qed_path = self._photons_qed_path / f"photon_theoryID_{theoryID.id}_fit_{luxset}"
         if not photon_qed_path.exists():
             raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {theoryID}")
         return photon_qed_path
@@ -922,7 +921,7 @@ def remote_ekos(self):
     @property
     @functools.lru_cache
     def remote_photons_qed(self):
-        token = 'photon_qed_'
+        token = 'photon_'
         rt = self.remote_files(self.photon_qed_urls, self.photon_qed_index, thing="photons_qed")
         return {k[len(token) :]: v for k, v in rt.items()}
 
@@ -1148,15 +1147,15 @@ def download_eko(self, thid):
 
     def download_photonQED(self, thid, luxset: str):
         """Download the Photon set for a given theory ID"""
-        # thid = str(thid)
-        # remote = self.remote_photons_qed
-        # key = f"{thid}_{luxset}"
-        # if key not in remote:
-        #     raise PhotonQEDNotFound(f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server")
-        # # Check that we have the theory we need
-        # target_path = self._photon_qed_path / f"photon_qed_{int(thid)}_{luxset}.tar"
-        # download_file(remote[key], target_path)
-        log.warning("Downloading Photon QED sets is not implemented yet.")
+        thid = thid.id
+        remote = self.remote_photons_qed
+        key = f"theoryID_{thid}_fit_{luxset}"
+        if key not in remote:
+            raise PhotonQEDNotFound(f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server")
+        # Check that we have the theory we need
+        target_path = self._photons_qed_path
+        download_and_extract(remote[key], target_path)
+        
     
     def download_vp_output_file(self, filename, **kwargs):
         try:
diff --git a/validphys2/src/validphys/nnprofile_default.yaml b/validphys2/src/validphys/nnprofile_default.yaml
index c2d7e77965..f09ffe9559 100644
--- a/validphys2/src/validphys/nnprofile_default.yaml
+++ b/validphys2/src/validphys/nnprofile_default.yaml
@@ -63,7 +63,7 @@ eko_urls:
 eko_index: 'ekodata.json'
 
 photon_qed_urls:
-    - 'https://nnpdf.nikhef.nl/nnpdf/photons_qed/'
+    - 'https://data.nnpdf.science/photons/'
 
 photon_qed_index: 'photondata.json'
 

From 72d419a449868706d5f6610bd5a56cc59e5dd3ae Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 4 Nov 2025 18:44:57 +0000
Subject: [PATCH 07/30] Restore check for photon [skip ci]

---
 n3fit/src/n3fit/checks.py     | 5 +++--
 n3fit/src/n3fit/performfit.py | 1 +
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index 18727ce8d2..0ecd92beed 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -467,8 +467,9 @@ def check_photonQED_exists(theoryid, fiatlux):
       luxset = fiatlux['luxset']
       try:
           _ = FallbackLoader().check_photonQED(theoryid, luxset)
-      except FileNotFoundError as e:
-          raise CheckError(f"Photon QED set for TheoryID {theoryid.id} and luxset {luxset} not found.") from e
+      except FileNotFoundError:
+          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
+                      "will be compute using FiatLux. This may impact performance.")
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py
index d155b7431a..312885de4c 100644
--- a/n3fit/src/n3fit/performfit.py
+++ b/n3fit/src/n3fit/performfit.py
@@ -13,6 +13,7 @@
 
 # Action to be called by validphys
 # All information defining the NN should come here in the "parameters" dict
+@n3fit.checks.check_photonQED_exists
 @n3fit.checks.check_polarized_configs
 def performfit(
     *,

From bdc7042d5463a711678caff404272b18178a4c74 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 4 Nov 2025 19:18:20 +0000
Subject: [PATCH 08/30] Fix compute + index photon

---
 validphys2/serverscripts/index-photon.py   | 20 ++++++++++++++++++++
 validphys2/src/validphys/photon/compute.py |  2 +-
 2 files changed, 21 insertions(+), 1 deletion(-)
 create mode 100755 validphys2/serverscripts/index-photon.py

diff --git a/validphys2/serverscripts/index-photon.py b/validphys2/serverscripts/index-photon.py
new file mode 100755
index 0000000000..927fdb5d2d
--- /dev/null
+++ b/validphys2/serverscripts/index-photon.py
@@ -0,0 +1,20 @@
+# -*- coding: utf-8 -*-
+"""
+Generate an index with the existing internal or unpublished PDFS.
+"""
+
+import pathlib
+import json
+
+root = '/home/nnpdf/WEB/photons'
+
+glob = '*.tar'
+
+indexname = 'photondata.json'
+
+if __name__ == '__main__':
+    p = pathlib.Path(root)
+    files = p.glob(glob)
+    files = [f.name for f in files]
+    with (p/indexname).open('w') as f:
+        json.dump({'files':files}, f, separators=(',',':'))
\ No newline at end of file
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 467dc259e9..29b05e9e6e 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -169,7 +169,7 @@ def evaluate_at_x(x):
             with ThreadPoolExecutor() as executor:
                 photon_qin = np.array(list(executor.map(evaluate_at_x, XGRID)))
 
-            # photon_qin += self.generate_errors(replica)
+            photon_qin += self.generate_errors(replica)
 
             # fiatlux computes x * gamma(x)
             photon_qin /= XGRID

From a731ac6e20ae442d9f14966037f2ae5209c1e19f Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 4 Nov 2025 19:20:03 +0000
Subject: [PATCH 09/30] Change name when photon is computed

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

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 29b05e9e6e..8ad13786f7 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -201,7 +201,7 @@ def evaluate_at_x(x):
             photon_array = XGRID * photon_Q0
 
             if self.save_to_disk:
-                path_to_photon = loader._photons_qed_path / f"photon_qed_{self.theoryid.id}_{self.luxpdfset._name}"
+                path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
                 path_to_photon.mkdir(parents=True, exist_ok=True)
                 np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
                 log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")

From 3a0c2098e3cb10c47ac5e6ecdd32a57dae8d3c7a Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 5 Nov 2025 17:05:37 +0000
Subject: [PATCH 10/30] Addtional errors added in fiatlux exec

---
 n3fit/src/n3fit/scripts/fiatlux_exec.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/n3fit/src/n3fit/scripts/fiatlux_exec.py b/n3fit/src/n3fit/scripts/fiatlux_exec.py
index 735bcf3957..e0ca202051 100755
--- a/n3fit/src/n3fit/scripts/fiatlux_exec.py
+++ b/n3fit/src/n3fit/scripts/fiatlux_exec.py
@@ -39,7 +39,8 @@ def main():
     theoryID = loader.check_theoryID(runcard["theory"]["theoryid"])
     fiatlux_params = runcard.get("fiatlux", None)
     fiatlux_params['luxset'] = loader.check_pdf(fiatlux_params['luxset'])
-
+    if fiatlux_params.get('additional_errors', False):
+        fiatlux_params['additional_errors'] = loader.check_pdf("LUXqed17_plus_PDF4LHC15_nnlo_100")
     Photon(theoryID, fiatlux_params, replicas=replicas, save_to_disk=True, force_computation=True)
 
 if __name__ == "__main__":

From 541bbd186ecff005899162218fde08d5455a5fbc Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 6 Nov 2025 09:40:53 +0000
Subject: [PATCH 11/30] Fixing names for photon resource

---
 n3fit/src/n3fit/checks.py                  |  3 ++-
 validphys2/src/validphys/loader.py         | 17 ++++++++---------
 validphys2/src/validphys/photon/compute.py |  2 +-
 3 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index 0ecd92beed..05282fc79b 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -466,7 +466,8 @@ def check_photonQED_exists(theoryid, fiatlux):
     if fiatlux is not None:
       luxset = fiatlux['luxset']
       try:
-          _ = FallbackLoader().check_photonQED(theoryid, luxset)
+          _ = FallbackLoader().check_photonQED(theoryid.id, luxset)
+          log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
       except FileNotFoundError:
           log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
                       "will be compute using FiatLux. This may impact performance.")
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index 3111706527..1aadfcbc6b 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -218,7 +218,7 @@ def available_ekos(self):
         }
     
     @functools.cached_property
-    def available_photons_qed(self):
+    def available_photons(self):
         """Return a string token for each of the available theories"""
         return {
             photon_path.name.split("photon_")[1] for photon_path in self._photons_qed_path.glob("photon_*")
@@ -323,9 +323,9 @@ def check_eko(self, theoryID):
     @functools.lru_cache
     def check_photonQED(self, theoryID, luxset):
         """Check the Photon QED set exists and return the path to it"""
-        photon_qed_path = self._photons_qed_path / f"photon_theoryID_{theoryID.id}_fit_{luxset}"
+        photon_qed_path = self._photons_qed_path / f"photon_theoryID_{int(theoryID)}_fit_{luxset}"
         if not photon_qed_path.exists():
-            raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {theoryID}")
+            raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {int(theoryID)}")
         return photon_qed_path
 
     @property
@@ -920,9 +920,9 @@ def remote_ekos(self):
     
     @property
     @functools.lru_cache
-    def remote_photons_qed(self):
+    def remote_photons(self):
         token = 'photon_'
-        rt = self.remote_files(self.photon_qed_urls, self.photon_qed_index, thing="photons_qed")
+        rt = self.remote_files(self.photon_qed_urls, self.photon_qed_index, thing="photons")
         return {k[len(token) :]: v for k, v in rt.items()}
 
     @property
@@ -960,8 +960,8 @@ def downloadable_ekos(self):
         return list(self.remote_ekos)
     
     @property
-    def downloadable_photonsQED(self):
-        return list(self.remote_photons_qed)
+    def downloadable_photons(self):
+        return list(self.remote_photons)
 
     @property
     def lhapdf_pdfs(self):
@@ -1147,8 +1147,7 @@ def download_eko(self, thid):
 
     def download_photonQED(self, thid, luxset: str):
         """Download the Photon set for a given theory ID"""
-        thid = thid.id
-        remote = self.remote_photons_qed
+        remote = self.remote_photons
         key = f"theoryID_{thid}_fit_{luxset}"
         if key not in remote:
             raise PhotonQEDNotFound(f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server")
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 8ad13786f7..9120de5edf 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -268,7 +268,7 @@ def generate_errors(self, replica):
     
     def load_photon(self):
       """Load the photon resource using the Loader class."""
-      path_to_photon = loader.check_photonQED(self.theoryid, self.luxpdfset._name)
+      path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
       log.info(f"Loading photon QED set from {path_to_photon}")
 
       interpolator = []

From 1dc141ab0320675e23a04029f4d9e92958c14c4e Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 6 Nov 2025 10:40:02 +0000
Subject: [PATCH 12/30] Make photon downloadable from vp-setupfit only

---
 n3fit/src/n3fit/checks.py                  | 15 +++++----------
 n3fit/src/n3fit/scripts/vp_setupfit.py     |  2 +-
 validphys2/src/validphys/filters.py        | 18 ++++++++++++++++++
 validphys2/src/validphys/photon/compute.py |  1 -
 4 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index 05282fc79b..6307761bae 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -9,7 +9,7 @@
 from n3fit.hyper_optimization import penalties as penalties_module
 from n3fit.hyper_optimization.rewards import IMPLEMENTED_LOSSES, IMPLEMENTED_STATS
 from reportengine.checks import CheckError, make_argcheck
-from validphys.loader import FallbackLoader
+from validphys.loader import FallbackLoader, Loader
 from validphys.pdfbases import check_basis
 
 log = logging.getLogger(__name__)
@@ -452,25 +452,20 @@ def check_deprecated_options(fitting):
     for option in nnfit_options:
         if option in fitting:
             log.warning("'fitting::%s' is an nnfit-only key, it will be ignored", option)
-
-
-@make_argcheck
-def check_multireplica_qed(replicas, fiatlux):
-    if fiatlux is not None:
-        if len(replicas) > 1:
-            raise CheckError("At the moment, running a multireplica QED fits is not allowed.")
         
+
 @make_argcheck
 def check_photonQED_exists(theoryid, fiatlux):
     """Check that the Photon QED set for this theoryid and luxset exists"""
     if fiatlux is not None:
       luxset = fiatlux['luxset']
       try:
-          _ = FallbackLoader().check_photonQED(theoryid.id, luxset)
+          _ = Loader().check_photonQED(theoryid.id, luxset)
           log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
       except FileNotFoundError:
           log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
-                      "will be compute using FiatLux. This may impact performance.")
+                      "will be compute using FiatLux. This may impact performance. Make"
+                      "sure vp-setupfit has been run prior to the fit to download necessary resources.")
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index d63429cc70..0fa7094d48 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -184,7 +184,7 @@ def from_yaml(cls, o, *args, **kwargs):
         # Check fiatlux configuration
         fiatlux = file_content.get('fiatlux')
         if fiatlux is not None:
-            SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset')
+            SETUPFIT_FIXED_CONFIG['actions_'].append('theory::fiatlux check_photonQED_exists')
             if fiatlux.get("additional_errors"):
                 SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
 
diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py
index add8a3bf4b..35b87195c9 100644
--- a/validphys2/src/validphys/filters.py
+++ b/validphys2/src/validphys/filters.py
@@ -479,6 +479,24 @@ def check_luxset(luxset):
     log.info(f'{luxset} Lux pdf checked.')
 
 
+def check_photonQED_exists(theoryid, fiatlux):
+    """Check that the Photon QED set for this theoryid and luxset exists"""
+    if fiatlux is not None:
+      from validphys.loader import FallbackLoader
+      luxset = fiatlux['luxset']
+      try:
+          _ = FallbackLoader().check_photonQED(theoryid.id, luxset)
+          log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
+      except FileNotFoundError:
+          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
+                      "will be compute using FiatLux. This may impact performance.")
+          
+          # Since the photon is missing and will be computed on the fly, check
+          # the luxset exists
+          luxset.load()
+          log.info(f'{luxset} Lux pdf checked.')
+
+
 def check_unpolarized_bc(unpolarized_bc):
     """Check that unpolarized PDF bound can be loaded normally."""
     unpolarized_bc.load()
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 9120de5edf..e3f46a1bda 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -97,7 +97,6 @@ def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False, forc
             self.compute_photon_set()
             return
             
-
         try:
           self.load_photon()
         except PhotonQEDNotFound:

From 5eaa6687864b1d98aaf3092c6af599ed7a626474 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 12 Nov 2025 12:31:13 +0000
Subject: [PATCH 13/30] Moving photon computation into vp-setupfit

---
 n3fit/src/n3fit/scripts/vp_setupfit.py     | 29 +++++++++++----
 validphys2/src/validphys/filters.py        | 16 ++-------
 validphys2/src/validphys/photon/compute.py | 42 +++++++++++++++-------
 3 files changed, 56 insertions(+), 31 deletions(-)

diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index 0fa7094d48..a365a1588d 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -37,10 +37,10 @@
 from reportengine import colors
 from validphys.app import App
 from validphys.config import Config, ConfigError, Environment, EnvironmentError_
-from validphys.loader import Loader, TheoryNotFound
+from validphys.loader import FallbackLoader, TheoryNotFound, PhotonQEDNotFound
 from validphys.utils import yaml_safe
 
-l = Loader()
+loader = FallbackLoader()
 
 SETUPFIT_FIXED_CONFIG = dict(
     actions_=[
@@ -57,6 +57,7 @@
     'validphys.filters',
     'validphys.results',
     'validphys.theorycovariance.construction',
+    'validphys.photon.compute'
 ]
 
 SETUPFIT_DEFAULTS = dict(use_cuts='internal')
@@ -155,7 +156,7 @@ def from_yaml(cls, o, *args, **kwargs):
                 SETUPFIT_FIXED_CONFIG['theory']['theoryid'] = closuredict['faketheoryid']
                 # download theoryid since it will be used in the fit
                 try:
-                    l.check_theoryID(file_content['theory']['theoryid'])
+                    loader.check_theoryID(file_content['theory']['theoryid'])
                 except TheoryNotFound as e:
                     log.warning(e)
             filter_action = 'datacuts::closuretest::theory::fitting filter'
@@ -184,9 +185,25 @@ def from_yaml(cls, o, *args, **kwargs):
         # Check fiatlux configuration
         fiatlux = file_content.get('fiatlux')
         if fiatlux is not None:
-            SETUPFIT_FIXED_CONFIG['actions_'].append('theory::fiatlux check_photonQED_exists')
-            if fiatlux.get("additional_errors"):
-                SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
+            luxset = fiatlux['luxset']
+            theoryid = file_content['theory']['theoryid']
+            force_compute = fiatlux.get('force_computation', False)
+            try:
+                _ = loader.check_photonQED(theoryid, luxset)
+                log.info(f"Photon QED set found for {theoryid} with luxset {luxset}.")
+            except PhotonQEDNotFound:
+                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. It "\
+                            "will be computed using FiatLux. This may impact performance.")
+                force_compute = True
+
+            if force_compute:
+              log.info("Forcing photon computation with FiatLux.")
+              # Since the photon will be computed, check that the luxset and additional_errors exist
+              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset_exists')
+              if fiatlux.get("additional_errors"):
+                  SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
+              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux::theory compute_photon')
+
 
         # Check positivity bound
         if file_content.get('positivity_bound') is not None:
diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py
index 35b87195c9..7cedcea9a0 100644
--- a/validphys2/src/validphys/filters.py
+++ b/validphys2/src/validphys/filters.py
@@ -479,22 +479,12 @@ def check_luxset(luxset):
     log.info(f'{luxset} Lux pdf checked.')
 
 
-def check_photonQED_exists(theoryid, fiatlux):
+def check_luxset_exists(fiatlux):
     """Check that the Photon QED set for this theoryid and luxset exists"""
     if fiatlux is not None:
-      from validphys.loader import FallbackLoader
       luxset = fiatlux['luxset']
-      try:
-          _ = FallbackLoader().check_photonQED(theoryid.id, luxset)
-          log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
-      except FileNotFoundError:
-          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
-                      "will be compute using FiatLux. This may impact performance.")
-          
-          # Since the photon is missing and will be computed on the fly, check
-          # the luxset exists
-          luxset.load()
-          log.info(f'{luxset} Lux pdf checked.')
+      luxset.load()
+      log.info(f'{luxset} Lux pdf checked.')
 
 
 def check_unpolarized_bc(unpolarized_bc):
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index e3f46a1bda..0185bfa2bf 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -3,6 +3,7 @@
 import logging
 import tempfile
 from concurrent.futures import ThreadPoolExecutor
+from joblib import Parallel, delayed, effective_n_jobs
 
 import numpy as np
 from scipy.integrate import trapezoid
@@ -125,10 +126,8 @@ def compute_photon_set(self):
         # set fiatlux
         mb_thr = theory["kbThr"] * theory["mb"]
         mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
-        interpolator = []
-        integral = []
 
-        for replica in self.replicas:
+        def process_single_replica(replica):
             # As input replica for the photon computation we take the MOD of the luxset_members to
             # avoid failing due to limited number of replicas in the luxset
             photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
@@ -198,19 +197,29 @@ def evaluate_at_x(x):
 
             photon_Q0 = pdfs_final[ph_id]
             photon_array = XGRID * photon_Q0
+            return (photon_array, photonreplica)
+        
 
-            if self.save_to_disk:
-                path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
-                path_to_photon.mkdir(parents=True, exist_ok=True)
-                np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
-                log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
+        log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
+        photon_tuples = Parallel(n_jobs=-1)(delayed(process_single_replica(replica)) for replica in self.replicas)
 
-            interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
-            integral.append(trapezoid(photon_array, XGRID))
+        interpolator = []
+        integral = []
 
-        integral = np.stack(integral, axis=-1)
+        # If saving to disk, create the directory
+        if self.save_to_disk:
+          path_to_photon.mkdir(parents=True, exist_ok=True)
 
+        for photon_array, photonreplica in photon_tuples:
+          if self.save_to_disk:
+              path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
+              np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
+              log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
 
+          interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
+          integral.append(trapezoid(photon_array, XGRID))
+
+        integral = np.stack(integral, axis=-1)
         self.integral = integral
         self.interpolator = interpolator
 
@@ -285,4 +294,13 @@ def load_photon(self):
       
       self.interpolator = interpolator
       self.integral = np.stack(integral, axis=-1)
-      return
\ No newline at end of file
+      return
+    
+def compute_photon(theoryid, fiatlux, force_fiatlux=False):
+    """Function to compute the photon PDF set.""" 
+    luxset = fiatlux['luxset'].load()
+    replicas = list(range(1, luxset.n_members))
+    Photon(theoryid, fiatlux, replicas=replicas, save_to_disk=True, force_computation=force_fiatlux)
+
+
+    
\ No newline at end of file

From 9efa4e2c64747896702ab8b6145602b9f4611039 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 12 Nov 2025 21:52:47 +0000
Subject: [PATCH 14/30] force_compute only if specified by the user

---
 n3fit/src/n3fit/scripts/vp_setupfit.py     | 8 ++++----
 validphys2/src/validphys/photon/compute.py | 5 +++--
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index a365a1588d..c92159a11c 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -187,14 +187,14 @@ def from_yaml(cls, o, *args, **kwargs):
         if fiatlux is not None:
             luxset = fiatlux['luxset']
             theoryid = file_content['theory']['theoryid']
-            force_compute = fiatlux.get('force_computation', False)
+            force_compute = fiatlux.get('compute_in_setupfit', False)
             try:
                 _ = loader.check_photonQED(theoryid, luxset)
                 log.info(f"Photon QED set found for {theoryid} with luxset {luxset}.")
             except PhotonQEDNotFound:
-                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. It "\
-                            "will be computed using FiatLux. This may impact performance.")
-                force_compute = True
+                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Set "\
+                            "compute_in_setupfit to true in the runcard. Otherwise n3fit" \
+                            "will take care of the photon computation. May impact performance.")
 
             if force_compute:
               log.info("Forcing photon computation with FiatLux.")
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 0185bfa2bf..88112beb77 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -296,11 +296,12 @@ def load_photon(self):
       self.integral = np.stack(integral, axis=-1)
       return
     
-def compute_photon(theoryid, fiatlux, force_fiatlux=False):
+def compute_photon(theoryid, fiatlux):
     """Function to compute the photon PDF set.""" 
     luxset = fiatlux['luxset'].load()
+    force_compute = fiatlux.get('compute_in_setupfit', False)
     replicas = list(range(1, luxset.n_members))
-    Photon(theoryid, fiatlux, replicas=replicas, save_to_disk=True, force_computation=force_fiatlux)
+    Photon(theoryid, fiatlux, replicas=replicas, save_to_disk=True, force_computation=force_compute)
 
 
     
\ No newline at end of file

From b5e51806293f15184fc4a751862610e7f9297a11 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 12 Nov 2025 22:11:57 +0000
Subject: [PATCH 15/30] Correct path to photon

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

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 88112beb77..1b600ed339 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -208,11 +208,11 @@ def evaluate_at_x(x):
 
         # If saving to disk, create the directory
         if self.save_to_disk:
+          path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
           path_to_photon.mkdir(parents=True, exist_ok=True)
 
         for photon_array, photonreplica in photon_tuples:
           if self.save_to_disk:
-              path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
               np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
               log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
 

From 5b59615a269ebfecedee42057dacf693a1be10ba Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 12 Nov 2025 22:40:38 +0000
Subject: [PATCH 16/30] Add photon path to profile utils

---
 nnpdf_data/nnpdf_data/utils.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/nnpdf_data/nnpdf_data/utils.py b/nnpdf_data/nnpdf_data/utils.py
index fea46e285a..7013cfdbf3 100644
--- a/nnpdf_data/nnpdf_data/utils.py
+++ b/nnpdf_data/nnpdf_data/utils.py
@@ -102,6 +102,7 @@ def get_nnpdf_profile(profile_path=None):
         "validphys_cache_path",
         "hyperscan_path",
         "ekos_path",
+        "photons_qed_path",
     ]:
         # if there are any problems setting or getting these variable erroring out is more than justified
         absolute_var = nnpdf_share / pathlib.Path(profile_dict[var]).expanduser()

From 169f6c498c5a3ba36970a05e77f1b75e549ca7ba Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 13 Nov 2025 09:01:32 +0000
Subject: [PATCH 17/30] Move joblib outside the Photon class

---
 n3fit/src/n3fit/checks.py                  |   6 +-
 n3fit/src/n3fit/scripts/vp_setupfit.py     |   5 +-
 validphys2/src/validphys/photon/compute.py | 147 ++++++++++-----------
 3 files changed, 75 insertions(+), 83 deletions(-)

diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index 6307761bae..8424181c5f 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -463,9 +463,9 @@ def check_photonQED_exists(theoryid, fiatlux):
           _ = Loader().check_photonQED(theoryid.id, luxset)
           log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
       except FileNotFoundError:
-          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset} and "\
-                      "will be compute using FiatLux. This may impact performance. Make"
-                      "sure vp-setupfit has been run prior to the fit to download necessary resources.")
+          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset}. It "\
+                      "will be compute using FiatLux. This may impact performance. Make "
+                      "sure vp-setupfit has been run prior to n3fit to download necessary resources.")
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index c92159a11c..6234c9ba82 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -192,8 +192,9 @@ def from_yaml(cls, o, *args, **kwargs):
                 _ = loader.check_photonQED(theoryid, luxset)
                 log.info(f"Photon QED set found for {theoryid} with luxset {luxset}.")
             except PhotonQEDNotFound:
-                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Set "\
-                            "compute_in_setupfit to true in the runcard. Otherwise n3fit" \
+                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Consider "\
+                            "using `compute_in_setupfit` in the runcard to compute all replicas of the photon "
+                            "in vp-setupfit. Otherwise n3fit " \
                             "will take care of the photon computation. May impact performance.")
 
             if force_compute:
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 1b600ed339..2aa188d851 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -127,82 +127,6 @@ def compute_photon_set(self):
         mb_thr = theory["kbThr"] * theory["mb"]
         mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
 
-        def process_single_replica(replica):
-            # As input replica for the photon computation we take the MOD of the luxset_members to
-            # avoid failing due to limited number of replicas in the luxset
-            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
-
-            f2lo = sf.F2LO(self.luxpdfset.members[photonreplica], theory)
-
-            if theory["PTO"] > 0:
-                f2 = sf.InterpStructureFunction(path_to_F2, self.luxpdfset.members[photonreplica])
-                fl = sf.InterpStructureFunction(path_to_FL, self.luxpdfset.members[photonreplica])
-                if not np.isclose(f2.q2_max, fl.q2_max):
-                    log.error(
-                        "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
-                    )
-                self.fiatlux_runcard["q2_max"] = float(f2.q2_max)
-            else:
-                f2 = f2lo
-                fl = sf.FLLO()
-                # using a default value for q2_max
-                self.fiatlux_runcard["q2_max"] = 1e8
-
-            alpha = Alpha(theory, self.fiatlux_runcard["q2_max"])
-            with tempfile.NamedTemporaryFile(mode="w") as tmp:
-                yaml.dump(self.fiatlux_runcard, tmp)
-                lux = fiatlux.FiatLux(tmp.name)
-
-            # we have a dict but fiatlux wants a yaml file
-            # TODO : once that fiatlux will allow dictionaries
-            # pass directly fiatlux_runcard
-            lux.PlugAlphaQED(alpha.alpha_em, alpha.qref)
-            lux.InsertInelasticSplitQ([mb_thr, mt_thr])
-            lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
-
-            # Evaluate photon for every point in the grid xgrid
-            def evaluate_at_x(x):
-              return lux.EvaluatePhoton(x, self.q_in**2).total
-            
-            with ThreadPoolExecutor() as executor:
-                photon_qin = np.array(list(executor.map(evaluate_at_x, XGRID)))
-
-            photon_qin += self.generate_errors(replica)
-
-            # fiatlux computes x * gamma(x)
-            photon_qin /= XGRID
-
-            # Load eko and reshape it
-            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
-                
-                # 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
-                        else:
-                            pdfs_init[j] = np.array(
-                                [self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x for x in XGRID]
-                            )
-
-                    pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
-
-            photon_Q0 = pdfs_final[ph_id]
-            photon_array = XGRID * photon_Q0
-            return (photon_array, photonreplica)
-        
-
-        log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
-        photon_tuples = Parallel(n_jobs=-1)(delayed(process_single_replica(replica)) for replica in self.replicas)
-
         interpolator = []
         integral = []
 
@@ -211,7 +135,71 @@ def evaluate_at_x(x):
           path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
           path_to_photon.mkdir(parents=True, exist_ok=True)
 
-        for photon_array, photonreplica in photon_tuples:
+        for replica in self.replicas:
+          # As input replica for the photon computation we take the MOD of the luxset_members to
+          # avoid failing due to limited number of replicas in the luxset
+          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+
+          f2lo = sf.F2LO(self.luxpdfset.members[photonreplica], theory)
+
+          if theory["PTO"] > 0:
+              f2 = sf.InterpStructureFunction(path_to_F2, self.luxpdfset.members[photonreplica])
+              fl = sf.InterpStructureFunction(path_to_FL, self.luxpdfset.members[photonreplica])
+              if not np.isclose(f2.q2_max, fl.q2_max):
+                  log.error(
+                      "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
+                  )
+              self.fiatlux_runcard["q2_max"] = float(f2.q2_max)
+          else:
+              f2 = f2lo
+              fl = sf.FLLO()
+              # using a default value for q2_max
+              self.fiatlux_runcard["q2_max"] = 1e8
+
+          alpha = Alpha(theory, self.fiatlux_runcard["q2_max"])
+          with tempfile.NamedTemporaryFile(mode="w") as tmp:
+              yaml.dump(self.fiatlux_runcard, tmp)
+              lux = fiatlux.FiatLux(tmp.name)
+
+          # we have a dict but fiatlux wants a yaml file
+          # TODO : once that fiatlux will allow dictionaries
+          # pass directly fiatlux_runcard
+          lux.PlugAlphaQED(alpha.alpha_em, alpha.qref)
+          lux.InsertInelasticSplitQ([mb_thr, mt_thr])
+          lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
+
+          photon_qin = np.array(lux.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
+
+          # Load eko and reshape it
+          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
+              
+              # 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
+                      else:
+                          pdfs_init[j] = np.array(
+                              [self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x for x in XGRID]
+                          )
+
+                  pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
+
+          photon_Q0 = pdfs_final[ph_id]
+          photon_array = XGRID * photon_Q0
+
           if self.save_to_disk:
               np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
               log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
@@ -301,7 +289,10 @@ def compute_photon(theoryid, fiatlux):
     luxset = fiatlux['luxset'].load()
     force_compute = fiatlux.get('compute_in_setupfit', False)
     replicas = list(range(1, luxset.n_members))
-    Photon(theoryid, fiatlux, replicas=replicas, save_to_disk=True, force_computation=force_compute)
+
+    log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
+    _ = Parallel(n_jobs=-1)(delayed(Photon)(theoryid, fiatlux, replicas=[replica], save_to_disk=True, force_computation=force_compute) for replica in replicas)
+    
 
 
     
\ No newline at end of file

From 9464d0bfa06638c029fc65ac65021a79eef980d7 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 13 Nov 2025 09:10:01 +0000
Subject: [PATCH 18/30] Fix shape in photon computation

---
 n3fit/src/n3fit/scripts/fiatlux_exec.py    | 47 ----------------------
 validphys2/src/validphys/photon/compute.py |  4 +-
 2 files changed, 3 insertions(+), 48 deletions(-)
 delete mode 100755 n3fit/src/n3fit/scripts/fiatlux_exec.py

diff --git a/n3fit/src/n3fit/scripts/fiatlux_exec.py b/n3fit/src/n3fit/scripts/fiatlux_exec.py
deleted file mode 100755
index e0ca202051..0000000000
--- a/n3fit/src/n3fit/scripts/fiatlux_exec.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import argparse
-import pathlib
-import logging
-
-from validphys.loader import FallbackLoader as Loader
-from validphys.utils import yaml_safe
-from validphys.photon.compute import Photon
-
-def check_positive(value):
-    ivalue = int(value)
-    if ivalue <= 0:
-        raise argparse.ArgumentTypeError("%s is an invalid positive int value." % value)
-    return ivalue
-
-log = logging.getLogger(__name__)
-loader = Loader()
-
-
-
-def main():
-    parser = argparse.ArgumentParser(description="FitLux-exec - compute the photon PDF using FiatLux")
-    parser.add_argument("config_yml", help="Path to the configuration file", type=pathlib.Path)
-    parser.add_argument("replica", help="MC replica number", type=check_positive)
-    parser.add_argument(
-        "-r",
-        "--replica_range",
-        help="End of the range of replicas to compute",
-        type=check_positive,
-    )
-    args = parser.parse_args()
-    config_path = args.config_yml.absolute()
-    replica = args.replica
-    if args.replica_range:
-        replicas = list(range(replica, args.replica_range + 1))
-    else:
-        replicas = [replica]
-
-    runcard = yaml_safe.load(config_path.read_text(encoding="UTF-8"))
-    theoryID = loader.check_theoryID(runcard["theory"]["theoryid"])
-    fiatlux_params = runcard.get("fiatlux", None)
-    fiatlux_params['luxset'] = loader.check_pdf(fiatlux_params['luxset'])
-    if fiatlux_params.get('additional_errors', False):
-        fiatlux_params['additional_errors'] = loader.check_pdf("LUXqed17_plus_PDF4LHC15_nnlo_100")
-    Photon(theoryID, fiatlux_params, replicas=replicas, save_to_disk=True, force_computation=True)
-
-if __name__ == "__main__":
-    main()
\ No newline at end of file
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 2aa188d851..be2111cdcc 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -168,7 +168,9 @@ def compute_photon_set(self):
           lux.InsertInelasticSplitQ([mb_thr, mt_thr])
           lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
 
-          photon_qin = np.array(lux.EvaluatePhoton(x, self.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)

From 32ec19e366b750fbc1e13538fc70cac54dad9685 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 13 Nov 2025 09:28:42 +0000
Subject: [PATCH 19/30] Fixing self.lux -> lux

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

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index be2111cdcc..7f9f48bb20 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -169,7 +169,7 @@ def compute_photon_set(self):
           lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
 
           photon_qin = np.array(
-            [self.lux[replica].EvaluatePhoton(x, self.q_in**2).total for x in XGRID]
+            [lux.EvaluatePhoton(x, self.q_in**2).total for x in XGRID]
         )
           photon_qin += self.generate_errors(replica)
 

From ff13285600e795406bea7e9b42055782c1842dca Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 13 Nov 2025 10:17:51 +0000
Subject: [PATCH 20/30] Remove unused imports

---
 validphys2/src/validphys/photon/compute.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 7f9f48bb20..aef4fcbb1c 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -2,7 +2,6 @@
 
 import logging
 import tempfile
-from concurrent.futures import ThreadPoolExecutor
 from joblib import Parallel, delayed, effective_n_jobs
 
 import numpy as np

From 1fb51c833ba3db8668958b19e8476dd340f80453 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 13 Nov 2025 11:22:23 +0000
Subject: [PATCH 21/30] Avoid pickling issues with joblib

---
 validphys2/src/validphys/photon/compute.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index aef4fcbb1c..34724b2b42 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -291,9 +291,10 @@ def compute_photon(theoryid, fiatlux):
     force_compute = fiatlux.get('compute_in_setupfit', False)
     replicas = list(range(1, luxset.n_members))
 
-    log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
-    _ = Parallel(n_jobs=-1)(delayed(Photon)(theoryid, fiatlux, replicas=[replica], save_to_disk=True, force_computation=force_compute) for replica in replicas)
+    # Return None and avoid pickling issues with the phoyon class.
+    def wrapper_fn(replica, theoryid, fiatlux, force_compute):
+        Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=True, force_computation=force_compute)
+        return None
     
-
-
-    
\ No newline at end of file
+    log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
+    _ = Parallel(n_jobs=-1)(delayed(wrapper_fn)(replica=replica, theoryid=theoryid, fiatlux=fiatlux, force_compute=force_compute) for replica in replicas)
\ No newline at end of file

From 38dfd1bc99574eef9ceae2d8c1209ee55304476c Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 18 Nov 2025 11:28:49 +0000
Subject: [PATCH 22/30] Remove fiatlux-exec from pyproject

---
 pyproject.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index acb0f039da..5bdc162905 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -55,8 +55,6 @@ vp-list = "validphys.scripts.vp_list:main"
 vp-nextfitruncard = "validphys.scripts.vp_nextfitruncard:main"
 vp-hyperoptplot = "validphys.scripts.vp_hyperoptplot:main"
 vp-deltachi2 = "validphys.scripts.vp_deltachi2:main"
-# FiatLux scripts
-fiatlux-run = "n3fit.scripts.fiatlux_exec:main"
 
 [tool.poetry.dependencies]
 # Generic dependencies (i.e., validphys)

From d969ca760ea3bd9067ab82579d55956a1d6cc3be Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 18 Nov 2025 11:37:31 +0000
Subject: [PATCH 23/30] Move loader initialization inside function calls

---
 validphys2/src/validphys/photon/compute.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 34724b2b42..11f07307f0 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -19,7 +19,6 @@
 from .alpha import Alpha
 
 log = logging.getLogger(__name__)
-loader = Loader()
 
 # not the complete fiatlux runcard since some parameters are set in the code
 FIATLUX_DEFAULT = {
@@ -131,6 +130,7 @@ def compute_photon_set(self):
 
         # If saving to disk, create the directory
         if self.save_to_disk:
+          loader = Loader()
           path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
           path_to_photon.mkdir(parents=True, exist_ok=True)
 
@@ -265,6 +265,7 @@ def generate_errors(self, replica):
     
     def load_photon(self):
       """Load the photon resource using the Loader class."""
+      loader = Loader()
       path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
       log.info(f"Loading photon QED set from {path_to_photon}")
 

From 4860d8e3a22fd46c20716d6760c7f49aa6ccca62 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 18 Nov 2025 11:42:35 +0000
Subject: [PATCH 24/30] Add check for non-photon resources

---
 validphys2/src/validphys/scripts/vp_get.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/validphys2/src/validphys/scripts/vp_get.py b/validphys2/src/validphys/scripts/vp_get.py
index 70f06a87b0..2323db834d 100644
--- a/validphys2/src/validphys/scripts/vp_get.py
+++ b/validphys2/src/validphys/scripts/vp_get.py
@@ -57,6 +57,9 @@ def main():
     tp = args.resource_type
     name = args.resource_name
 
+    if len(name) > 1 and tp != 'photonQED':
+        sys.exit("Only 'photonQED' resource type requires theoryID and luxset.")
+
     try:
         f = getattr(l, f'check_{tp}')
     except AttributeError as e:

From fe7330a4078349b233766ef93b88d021d86ebd26 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 18 Nov 2025 12:02:38 +0000
Subject: [PATCH 25/30] Use maxcores + rename function that computes the photon

---
 n3fit/src/n3fit/scripts/vp_setupfit.py     | 2 +-
 validphys2/src/validphys/photon/compute.py | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index 6234c9ba82..2aee06f900 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -203,7 +203,7 @@ def from_yaml(cls, o, *args, **kwargs):
               SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset_exists')
               if fiatlux.get("additional_errors"):
                   SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
-              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux::theory compute_photon')
+              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux::theory compute_photon_to_disk')
 
 
         # Check positivity bound
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 11f07307f0..dc091f8a27 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -286,7 +286,7 @@ def load_photon(self):
       self.integral = np.stack(integral, axis=-1)
       return
     
-def compute_photon(theoryid, fiatlux):
+def compute_photon_to_disk(theoryid, fiatlux, maxcores):
     """Function to compute the photon PDF set.""" 
     luxset = fiatlux['luxset'].load()
     force_compute = fiatlux.get('compute_in_setupfit', False)
@@ -294,8 +294,8 @@ def compute_photon(theoryid, fiatlux):
 
     # Return None and avoid pickling issues with the phoyon class.
     def wrapper_fn(replica, theoryid, fiatlux, force_compute):
-        Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=True, force_computation=force_compute)
+        Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=force_compute, force_computation=force_compute)
         return None
     
-    log.info(f"Starting computation of the photon using {effective_n_jobs(-1)} effective cores...")
-    _ = Parallel(n_jobs=-1)(delayed(wrapper_fn)(replica=replica, theoryid=theoryid, fiatlux=fiatlux, force_compute=force_compute) for replica in replicas)
\ No newline at end of file
+    log.info(f"Starting computation of the photon using {maxcores} effective cores...")
+    _ = Parallel(n_jobs=maxcores)(delayed(wrapper_fn)(replica=replica, theoryid=theoryid, fiatlux=fiatlux, force_compute=force_compute) for replica in replicas)
\ No newline at end of file

From 1f628b7c840b75177af87396ae89b8c1e2f3c6a9 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 18 Nov 2025 13:25:28 +0000
Subject: [PATCH 26/30] Change log message depending on flag

---
 n3fit/src/n3fit/scripts/vp_setupfit.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index 2aee06f900..6b00419e3d 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -187,17 +187,21 @@ def from_yaml(cls, o, *args, **kwargs):
         if fiatlux is not None:
             luxset = fiatlux['luxset']
             theoryid = file_content['theory']['theoryid']
-            force_compute = fiatlux.get('compute_in_setupfit', False)
+            compute_in_setupfit = fiatlux.get('compute_in_setupfit', False)
             try:
                 _ = loader.check_photonQED(theoryid, luxset)
                 log.info(f"Photon QED set found for {theoryid} with luxset {luxset}.")
             except PhotonQEDNotFound:
-                log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Consider "\
-                            "using `compute_in_setupfit` in the runcard to compute all replicas of the photon "
-                            "in vp-setupfit. Otherwise n3fit " \
-                            "will take care of the photon computation. May impact performance.")
-
-            if force_compute:
+                if compute_in_setupfit:
+                  log.warning(f"Photon QED set for {theoryid} with luxset {luxset} not found. " \
+                              "It will be computed in vp-setupfit.")
+                else:
+                  log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Consider "\
+                              "using `compute_in_setupfit` in the runcard to compute all replicas of the photon "
+                              "in vp-setupfit. Otherwise n3fit " \
+                              "will take care of the photon computation. May impact performance.")
+
+            if compute_in_setupfit:
               log.info("Forcing photon computation with FiatLux.")
               # Since the photon will be computed, check that the luxset and additional_errors exist
               SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset_exists')

From 0a8508bb472724e860d8e759688c7e1ec399f88c Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 25 Nov 2025 17:37:53 +0000
Subject: [PATCH 27/30] Refactoring Photon + one EKO for all replicas

---
 validphys2/src/validphys/photon/compute.py | 335 +++++++++++----------
 1 file changed, 173 insertions(+), 162 deletions(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index dc091f8a27..57ceb616b0 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -2,7 +2,8 @@
 
 import logging
 import tempfile
-from joblib import Parallel, delayed, effective_n_jobs
+from joblib import Parallel, delayed
+from functools import lru_cache
 
 import numpy as np
 from scipy.integrate import trapezoid
@@ -14,6 +15,7 @@
 from n3fit.io.writer import XGRID
 from validphys.n3fit_data import replica_luxseed
 from validphys.loader import Loader, PhotonQEDNotFound
+from validphys.core import FKTableSpec
 
 from . import structure_functions as sf
 from .alpha import Alpha
@@ -23,6 +25,7 @@
 # not the complete fiatlux runcard since some parameters are set in the code
 FIATLUX_DEFAULT = {
     "apfel": False,
+    "eps_base": 1e-5,
     "eps_rel": 1e-1,  # extra precision on any single integration.
     "mum_proton": 2.792847356,  # proton magnetic moment, from
     # http://pdglive.lbl.gov/DataBlock.action?node=S016MM which itself
@@ -51,113 +54,97 @@
 
 
 class Photon:
-    """Photon class computing the photon array with the LuxQED approach.
-    
-    Parameters
-    ----------
-    theoryid : validphys.core.TheoryIDSpec
-        TheoryIDSpec object describing the theory to be used as
-        specified in the runcard.
-    lux_params : dict
-        Dictionary containing the LuxQED parameters as specified
-        in the runcard.
-    replica_list: list[int], optional
-        List of replica ids to be computed. If None, all replicas
-        will be computed based on the luxqed pdf set.
-    """
-    def __init__(self, theoryid, lux_params, replicas=None, save_to_disk=False, force_computation=False):
+    def __init__(self, theoryid, lux_params, replicas: list[int], save_to_disk=False, force_computation=False, q_in=100.0):
         self.theoryid = theoryid
         self.lux_params = lux_params
         self.replicas = replicas
         self.save_to_disk = save_to_disk
+        self.q_in = q_in
+        self.luxpdfset = lux_params["luxset"].load()
+        self.luxpdfset_members = self.luxpdfset.n_members - 1 #
 
-        fiatlux_runcard = FIATLUX_DEFAULT
-        # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
-        # This may be changed in the future in favor of a bool em_running in the runcard
-        fiatlux_runcard["qed_running"] = True
-        fiatlux_runcard["mproton"] = float(theoryid.get_description()["MP"])
-
-        # precision on final integration of double integral
-        if "eps_base" in lux_params:
-            fiatlux_runcard["eps_base"] = lux_params["eps_base"]
-            log.warning(f"Using fiatlux parameter eps_base from runcard")
+        # Compute or load photon_qin
+        if force_computation:
+          photon_in = self._compute_photon_set()
         else:
-            fiatlux_runcard["eps_base"] = 1e-5
-            log.info(f"Using default value for fiatlux parameter eps_base")
+          try:
+            photon_in = self._load_photon()
+          except PhotonQEDNotFound:
+            photon_in = self._compute_photon_set()
+        self.photon_qin = photon_in
 
-        self.fiatlux_runcard = fiatlux_runcard
-        # Metadata for the photon ste
-        self.luxpdfset = lux_params["luxset"].load()
-        self.additional_errors = lux_params["additional_errors"]
-        self.luxseed = lux_params["luxseed"]
-        self.luxpdfset_members = self.luxpdfset.n_members - 1 # Remove replica 0
+    def _load_photon(self):
+      """Load the photon resource using the Loader class."""
+      loader = Loader()
+      path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
+      log.info(f"Loading photon QED set from {path_to_photon}")
 
-        if force_computation:
-            self.compute_photon_set()
-            return
-            
-        try:
-          self.load_photon()
-        except PhotonQEDNotFound:
-          log.info(f"Photon set for theory ID {self.theoryid.id} and luxset {self.luxpdfset._name} not found. Computing it now...")
-          self.compute_photon_set()
-
-    def compute_photon_set(self):
-        """Compute the photon set for the desired replicas."""
-        # load fiatlux
-        try:
+      # Load the needed replicas
+      photon_qin_array = []
+      for replica in self.replicas:
+          # As input replica for the photon computation we take the MOD of the luxset_members to
+          # avoid failing due to limited number of replicas in the luxset
+          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+          photon_qin = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_qin"]
+          photon_qin_array.append(photon_qin)
+
+      photon_qin_array = np.stack(photon_qin_array, axis=0)
+      return photon_qin_array
+    
+    def _compute_photon_set(self):
+       try:
           import fiatlux
-        except ModuleNotFoundError as e:
+       except ModuleNotFoundError as e:              
           log.error("fiatlux not found, please install fiatlux")
           raise ModuleNotFoundError("Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`") from e
-        
-        theory = self.theoryid.get_description()
-
-        if theory["PTO"] > 0:
+       
+       #########
+       replicas = self.replicas
+       luxset = self.lux_params['luxset'].load()
+       fiatlux_runcard = self._setup_fiatlux_runcard()
+       #########
+
+       # Set fiatlux parameters
+       theory = self.theoryid.get_description()
+       mb_thr = theory["kbThr"] * theory["mb"]
+       mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
+       f2lo_func = lambda member: sf.F2LO(member, theory)
+       if theory["PTO"] > 0:
             path_to_F2 = self.theoryid.path / "fastkernel/FIATLUX_DIS_F2.pineappl.lz4"
             path_to_FL = self.theoryid.path / "fastkernel/FIATLUX_DIS_FL.pineappl.lz4"
-
-        self.path_to_eko_photon = self.theoryid.path / "eko_photon.tar"
-        with EKO.read(self.path_to_eko_photon) as eko:
-            self.q_in = np.sqrt(eko.mu20)
-
-        # set fiatlux
-        mb_thr = theory["kbThr"] * theory["mb"]
-        mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
-
-        interpolator = []
-        integral = []
-
-        # If saving to disk, create the directory
-        if self.save_to_disk:
+            f2_func = lambda member: sf.InterpStructureFunction(path_to_F2, member)
+            fl_func = lambda member: sf.InterpStructureFunction(path_to_FL, member)
+            
+            # Use central replica to check q2_max
+            f2 = f2_func(luxset.central_member)
+            fl = fl_func(luxset.central_member)
+            if not np.isclose(f2.q2_max, fl.q2_max):
+                log.error(
+                    "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
+                )
+                raise ValueError("FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max")
+            fiatlux_runcard["q2_max"] = float(f2.q2_max)
+       else:
+            f2_func = f2lo_func
+            fl_func = lambda _: sf.FLLO()
+
+       alpha = Alpha(theory, fiatlux_runcard["q2_max"])
+
+       if self.save_to_disk:
           loader = Loader()
-          path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{self.luxpdfset._name}"
+          path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{luxset._name}"
           path_to_photon.mkdir(parents=True, exist_ok=True)
 
-        for replica in self.replicas:
-          # As input replica for the photon computation we take the MOD of the luxset_members to
-          # avoid failing due to limited number of replicas in the luxset
+       photon_qin_array = []
+       for replica in replicas:
+          # Avoid failing due to limited number of replicas in the luxset
           photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+          f2lo = f2lo_func(luxset.members[photonreplica])
+          f2 = f2_func(luxset.members[photonreplica])
+          fl = fl_func(luxset.members[photonreplica])
 
-          f2lo = sf.F2LO(self.luxpdfset.members[photonreplica], theory)
-
-          if theory["PTO"] > 0:
-              f2 = sf.InterpStructureFunction(path_to_F2, self.luxpdfset.members[photonreplica])
-              fl = sf.InterpStructureFunction(path_to_FL, self.luxpdfset.members[photonreplica])
-              if not np.isclose(f2.q2_max, fl.q2_max):
-                  log.error(
-                      "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
-                  )
-              self.fiatlux_runcard["q2_max"] = float(f2.q2_max)
-          else:
-              f2 = f2lo
-              fl = sf.FLLO()
-              # using a default value for q2_max
-              self.fiatlux_runcard["q2_max"] = 1e8
-
-          alpha = Alpha(theory, self.fiatlux_runcard["q2_max"])
           with tempfile.NamedTemporaryFile(mode="w") as tmp:
-              yaml.dump(self.fiatlux_runcard, tmp)
+              yaml.dump(fiatlux_runcard, tmp)
               lux = fiatlux.FiatLux(tmp.name)
 
           # we have a dict but fiatlux wants a yaml file
@@ -169,17 +156,87 @@ def compute_photon_set(self):
 
           photon_qin = np.array(
             [lux.EvaluatePhoton(x, self.q_in**2).total for x in XGRID]
-        )
-          photon_qin += self.generate_errors(replica)
+          )
+          photon_qin += self._generate_errors(replica)
 
           # fiatlux computes x * gamma(x)
           photon_qin /= XGRID
 
-          # Load eko and reshape it
-          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
-              
+          if self.save_to_disk:
+              np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz", photon_qin=photon_qin)
+              log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
+
+          photon_qin_array.append(photon_qin)
+
+       photon_qin_array = np.stack(photon_qin_array, axis=0)
+       return photon_qin_array
+       
+    def _setup_fiatlux_runcard(self):
+      fiatlux_runcard = FIATLUX_DEFAULT
+
+      # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
+      # This may be changed in the future in favor of a bool em_running in the runcard
+      fiatlux_runcard["qed_running"] = True
+      fiatlux_runcard["mproton"] = float(self.theoryid.get_description()["MP"])
+
+      # precision on final integration of double integral
+      if "eps_base" in self.lux_params:
+          fiatlux_runcard["eps_base"] = self.lux_params["eps_base"]
+          log.warning(f"Using fiatlux parameter eps_base from runcard")
+      else:
+          log.info(f"Using default value for fiatlux parameter eps_base = {fiatlux_runcard['eps_base']}")
+      return fiatlux_runcard
+    
+    @property
+    def error_matrix(self):
+        """Generate error matrix to be used in generate_errors."""
+        if "additional_errors" not in self.lux_params:
+            return None
+        extra_set = self.lux_params["additional_errors"].load()
+        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):
+            tmp = np.array(extra_set.members[idx_member].xfxQ(22, XGRID, qs))
+            res.append(tmp - res_central)
+        # first index must be x, while second one must be replica index
+        return np.stack(res, axis=1)
+    
+    def _generate_errors(self, replica):
+        """
+        Generate LUX additional errors according to the procedure
+        described in sec. 2.5 of https://arxiv.org/pdf/1712.07053.pdf
+        """
+        if self.error_matrix is None:
+            return np.zeros_like(XGRID)
+        log.info(f"Generating photon additional errors")
+        luxseed = self.lux_params.get("luxseed", None)
+        seed = replica_luxseed(replica, luxseed)
+        rng = np.random.default_rng(seed=seed)
+        u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False)
+        errors = u @ (s * rng.normal(size=7))
+        return errors
+
+    @lru_cache(maxsize=None)
+    def _evolve(self):
+        """Perform the EKOs to evolve the photon from q_in to Q0."""
+        log.info(f"Evolving photon from q_in={self.q_in} GeV to Q0 using EKO...")
+        photon_qin_array = self.photon_qin
+        interpolator = []
+        integral = []
+
+        path_to_eko_photon = self.theoryid.path / "eko_photon.tar"
+        with EKO.read(path_to_eko_photon) as eko_photon:
+            # Check that qin mathces with the one in the EKO
+            if not np.isclose(self.q_in, np.sqrt(eko_photon.mu20)):
+                log.error(f"Photon q_in {self.q_in} does not match the one in the EKO {np.sqrt(eko_photon.mu20)}")
+                raise ValueError("Photon q_in does not match the one in the EKO")
+
+            # TODO : if the eko has not the correct grid we have to reshape it
+            # it has to be done inside vp-setupfit
+            for replica in self.replicas:
+              photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+            
               # NB: the eko should contain a single operator
               for _, elem in eko_photon.items():
                   eko_op = elem.operator
@@ -187,7 +244,7 @@ def compute_photon_set(self):
                   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
+                          pdfs_init[j] = photon_qin_array[replica-1]
                           ph_id = j
                       elif pid not in self.luxpdfset.flavors:
                           continue
@@ -198,20 +255,26 @@ def compute_photon_set(self):
 
                   pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
 
-          photon_Q0 = pdfs_final[ph_id]
-          photon_array = XGRID * photon_Q0
-
-          if self.save_to_disk:
-              np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz",photon_array=photon_array)
-              log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
-
-          interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
-          integral.append(trapezoid(photon_array, XGRID))
+              photon_Q0 = pdfs_final[ph_id]
+              photon_array = XGRID * photon_Q0
+              interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
+              integral.append(trapezoid(photon_array, XGRID))
 
         integral = np.stack(integral, axis=-1)
-        self.integral = integral
-        self.interpolator = interpolator
+        return integral, interpolator
+    
+    @property
+    def integral(self):
+        """Return the integral values."""
+        integral, _ = self._evolve()
+        return integral
 
+    @property
+    def interpolator(self):
+        """Return the interpolator functions."""
+        _, interpolator = self._evolve()
+        return interpolator
+    
     def __call__(self, xgrid):
         """
         Compute the photon interpolating the values of self.photon_array.
@@ -233,58 +296,6 @@ def __call__(self, xgrid):
             ],
             axis=1,
         )
-
-    @property
-    def error_matrix(self):
-        """Generate error matrix to be used in generate_errors."""
-        if not self.additional_errors:
-            return None
-        extra_set = self.additional_errors.load()
-        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):
-            tmp = np.array(extra_set.members[idx_member].xfxQ(22, XGRID, qs))
-            res.append(tmp - res_central)
-        # first index must be x, while second one must be replica index
-        return np.stack(res, axis=1)
-
-    def generate_errors(self, replica):
-        """
-        Generate LUX additional errors according to the procedure
-        described in sec. 2.5 of https://arxiv.org/pdf/1712.07053.pdf
-        """
-        if self.error_matrix is None:
-            return np.zeros_like(XGRID)
-        log.info(f"Generating photon additional errors")
-        seed = replica_luxseed(replica, self.luxseed)
-        rng = np.random.default_rng(seed=seed)
-        u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False)
-        errors = u @ (s * rng.normal(size=7))
-        return errors
-    
-    def load_photon(self):
-      """Load the photon resource using the Loader class."""
-      loader = Loader()
-      path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
-      log.info(f"Loading photon QED set from {path_to_photon}")
-
-      interpolator = []
-      integral = []
-
-      # Load the needed replicas
-      for replica in self.replicas:
-          # As input replica for the photon computation we take the MOD of the luxset_members to
-          # avoid failing due to limited number of replicas in the luxset
-          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
-
-          photon_array = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_array"]
-          interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
-          integral.append(trapezoid(photon_array, XGRID))
-      
-      self.interpolator = interpolator
-      self.integral = np.stack(integral, axis=-1)
-      return
     
 def compute_photon_to_disk(theoryid, fiatlux, maxcores):
     """Function to compute the photon PDF set.""" 
@@ -292,9 +303,9 @@ def compute_photon_to_disk(theoryid, fiatlux, maxcores):
     force_compute = fiatlux.get('compute_in_setupfit', False)
     replicas = list(range(1, luxset.n_members))
 
-    # Return None and avoid pickling issues with the phoyon class.
+    # Return None and avoid pickling issues with the photon class.
     def wrapper_fn(replica, theoryid, fiatlux, force_compute):
-        Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=force_compute, force_computation=force_compute)
+        _ = Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=force_compute, force_computation=force_compute)
         return None
     
     log.info(f"Starting computation of the photon using {maxcores} effective cores...")

From d0a21cf5eedd2d9a597ea7eaee3efab5905a2284 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Tue, 25 Nov 2025 17:43:57 +0000
Subject: [PATCH 28/30] Pre-commit

---
 validphys2/src/validphys/photon/compute.py | 315 ++++++++++++---------
 1 file changed, 175 insertions(+), 140 deletions(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 57ceb616b0..47af4fc602 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -1,10 +1,10 @@
 """Module that calls fiatlux to add the photon PDF."""
 
+from functools import lru_cache
 import logging
 import tempfile
-from joblib import Parallel, delayed
-from functools import lru_cache
 
+from joblib import Parallel, delayed
 import numpy as np
 from scipy.integrate import trapezoid
 from scipy.interpolate import interp1d
@@ -13,9 +13,9 @@
 from eko import basis_rotation
 from eko.io import EKO
 from n3fit.io.writer import XGRID
-from validphys.n3fit_data import replica_luxseed
-from validphys.loader import Loader, PhotonQEDNotFound
 from validphys.core import FKTableSpec
+from validphys.loader import Loader, PhotonQEDNotFound
+from validphys.n3fit_data import replica_luxseed
 
 from . import structure_functions as sf
 from .alpha import Alpha
@@ -54,67 +54,77 @@
 
 
 class Photon:
-    def __init__(self, theoryid, lux_params, replicas: list[int], save_to_disk=False, force_computation=False, q_in=100.0):
+    def __init__(
+        self,
+        theoryid,
+        lux_params,
+        replicas: list[int],
+        save_to_disk=False,
+        force_computation=False,
+        q_in=100.0,
+    ):
         self.theoryid = theoryid
         self.lux_params = lux_params
         self.replicas = replicas
         self.save_to_disk = save_to_disk
         self.q_in = q_in
         self.luxpdfset = lux_params["luxset"].load()
-        self.luxpdfset_members = self.luxpdfset.n_members - 1 #
+        self.luxpdfset_members = self.luxpdfset.n_members - 1  #
 
         # Compute or load photon_qin
         if force_computation:
-          photon_in = self._compute_photon_set()
-        else:
-          try:
-            photon_in = self._load_photon()
-          except PhotonQEDNotFound:
             photon_in = self._compute_photon_set()
+        else:
+            try:
+                photon_in = self._load_photon()
+            except PhotonQEDNotFound:
+                photon_in = self._compute_photon_set()
         self.photon_qin = photon_in
 
     def _load_photon(self):
-      """Load the photon resource using the Loader class."""
-      loader = Loader()
-      path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
-      log.info(f"Loading photon QED set from {path_to_photon}")
-
-      # Load the needed replicas
-      photon_qin_array = []
-      for replica in self.replicas:
-          # As input replica for the photon computation we take the MOD of the luxset_members to
-          # avoid failing due to limited number of replicas in the luxset
-          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
-          photon_qin = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_qin"]
-          photon_qin_array.append(photon_qin)
-
-      photon_qin_array = np.stack(photon_qin_array, axis=0)
-      return photon_qin_array
-    
+        """Load the photon resource using the Loader class."""
+        loader = Loader()
+        path_to_photon = loader.check_photonQED(self.theoryid.id, self.luxpdfset._name)
+        log.info(f"Loading photon QED set from {path_to_photon}")
+
+        # Load the needed replicas
+        photon_qin_array = []
+        for replica in self.replicas:
+            # As input replica for the photon computation we take the MOD of the luxset_members to
+            # avoid failing due to limited number of replicas in the luxset
+            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+            photon_qin = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_qin"]
+            photon_qin_array.append(photon_qin)
+
+        photon_qin_array = np.stack(photon_qin_array, axis=0)
+        return photon_qin_array
+
     def _compute_photon_set(self):
-       try:
-          import fiatlux
-       except ModuleNotFoundError as e:              
-          log.error("fiatlux not found, please install fiatlux")
-          raise ModuleNotFoundError("Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`") from e
-       
-       #########
-       replicas = self.replicas
-       luxset = self.lux_params['luxset'].load()
-       fiatlux_runcard = self._setup_fiatlux_runcard()
-       #########
-
-       # Set fiatlux parameters
-       theory = self.theoryid.get_description()
-       mb_thr = theory["kbThr"] * theory["mb"]
-       mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
-       f2lo_func = lambda member: sf.F2LO(member, theory)
-       if theory["PTO"] > 0:
+        try:
+            import fiatlux
+        except ModuleNotFoundError as e:
+            log.error("fiatlux not found, please install fiatlux")
+            raise ModuleNotFoundError(
+                "Please install fiatlux: `pip install nnpdf[qed]` or `pip install fiatlux`"
+            ) from e
+
+        #########
+        replicas = self.replicas
+        luxset = self.lux_params['luxset'].load()
+        fiatlux_runcard = self._setup_fiatlux_runcard()
+        #########
+
+        # Set fiatlux parameters
+        theory = self.theoryid.get_description()
+        mb_thr = theory["kbThr"] * theory["mb"]
+        mt_thr = theory["ktThr"] * theory["mt"] if theory["MaxNfPdf"] == 6 else 1e100
+        f2lo_func = lambda member: sf.F2LO(member, theory)
+        if theory["PTO"] > 0:
             path_to_F2 = self.theoryid.path / "fastkernel/FIATLUX_DIS_F2.pineappl.lz4"
             path_to_FL = self.theoryid.path / "fastkernel/FIATLUX_DIS_FL.pineappl.lz4"
             f2_func = lambda member: sf.InterpStructureFunction(path_to_F2, member)
             fl_func = lambda member: sf.InterpStructureFunction(path_to_FL, member)
-            
+
             # Use central replica to check q2_max
             f2 = f2_func(luxset.central_member)
             fl = fl_func(luxset.central_member)
@@ -122,71 +132,77 @@ def _compute_photon_set(self):
                 log.error(
                     "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
                 )
-                raise ValueError("FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max")
+                raise ValueError(
+                    "FKtables for FIATLUX_DIS_F2 and FIATLUX_DIS_FL have two different q2_max"
+                )
             fiatlux_runcard["q2_max"] = float(f2.q2_max)
-       else:
+        else:
             f2_func = f2lo_func
             fl_func = lambda _: sf.FLLO()
 
-       alpha = Alpha(theory, fiatlux_runcard["q2_max"])
-
-       if self.save_to_disk:
-          loader = Loader()
-          path_to_photon = loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{luxset._name}"
-          path_to_photon.mkdir(parents=True, exist_ok=True)
-
-       photon_qin_array = []
-       for replica in replicas:
-          # Avoid failing due to limited number of replicas in the luxset
-          photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
-          f2lo = f2lo_func(luxset.members[photonreplica])
-          f2 = f2_func(luxset.members[photonreplica])
-          fl = fl_func(luxset.members[photonreplica])
-
-          with tempfile.NamedTemporaryFile(mode="w") as tmp:
-              yaml.dump(fiatlux_runcard, tmp)
-              lux = fiatlux.FiatLux(tmp.name)
-
-          # we have a dict but fiatlux wants a yaml file
-          # TODO : once that fiatlux will allow dictionaries
-          # pass directly fiatlux_runcard
-          lux.PlugAlphaQED(alpha.alpha_em, alpha.qref)
-          lux.InsertInelasticSplitQ([mb_thr, mt_thr])
-          lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
-
-          photon_qin = np.array(
-            [lux.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
-
-          if self.save_to_disk:
-              np.savez_compressed(path_to_photon / f"replica_{photonreplica}.npz", photon_qin=photon_qin)
-              log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
-
-          photon_qin_array.append(photon_qin)
-
-       photon_qin_array = np.stack(photon_qin_array, axis=0)
-       return photon_qin_array
-       
+        alpha = Alpha(theory, fiatlux_runcard["q2_max"])
+
+        if self.save_to_disk:
+            loader = Loader()
+            path_to_photon = (
+                loader._photons_qed_path / f"photon_theoryID_{self.theoryid.id}_fit_{luxset._name}"
+            )
+            path_to_photon.mkdir(parents=True, exist_ok=True)
+
+        photon_qin_array = []
+        for replica in replicas:
+            # Avoid failing due to limited number of replicas in the luxset
+            photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+            f2lo = f2lo_func(luxset.members[photonreplica])
+            f2 = f2_func(luxset.members[photonreplica])
+            fl = fl_func(luxset.members[photonreplica])
+
+            with tempfile.NamedTemporaryFile(mode="w") as tmp:
+                yaml.dump(fiatlux_runcard, tmp)
+                lux = fiatlux.FiatLux(tmp.name)
+
+            # we have a dict but fiatlux wants a yaml file
+            # TODO : once that fiatlux will allow dictionaries
+            # pass directly fiatlux_runcard
+            lux.PlugAlphaQED(alpha.alpha_em, alpha.qref)
+            lux.InsertInelasticSplitQ([mb_thr, mt_thr])
+            lux.PlugStructureFunctions(f2.fxq, fl.fxq, f2lo.fxq)
+
+            photon_qin = np.array([lux.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
+
+            if self.save_to_disk:
+                np.savez_compressed(
+                    path_to_photon / f"replica_{photonreplica}.npz", photon_qin=photon_qin
+                )
+                log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
+
+            photon_qin_array.append(photon_qin)
+
+        photon_qin_array = np.stack(photon_qin_array, axis=0)
+        return photon_qin_array
+
     def _setup_fiatlux_runcard(self):
-      fiatlux_runcard = FIATLUX_DEFAULT
-
-      # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
-      # This may be changed in the future in favor of a bool em_running in the runcard
-      fiatlux_runcard["qed_running"] = True
-      fiatlux_runcard["mproton"] = float(self.theoryid.get_description()["MP"])
-
-      # precision on final integration of double integral
-      if "eps_base" in self.lux_params:
-          fiatlux_runcard["eps_base"] = self.lux_params["eps_base"]
-          log.warning(f"Using fiatlux parameter eps_base from runcard")
-      else:
-          log.info(f"Using default value for fiatlux parameter eps_base = {fiatlux_runcard['eps_base']}")
-      return fiatlux_runcard
-    
+        fiatlux_runcard = FIATLUX_DEFAULT
+
+        # TODO: for the time being, Qedref=Qref and so alphaem running will always trigger
+        # This may be changed in the future in favor of a bool em_running in the runcard
+        fiatlux_runcard["qed_running"] = True
+        fiatlux_runcard["mproton"] = float(self.theoryid.get_description()["MP"])
+
+        # precision on final integration of double integral
+        if "eps_base" in self.lux_params:
+            fiatlux_runcard["eps_base"] = self.lux_params["eps_base"]
+            log.warning(f"Using fiatlux parameter eps_base from runcard")
+        else:
+            log.info(
+                f"Using default value for fiatlux parameter eps_base = {fiatlux_runcard['eps_base']}"
+            )
+        return fiatlux_runcard
+
     @property
     def error_matrix(self):
         """Generate error matrix to be used in generate_errors."""
@@ -201,7 +217,7 @@ def error_matrix(self):
             res.append(tmp - res_central)
         # first index must be x, while second one must be replica index
         return np.stack(res, axis=1)
-    
+
     def _generate_errors(self, replica):
         """
         Generate LUX additional errors according to the procedure
@@ -229,40 +245,47 @@ def _evolve(self):
         with EKO.read(path_to_eko_photon) as eko_photon:
             # Check that qin mathces with the one in the EKO
             if not np.isclose(self.q_in, np.sqrt(eko_photon.mu20)):
-                log.error(f"Photon q_in {self.q_in} does not match the one in the EKO {np.sqrt(eko_photon.mu20)}")
+                log.error(
+                    f"Photon q_in {self.q_in} does not match the one in the EKO {np.sqrt(eko_photon.mu20)}"
+                )
                 raise ValueError("Photon q_in does not match the one in the EKO")
 
             # TODO : if the eko has not the correct grid we have to reshape it
             # it has to be done inside vp-setupfit
             for replica in self.replicas:
-              photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
-            
-              # 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_array[replica-1]
-                          ph_id = j
-                      elif pid not in self.luxpdfset.flavors:
-                          continue
-                      else:
-                          pdfs_init[j] = np.array(
-                              [self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x for x in XGRID]
-                          )
-
-                  pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
-
-              photon_Q0 = pdfs_final[ph_id]
-              photon_array = XGRID * photon_Q0
-              interpolator.append(interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic"))
-              integral.append(trapezoid(photon_array, XGRID))
+                photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
+
+                # 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_array[replica - 1]
+                            ph_id = j
+                        elif pid not in self.luxpdfset.flavors:
+                            continue
+                        else:
+                            pdfs_init[j] = np.array(
+                                [
+                                    self.luxpdfset.xfxQ(x, self.q_in, photonreplica, pid) / x
+                                    for x in XGRID
+                                ]
+                            )
+
+                    pdfs_final = np.einsum("ajbk,bk", eko_op, pdfs_init)
+
+                photon_Q0 = pdfs_final[ph_id]
+                photon_array = XGRID * photon_Q0
+                interpolator.append(
+                    interp1d(XGRID, photon_array, fill_value="extrapolate", kind="cubic")
+                )
+                integral.append(trapezoid(photon_array, XGRID))
 
         integral = np.stack(integral, axis=-1)
         return integral, interpolator
-    
+
     @property
     def integral(self):
         """Return the integral values."""
@@ -274,7 +297,7 @@ def interpolator(self):
         """Return the interpolator functions."""
         _, interpolator = self._evolve()
         return interpolator
-    
+
     def __call__(self, xgrid):
         """
         Compute the photon interpolating the values of self.photon_array.
@@ -296,17 +319,29 @@ def __call__(self, xgrid):
             ],
             axis=1,
         )
-    
+
+
 def compute_photon_to_disk(theoryid, fiatlux, maxcores):
-    """Function to compute the photon PDF set.""" 
+    """Function to compute the photon PDF set."""
     luxset = fiatlux['luxset'].load()
     force_compute = fiatlux.get('compute_in_setupfit', False)
     replicas = list(range(1, luxset.n_members))
 
     # Return None and avoid pickling issues with the photon class.
     def wrapper_fn(replica, theoryid, fiatlux, force_compute):
-        _ = Photon(theoryid, fiatlux, replicas=[replica], save_to_disk=force_compute, force_computation=force_compute)
+        _ = Photon(
+            theoryid,
+            fiatlux,
+            replicas=[replica],
+            save_to_disk=force_compute,
+            force_computation=force_compute,
+        )
         return None
-    
+
     log.info(f"Starting computation of the photon using {maxcores} effective cores...")
-    _ = Parallel(n_jobs=maxcores)(delayed(wrapper_fn)(replica=replica, theoryid=theoryid, fiatlux=fiatlux, force_compute=force_compute) for replica in replicas)
\ No newline at end of file
+    _ = Parallel(n_jobs=maxcores)(
+        delayed(wrapper_fn)(
+            replica=replica, theoryid=theoryid, fiatlux=fiatlux, force_compute=force_compute
+        )
+        for replica in replicas
+    )

From b6770f4b00a2a62ecad073bb240318abaaaf5bce Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Wed, 26 Nov 2025 15:37:07 +0000
Subject: [PATCH 29/30] Storing photon replicas in a dictionary

---
 validphys2/src/validphys/photon/compute.py    | 20 +++++++++----------
 .../validphys/tests/photon/test_compute.py    | 14 +++++++++----
 2 files changed, 19 insertions(+), 15 deletions(-)

diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 47af4fc602..d3effa8ec1 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -69,7 +69,8 @@ def __init__(
         self.save_to_disk = save_to_disk
         self.q_in = q_in
         self.luxpdfset = lux_params["luxset"].load()
-        self.luxpdfset_members = self.luxpdfset.n_members - 1  #
+        self.luxpdfset_members = self.luxpdfset.n_members - 1
+        self.path_to_eko_photon = self.theoryid.path / "eko_photon.tar"
 
         # Compute or load photon_qin
         if force_computation:
@@ -88,15 +89,14 @@ def _load_photon(self):
         log.info(f"Loading photon QED set from {path_to_photon}")
 
         # Load the needed replicas
-        photon_qin_array = []
+        photon_qin_array = {}
         for replica in self.replicas:
             # As input replica for the photon computation we take the MOD of the luxset_members to
             # avoid failing due to limited number of replicas in the luxset
             photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
             photon_qin = np.load(path_to_photon / f"replica_{photonreplica}.npz")["photon_qin"]
-            photon_qin_array.append(photon_qin)
+            photon_qin_array[replica] = photon_qin
 
-        photon_qin_array = np.stack(photon_qin_array, axis=0)
         return photon_qin_array
 
     def _compute_photon_set(self):
@@ -149,7 +149,7 @@ def _compute_photon_set(self):
             )
             path_to_photon.mkdir(parents=True, exist_ok=True)
 
-        photon_qin_array = []
+        photon_qin_array = {}
         for replica in replicas:
             # Avoid failing due to limited number of replicas in the luxset
             photonreplica = (replica % self.luxpdfset_members) or self.luxpdfset_members
@@ -180,9 +180,8 @@ def _compute_photon_set(self):
                 )
                 log.info(f"Saved photon replica {photonreplica} to {path_to_photon}")
 
-            photon_qin_array.append(photon_qin)
+            photon_qin_array[replica] = photon_qin
 
-        photon_qin_array = np.stack(photon_qin_array, axis=0)
         return photon_qin_array
 
     def _setup_fiatlux_runcard(self):
@@ -206,7 +205,7 @@ def _setup_fiatlux_runcard(self):
     @property
     def error_matrix(self):
         """Generate error matrix to be used in generate_errors."""
-        if "additional_errors" not in self.lux_params:
+        if not self.lux_params["additional_errors"]:
             return None
         extra_set = self.lux_params["additional_errors"].load()
         qs = [self.q_in] * len(XGRID)
@@ -241,8 +240,7 @@ def _evolve(self):
         interpolator = []
         integral = []
 
-        path_to_eko_photon = self.theoryid.path / "eko_photon.tar"
-        with EKO.read(path_to_eko_photon) as eko_photon:
+        with EKO.read(self.path_to_eko_photon) as eko_photon:
             # Check that qin mathces with the one in the EKO
             if not np.isclose(self.q_in, np.sqrt(eko_photon.mu20)):
                 log.error(
@@ -262,7 +260,7 @@ def _evolve(self):
                     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_array[replica - 1]
+                            pdfs_init[j] = photon_qin_array[replica]
                             ph_id = j
                         elif pid not in self.luxpdfset.flavors:
                             continue
diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py
index 274f6673fa..db93bda3cb 100644
--- a/validphys2/src/validphys/tests/photon/test_compute.py
+++ b/validphys2/src/validphys/tests/photon/test_compute.py
@@ -37,12 +37,16 @@ def test_parameters_init():
     # we are not testing the photon here so we make it faster
     fiatlux_runcard['eps_base'] = 1e-1
 
-    photon = Photon(test_theory, fiatlux_runcard, [1, 2, 3])
+    photon = Photon(
+        test_theory, fiatlux_runcard, [1, 2, 3], save_to_disk=False, force_computation=True
+    )
 
     np.testing.assert_equal(photon.replicas, [1, 2, 3])
     np.testing.assert_equal(photon.luxpdfset._name, fiatlux_runcard["luxset"].name)
-    np.testing.assert_equal(photon.additional_errors.name, "LUXqed17_plus_PDF4LHC15_nnlo_100")
-    np.testing.assert_equal(photon.luxseed, fiatlux_runcard["luxseed"])
+    np.testing.assert_equal(
+        photon.lux_params["additional_errors"].name, "LUXqed17_plus_PDF4LHC15_nnlo_100"
+    )
+    np.testing.assert_equal(photon.lux_params["luxseed"], fiatlux_runcard["luxseed"])
     np.testing.assert_equal(photon.path_to_eko_photon, test_theory.path / "eko_photon.tar")
     np.testing.assert_equal(photon.q_in, 100.0)
 
@@ -58,7 +62,9 @@ def test_photon():
     theory = test_theory.get_description()
 
     for replica in [1, 2, 3]:
-        photon = Photon(test_theory, fiatlux_runcard, [replica])
+        photon = Photon(
+            test_theory, fiatlux_runcard, [replica], save_to_disk=False, force_computation=True
+        )
 
         # set up fiatlux
         path_to_F2 = test_theory.path / "fastkernel/FIATLUX_DIS_F2.pineappl.lz4"

From fab60cdbb6fd23dbe31fd0255c49f6d4d800f2d6 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 18 Dec 2025 13:34:00 +0100
Subject: [PATCH 30/30] Upload documentation

---
 doc/sphinx/source/tutorials/run-qed-fit.rst | 234 ++++++++++++++++++--
 doc/sphinx/source/vp/download.rst           |  44 +++-
 doc/sphinx/source/vp/nnprofile.rst          |  20 +-
 n3fit/src/n3fit/checks.py                   |  22 +-
 n3fit/src/n3fit/scripts/vp_setupfit.py      |  39 ++--
 validphys2/src/validphys/loader.py          |  26 ++-
 6 files changed, 316 insertions(+), 69 deletions(-)

diff --git a/doc/sphinx/source/tutorials/run-qed-fit.rst b/doc/sphinx/source/tutorials/run-qed-fit.rst
index 322b466c63..6b1432aee1 100644
--- a/doc/sphinx/source/tutorials/run-qed-fit.rst
+++ b/doc/sphinx/source/tutorials/run-qed-fit.rst
@@ -3,32 +3,232 @@
 ==========================
 How to run a QED fit
 ==========================
+This tutorial describes how to run a QED fit using the LuxQED approach, as
+described in `arXiv:1607.04266 <https://arxiv.org/abs/1607.04266>`_ and
+`arXiv:1708.01256 <https://arxiv.org/abs/1708.01256>`_.
 
-It is possible to perform a QED fit adding the key ``fiatlux`` to the runcard. In this way
-a photon PDF will be generated using the FiatLux public library that implements the LuxQED
-(see :cite:p:`Manohar:2016nzj` and :cite:p:`Manohar:2017eqh`) approach.
-The parameters to be added are the following:
+Setting up the runcard
+----------------------
+
+It is possible to perform a QED fit by adding a ``fiatlux`` block to the
+runcard. The following code snippet shows an example of a QED fit configuration:
 
 .. code-block:: yaml
 
     fiatlux:
-      luxset: NNPDF40_nnlo_as_01180
+      luxset: 251127-jcm-lh-qed-001
       additional_errors: true
       luxseed: 1234567890
 
-``luxset`` is the PDF set used to generate the photon PDF with `FiatLux <https://github.com/scarrazza/fiatlux/>`.
-The code will generate as many photon replicas as the number of replicas contained in the ``luxset``. Therefore, if the user
-tries to generate a replica with ID higher than the maximum ID of the ``luxset``, the code will
-raise an error. Moreover, being the LuxQED approach an iterated prcedure, and given that some replicas
-do not pass the ``postfit`` selection criteria, the user should make sure that the number of replicas in
-the ``luxset`` is high enough such that in the final iteration there will be a number of replicas
-higher than the final replicas desired.
-``additional_errors`` is a parameter that switches on and off the additional errors of the LuxQED  approach,
-while ``luxseed`` is the seed used to generate such errors.
-This errors should be switched on only in the very last iteration of the procedure.
+The parameters contained in the ``fiatlux`` block are:
+
+* ``luxset``
+      The name of the PDF set used to generate the photon PDF with `FiatLux
+      <https://github.com/scarrazza/fiatlux/>`_. The code will use as many
+      photon replicas as the number of replicas contained in the ``luxset``. If
+      the user tries to generate a replica with ID higher than the maximum ID of
+      the ``luxset``, the code will start reusing photon replica from the first.
+      Being the LuxQED approach an iterated procedure, and given that some
+      replicas do not pass the ``postfit`` selection criteria, the user should
+      make sure that the number of replicas in the ``luxset`` is high enough
+      such that in the final iteration there will be a number of replicas higher
+      than the final replicas desired.
+
+* ``additional_errors``
+      Boolean flag to switch on and off the additional errors of the LuxQED approach.
+
+      .. note::
+
+        The ``additional_errors`` flag should be switched on only in the very last
+        iteration of the procedure.
+
+* ``luxseed``
+      The seed used to generate the additional errors of the LuxQED as in ``additional_errors``.
+
+The code uses both the ``fiatlux`` block and the ``theoryid`` specified in the
+runcard to identify the photon PDF set. As explained below, the code searches
+for precomputed photon PDF sets using the pair of ``luxset`` and ``theoryid``
+parameters, first locally and then on the NNPDF server (see
+:ref:`photon-resources` for details). The photon PDF set will be named
+``photon_theoryID_<theoryid>_fit_<luxset>``.
+
+.. _photon-resources:
+
+Running with Photon PDF sets
+-----------------------------
+
+The generation of a photon PDF set can add significant overhead to the fitting
+procedure. Moreover, the same photon PDF set might be used in different fits. To
+minimize the overhead due to the generation of photon PDF sets and avoid
+redundant computations, the code looks for precomputed photon resources either
+locally or on the NNPDF server. If a desired photon PDF set does not exist in
+either of the two locations, it will be computed on the fly and stored locally.
+The following sections describe how to use existing photon PDF sets or generate
+new ones.
+
+Using available Photon PDF sets
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Before running a QED fit, it is strongly advised to prepare the fit using
+``vp-setupfit``, as explained in :ref:`this tutorial <run-n3fit-fit>`. This will
+ensure that all the required resources are available, including the photon PDF.
+The desired photon PDF is specified by the ``luxset`` parameter in the
+``fiatlux`` block and the ``theoryid``, as explained above. The code will first
+look for the photon PDF set in the local directory specified in the
+:ref:`profile file <nnprofile>`. If the set is not found locally, it will try to
+download it from the NNPDF server. The following is an example of running
+``vp-setupfit`` using the ``fiatlux`` block shown above:
+
+.. code-block:: bash
+
+    $ vp-setupfit qed_example_runcard.yml
+    [INFO]: Could not find a resource (photonQED): Could not find Photon QED set /user/envs/nnpdf/share/NNPDF/photons_qed/photon_theoryID_702_fit_251127-jcm-lh-qed-001 in theory: 702. Attempting to download it.
+    [INFO]: Downloading https://data.nnpdf.science/photons/photon_theoryID_702_fit_251127-jcm-lh-qed-001.tar.
+    [==================================================] (100%)
+    [INFO]: Extracting archive to /opt/homebrew/Caskroom/miniconda/base/envs/nnpdf/share/NNPDF/photons_qed
+    [INFO]: Photon QED set found for 702 with luxset 251127-jcm-lh-qed-001.
+    [WARNING]: Using q2min from runcard
+    [WARNING]: Using w2min from runcard
+    Using Keras backend
+    [INFO]: All requirements processed and checked successfully. Executing actions.
+    ...
+
+This will download and extract the photon PDF set in the local
+``photons_qed_path`` specified in the :ref:`profile file <nnprofile>`.
+
+The ``vp-list`` utility tool can be used to list all the available photon PDF
+sets locally and on the NNPDF server. To list the available photon PDF sets,
+just run:
+
+.. code-block:: bash
+
+   $ vp-list photons
+   [INFO]: The following photons are available locally:
+    - theoryID_702_fit_251127-jcm-lh-qed-001
+   [INFO]: The following photons are downloadable:
+    - theoryID_702_fit_251127-jcm-lh-qed-002
+
+In this example, there is one photon PDF set available locally, and one photon
+resource available on the NNPDF server precomputed with theory ID 702 and
+``251127-jcm-lh-qed-002`` as input PDF. The user can manually download a photon
+PDF set using the ``vp-get`` tool as explained :ref:`here <download>`. For
+example:
+
+.. code-block:: bash
+
+  $ vp-get photonQED 702 251127-jcm-lh-qed-002
+  [INFO]: Could not find a resource (photonQED): Could not find Photon QED set /user/envs/nnpdf/share/NNPDF/photons_qed/photon_theoryID_702_fit_251127-jcm-lh-qed-002 in theory: 702. Attempting to download it.
+  [INFO]: Downloading https://data.nnpdf.science/photons/photon_theoryID_702_fit_251127-jcm-lh-qed-002.tar.
+  [==================================================] (100%)
+  [INFO]: Extracting archive to /user/envs/nnpdf/share/NNPDF/photons_qed
+  PosixPath('/user/envs/nnpdf/share/NNPDF/photons_qed/photon_theoryID_702_fit_251127-jcm-lh-qed-002')
 
-Whenever the photon PDF is generated, it will remain constant during the fit and will be considered in the momentum sum rule.
+As in the case of ``vp-setupfit``, this will download and extract the photon PDF
+set in the local ``photons_qed_path`` specified in the :ref:`profile file <nnprofile>`.
+Once the photon PDF set is available locally, the user can proceed to run the
+fit with ``n3fit`` as usual.
 
 .. warning::
 
-   At the moment it is not possible to run QED fits in parallel, as the FiatLux library cannot run in parallel.
+   If ``vp-setupfit`` is not run before ``n3fit``, and the photon PDF set is not
+   available locally, the code will **not** attempt to download it from the server,
+   but will directly proceed to compute it on the fly. See the :ref:`next section <generating-photon-sets>` for
+   more details.
+
+.. _generating-photon-sets:
+Generating new Photon PDF sets
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+If the desired photon PDF set is not available locally nor on the NNPDF server, the
+code will generate the required photon PDF set replicas on the fly using `FiatLux <https://github.com/scarrazza/fiatlux/>`_.
+This can be done either during the ``vp-setupfit`` step, which precomputes all photon
+replicas before starting the fit, or during the fitting step with ``n3fit``, which
+generates the photon replica as needed for each replica being fitted.
+
+.. note::
+
+  Generating photon PDF sets on the fly can add significant overhead to the
+  fitting procedure. It is therefore strongly advised to precompute the photon
+  PDF set using ``vp-setupfit`` before running the fit with ``n3fit``.
+
+In either case, the generated photon PDF set will be stored locally in the
+``photons_qed_path`` specified in the :ref:`profile file <nnprofile>`, so that
+it can be reused in future fits. The folder containing the photon replicas will
+be named as ``photon_theoryID_<theoryid>_fit_<luxset>`` where ``<theoryid>`` and
+``<luxset>`` are the values specified in the runcard. The folder contains a ``npz``
+file for each photon replica, named as ``replica_<replica_id>.npz``. Each replica
+file contains the photon PDF grid computed with FiatLux at :math:`Q_{\rm init} = 100` GeV,
+prior to the evolution through EKO.
+
+.. important::
+
+  Automatic upload to the NNPDF server through ``vp-upload`` is **not**
+  supported at the moment. The user should manually create a ``tar`` archive
+  file containing the photon replicas and upload it to server. Refer to the
+  :ref:`profile file <nnprofile>` to find the remote URL where photon PDF sets are stored.
+
+
+Using ``vp-setupfit`` (preferred)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  In order to trigger the computation of the photon PDF set during ``vp-setupfit``,
+  the user needs to add the flag ``compute_in_setupfit: true`` in the ``fiatlux`` block
+  discussed above. The following is an example of running ``vp-setupfit`` with
+  this flag enabled:
+
+  .. code-block:: bash
+
+      $ vp-setupfit qed_example_runcard.yml
+      [INFO]: Could not find a resource (photonQED): Could not find Photon QED set /user/envs/nnpdf/share/NNPDF/photons_qed/photon_theoryID_703_fit_251127-jcm-lh-qed-001 in theory: 703. Attempting to download it.
+      [ERROR]: Resource not in the remote repository: Photon QED set for TheoryID 703 and luxset 251127-jcm-lh-qed-001 is not available in the remote server.
+      [WARNING]: Photon QED set for theory 703 with luxset 251127-jcm-lh-qed-001 not found. It will be computed in vp-setupfit.
+      [INFO]: Forcing photon computation with FiatLux during setupfit.
+      [WARNING]: Using q2min from runcard
+      [WARNING]: Using w2min from runcard
+      Using Keras backend
+      [INFO]: All requirements processed and checked successfully. Executing actions.
+      ...
+
+  In addition, the user can also specify the optional parameter ``eps_base`` to
+  the ``fiatlux`` block to set the base precision of the FiatLux computation,
+  which controls the precision on final integration of double integral. By
+  default, it is set to ``1e-5``. If the ``compute_in_setupfit`` flag is set to
+  ``true`` despite the photon PDF being available locally, the code will
+  recompute and overwrite the existing photon PDF set
+
+  .. tip::
+
+    During ``vp-setupfit``, the code tries to distribute the computation of the
+    photon replicas among the available CPU cores as specified in the
+    ``maxcores`` key of the runcard. This can significantly speed up the
+    computation of the photon PDF set. Make sure to set ``maxcores`` to a value
+    compatible with the available hardware. For example, using ``maxcores: 64``
+    will try to compute up to 64 photon replicas in parallel using 64 GB of RAM.
+
+  Once the preparation step is completed, and the photon PDF set is generated and stored
+  locally, the user can proceed to run the fit with ``n3fit`` as usual.
+
+Using ``n3fit`` (discouraged)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+  If the user prefers to compute the photon PDF set during the fitting step with
+  ``n3fit``, no additional flag is needed in the runcard and ``vp-setupfit`` is
+  not required beforehand (unless other resources are needed such as the
+  :ref:`theory covariance matrix <vptheorycov-index>`). The code will check for
+  the presence of the photon PDF set locally before starting the fit for each
+  replica. If it is not found, it will proceed to compute the photon replica as
+  needed for each replica being fitted. The following is an example of running
+  ``n3fit`` where the photon PDF set is computed during the fitting step:
+
+  .. code-block:: bash
+
+      $ n3fit qed_example_runcard.yml 1
+      [INFO]: Creating replica output folder in /user/runcards/qed_example_runcard/nnfit/replica_1
+      [WARNING]: Using q2min from runcard
+      [WARNING]: Using w2min from runcard
+      Using Keras backend
+      [WARNING]: No Photon QED set found for Theory 703 with luxset 251127-jcm-lh-qed-001. It will be computed using FiatLux. This may impact performance. It is recommended to precompute the photon set before running the fit. Refer to https://docs.nnpdf.science/tutorials/run-qed-fit.html for more details on precomputing photon PDF sets.
+      [INFO]: All requirements processed and checked successfully. Executing actions.
+
+  .. warning::
+
+    Computing the photon PDF set during `n3fit` with multiple replicas or using GPUs
+    has not been tested and may lead to unexpected behaviour. It is strongly advised to
+    precompute the photon PDF set using ``vp-setupfit`` before running the fit with ``n3fit``.
diff --git a/doc/sphinx/source/vp/download.rst b/doc/sphinx/source/vp/download.rst
index 133b019b25..493c3d642c 100644
--- a/doc/sphinx/source/vp/download.rst
+++ b/doc/sphinx/source/vp/download.rst
@@ -6,8 +6,8 @@ Downloading resources
 ``validphys`` is designed so that, by default, resources stored in known remote
 locations are downloaded automatically and seamlessly used where necessary.
 Available resources include PDF sets, completed fits, theories, and results of
-past ``validphys`` runs that have been :ref:`uploaded to the server <upload>`. 
-The ``vp-get`` tool, :ref:`described below <vp-get>`, 
+past ``validphys`` runs that have been :ref:`uploaded to the server <upload>`.
+The ``vp-get`` tool, :ref:`described below <vp-get>`,
 can be used to download the same items manually.
 
 Automatic operation
@@ -15,7 +15,7 @@ Automatic operation
 
 By default when some resource such as a PDF is required by ``validphys`` (or
 derived tools such as ``vp-setupfit``), the code will first look for it in some
-local directory specified in the [profile file](nnprofile). If it is not found
+local directory specified in the :ref:`profile file <nnprofile>`. If it is not found
 there, it will try to download it from some remote repository (also specified in
 the profile).
 
@@ -70,6 +70,21 @@ Theories
     Theories (specified by the ``theoryid`` key) are downloaded and
     uncompressed.
 
+Photon PDF sets
+    Photon PDF sets can be downloaded if they have been previously generated
+    using `FiatLux <https://github.com/scarrazza/fiatlux/>`_. A photon PDF set
+    is uniquely identified by the input PDF set used to generate the photon replicas
+    and the theory ID. See :ref:`this tutorial <run-qed-fit>` for more details on
+    how to set up a QED fit.
+
+EKOs
+    Evolution kernels (EKOs) can be downloaded if they have been previously
+    generated using EKO.
+
+Hyperparameter scans
+    Results of the hyperparameter scans generated in the hyperoptimization
+    procedure.
+
 ``validphys`` output files
     Files produced by ``validphys`` can be used as input to subsequent validphys
     analyses (for example χ² tables are used for αs fits). The user needs to
@@ -80,7 +95,7 @@ Theories
 
 .. _vp-get:
 
-The `vp-get` tool
+The ``vp-get`` tool
 -----------------
 
 The ``vp-get`` tool can be used to download resources manually, in the same way
@@ -99,8 +114,11 @@ They correspond to the resources described :ref:`above <what-can-be-downloaded>`
 
   $ vp-get --list
   Available resource types:
+  - eko
   - fit
+  - hyperscan
   - pdf
+  - photonQED
   - theoryID
   - vp_output_file
 
@@ -118,6 +136,15 @@ information on it and bail out:
   $ vp-get fit NNPDF31_nlo_as_0118_1000
   FitSpec(name='NNPDF31_nlo_as_0118_1000', path=PosixPath('/home/zah/anaconda3/envs/nnpdf-dev/share/NNPDF/results/NNPDF31_nlo_as_0118_1000'))
 
+.. note::
+
+  In order to download photon PDF sets, the user needs to specify both the input PDF
+  set and the theory ID. For example:
+
+  .. code:: bash
+
+    $ vp-get photonQED 40000000 NNPDF40_nnlo_as_01180
+
 Downloading resources in code (``validphys.loader``)
 ----------------------------------------------------
 
@@ -164,10 +191,10 @@ Conversely the ``Loader`` class will only search locally.
   ----> 1 l.check_dataset('NMC', theoryid=151)
 
   ~/nngit/nnpdf/validphys2/src/validphys/loader.py in check_dataset(self, name, rules, sysnum, theoryid, cfac, frac, cuts, use_fitcommondata, fit, weight)
-      416 
+      416
       417         if not isinstance(theoryid, TheoryIDSpec):
   --> 418             theoryid = self.check_theoryID(theoryid)
-      419 
+      419
       420         theoryno, _ = theoryid
 
   ~/nngit/nnpdf/validphys2/src/validphys/loader.py in check_theoryID(self, theoryID)
@@ -175,7 +202,7 @@ Conversely the ``Loader`` class will only search locally.
       289             raise TheoryNotFound(("Could not find theory %s. "
   --> 290                   "Folder '%s' not found") % (theoryID, theopath) )
       291         return TheoryIDSpec(theoryID, theopath)
-      292 
+      292
 
   TheoryNotFound: Could not find theory 151. Folder '/home/zah/anaconda3/share/NNPDF/data/theory_151' not found
 
@@ -184,10 +211,9 @@ Output files uploaded to the ``validphys`` can be retrieved specifying their pat
 (starting from the report ID). They will be either downloaded (when using
 ``FallbackLoader``) or retrieved from the cache:
 
-.. code:: python 
+.. code:: python
 
   from validphys.loader import FallbackLoader as Loader
   l = Loader()
   l.check_vp_output_file('qTpvLZLwS924oAsmpMzhFw==/figures/f_ns0_fitunderlyinglaw_plot_closure_pdf_histograms_0.pdf')
   PosixPath('/home/zah/anaconda3/share/NNPDF/vp-cache/qTpvLZLwS924oAsmpMzhFw==/figures/f_ns0_fitunderlyinglaw_plot_closure_pdf_histograms_0.pdf')
-
diff --git a/doc/sphinx/source/vp/nnprofile.rst b/doc/sphinx/source/vp/nnprofile.rst
index b917164914..e4242867be 100644
--- a/doc/sphinx/source/vp/nnprofile.rst
+++ b/doc/sphinx/source/vp/nnprofile.rst
@@ -39,16 +39,16 @@ the code. These should be specified in YAML format.
     It is possible to set the special key ``RELATIVE_TO_PYTHON``, in this case the code
     will use as share folder the share folder of the current environment (for instance ``${CONDA_PREFIX}/share/NNPDF``).
 
-``theories_path``
-    The path in the user's system where the theory files (FKtables and ekos)
-    are to be found, and stored when :ref:`downloaded <download>`.
-    Defaults to ``nnpdf_share/theories``.
-
 ``results_path``
     A path where completed fits are to be retrieved from,
     and stored when :ref:`downloaded <download>`.
     Defaults to ``nnpdf_share/results``.
 
+``theories_path``
+    The path in the user's system where the theory files (FKtables and ekos)
+    are to be found, and stored when :ref:`downloaded <download>`.
+    Defaults to ``nnpdf_share/theories``.
+
 ``data_path``
     List of paths where to read the data from.
     Regardless of the content of this variable, the installation path of the ``nnpdf_data`` package
@@ -58,6 +58,10 @@ the code. These should be specified in YAML format.
 ``validphys_cache_path``
     A path where to store downloaded validphys resources.
 
+``photons_qed_path``
+    A path where to store downloaded photon PDF sets generated with FiatLux. See
+    :ref:`this tutorial <run-qed-fit>` for more details.
+
 ``fit_urls``
     A list of URLs where to search completed fits from.
 
@@ -81,6 +85,12 @@ the code. These should be specified in YAML format.
 ``nnpdf_pdfs_index``
     The name of the remote PDF index. Shouldn't be changed.
 
+``photon_qed_urls``
+    A list of URLs pointing to repositories where photon PDF sets are stored.
+
+``photon_qed_index``
+    The name of the remote photon PDF index. Shouldn't be changed.
+
 ``upload_host``
     The SSH host (with user name as in ``user@host``) used to upload ``validphys`` reports and fits.
 
diff --git a/n3fit/src/n3fit/checks.py b/n3fit/src/n3fit/checks.py
index 8424181c5f..e8f17e2c70 100644
--- a/n3fit/src/n3fit/checks.py
+++ b/n3fit/src/n3fit/checks.py
@@ -452,20 +452,24 @@ def check_deprecated_options(fitting):
     for option in nnfit_options:
         if option in fitting:
             log.warning("'fitting::%s' is an nnfit-only key, it will be ignored", option)
-        
+
 
 @make_argcheck
 def check_photonQED_exists(theoryid, fiatlux):
     """Check that the Photon QED set for this theoryid and luxset exists"""
     if fiatlux is not None:
-      luxset = fiatlux['luxset']
-      try:
-          _ = Loader().check_photonQED(theoryid.id, luxset)
-          log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
-      except FileNotFoundError:
-          log.warning(f"No Photon QED set found for {theoryid} with luxset {luxset}. It "\
-                      "will be compute using FiatLux. This may impact performance. Make "
-                      "sure vp-setupfit has been run prior to n3fit to download necessary resources.")
+        luxset = fiatlux['luxset']
+        try:
+            _ = Loader().check_photonQED(theoryid.id, luxset)
+            log.info(f"Photon QED set found for {theoryid.id} with luxset {luxset}.")
+        except FileNotFoundError:
+            log.warning(
+                f"No Photon QED set found for {theoryid} with luxset {luxset}. It "
+                "will be computed using FiatLux. This may impact performance. It "
+                "is recommended to precompute the photon set before running the fit. "
+                "Refer to https://docs.nnpdf.science/tutorials/run-qed-fit.html for more details "
+                "on precomputing photon PDF sets."
+            )
 
 
 @make_argcheck
diff --git a/n3fit/src/n3fit/scripts/vp_setupfit.py b/n3fit/src/n3fit/scripts/vp_setupfit.py
index 6b00419e3d..849207ecc8 100644
--- a/n3fit/src/n3fit/scripts/vp_setupfit.py
+++ b/n3fit/src/n3fit/scripts/vp_setupfit.py
@@ -1,8 +1,8 @@
 #!/usr/bin/env python
 """
-    setup-fit - prepare and apply data cuts before fit
-    setup-fit constructs the fit [results] folder where data used by nnfit
-    will be stored.
+setup-fit - prepare and apply data cuts before fit
+setup-fit constructs the fit [results] folder where data used by nnfit
+will be stored.
 """
 
 # Implementation notes
@@ -37,7 +37,7 @@
 from reportengine import colors
 from validphys.app import App
 from validphys.config import Config, ConfigError, Environment, EnvironmentError_
-from validphys.loader import FallbackLoader, TheoryNotFound, PhotonQEDNotFound
+from validphys.loader import FallbackLoader, PhotonQEDNotFound, TheoryNotFound
 from validphys.utils import yaml_safe
 
 loader = FallbackLoader()
@@ -57,7 +57,7 @@
     'validphys.filters',
     'validphys.results',
     'validphys.theorycovariance.construction',
-    'validphys.photon.compute'
+    'validphys.photon.compute',
 ]
 
 SETUPFIT_DEFAULTS = dict(use_cuts='internal')
@@ -193,22 +193,25 @@ def from_yaml(cls, o, *args, **kwargs):
                 log.info(f"Photon QED set found for {theoryid} with luxset {luxset}.")
             except PhotonQEDNotFound:
                 if compute_in_setupfit:
-                  log.warning(f"Photon QED set for {theoryid} with luxset {luxset} not found. " \
-                              "It will be computed in vp-setupfit.")
+                    log.warning(
+                        f"Photon QED set for theory {theoryid} with luxset {luxset} not found. "
+                        "It will be computed in vp-setupfit."
+                    )
                 else:
-                  log.warning(f"No photon set found for {theoryid} with luxset {luxset}. Consider "\
-                              "using `compute_in_setupfit` in the runcard to compute all replicas of the photon "
-                              "in vp-setupfit. Otherwise n3fit " \
-                              "will take care of the photon computation. May impact performance.")
+                    log.warning(
+                        f"No photon set found for theory {theoryid} with luxset {luxset}. Consider "
+                        "using `compute_in_setupfit` in the runcard to compute all replicas of the photon "
+                        "in vp-setupfit. Otherwise n3fit "
+                        "will take care of the photon computation. May impact performance."
+                    )
 
             if compute_in_setupfit:
-              log.info("Forcing photon computation with FiatLux.")
-              # Since the photon will be computed, check that the luxset and additional_errors exist
-              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset_exists')
-              if fiatlux.get("additional_errors"):
-                  SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
-              SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux::theory compute_photon_to_disk')
-
+                log.info("Forcing photon computation with FiatLux during setupfit.")
+                # Since the photon will be computed, check that the luxset and additional_errors exist
+                SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_luxset_exists')
+                if fiatlux.get("additional_errors"):
+                    SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux check_additional_errors')
+                SETUPFIT_FIXED_CONFIG['actions_'].append('fiatlux::theory compute_photon_to_disk')
 
         # Check positivity bound
         if file_content.get('positivity_bound') is not None:
diff --git a/validphys2/src/validphys/loader.py b/validphys2/src/validphys/loader.py
index 1aadfcbc6b..b27e4222dc 100644
--- a/validphys2/src/validphys/loader.py
+++ b/validphys2/src/validphys/loader.py
@@ -216,12 +216,13 @@ def available_ekos(self):
         return {
             eko_path.parent.name.split("_")[1] for eko_path in self._theories_path.glob("*/eko.tar")
         }
-    
+
     @functools.cached_property
     def available_photons(self):
         """Return a string token for each of the available theories"""
         return {
-            photon_path.name.split("photon_")[1] for photon_path in self._photons_qed_path.glob("photon_*")
+            photon_path.name.split("photon_")[1]
+            for photon_path in self._photons_qed_path.glob("photon_*")
         }
 
     @property
@@ -319,13 +320,15 @@ def check_eko(self, theoryID):
         if not eko_path.exists():
             raise EkoNotFound(f"Could not find eko {eko_path} in theory: {theoryID}")
         return eko_path
-    
+
     @functools.lru_cache
     def check_photonQED(self, theoryID, luxset):
         """Check the Photon QED set exists and return the path to it"""
         photon_qed_path = self._photons_qed_path / f"photon_theoryID_{int(theoryID)}_fit_{luxset}"
         if not photon_qed_path.exists():
-            raise PhotonQEDNotFound(f"Could not find Photon QED set {photon_qed_path} in theory: {int(theoryID)}")
+            raise PhotonQEDNotFound(
+                f"Could not find Photon QED set {photon_qed_path} in theory: {int(theoryID)}"
+            )
         return photon_qed_path
 
     @property
@@ -837,12 +840,12 @@ def eko_index(self):
     @_key_or_loader_error
     def eko_urls(self):
         return self.nnprofile['eko_urls']
-    
+
     @property
     @_key_or_loader_error
     def photon_qed_index(self):
         return self.nnprofile['photon_qed_index']
-    
+
     @property
     @_key_or_loader_error
     def photon_qed_urls(self):
@@ -917,7 +920,7 @@ def remote_ekos(self):
         token = 'eko_'
         rt = self.remote_files(self.eko_urls, self.eko_index, thing="ekos")
         return {k[len(token) :]: v for k, v in rt.items()}
-    
+
     @property
     @functools.lru_cache
     def remote_photons(self):
@@ -958,7 +961,7 @@ def downloadable_theories(self):
     @property
     def downloadable_ekos(self):
         return list(self.remote_ekos)
-    
+
     @property
     def downloadable_photons(self):
         return list(self.remote_photons)
@@ -1150,12 +1153,13 @@ def download_photonQED(self, thid, luxset: str):
         remote = self.remote_photons
         key = f"theoryID_{thid}_fit_{luxset}"
         if key not in remote:
-            raise PhotonQEDNotFound(f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server")
+            raise PhotonQEDNotFound(
+                f"Photon QED set for TheoryID {thid} and luxset {luxset} is not available in the remote server."
+            )
         # Check that we have the theory we need
         target_path = self._photons_qed_path
         download_and_extract(remote[key], target_path)
-        
-    
+
     def download_vp_output_file(self, filename, **kwargs):
         try:
             root_url = self.nnprofile['reports_root_url']