Produce eko photon

n3fit/src/evolven3fit_new/eko_utils.py

@@ -39,6 +39,8 @@ def construct_eko_cards(
     x_grid,
     op_card_dict: Optional[Dict[str, Any]] = None,
     theory_card_dict: Optional[Dict[str, Any]] = None,
+    q_gamma = None,
+    is_eko_photon = False
 ):
     """
     Return the theory and operator cards used to construct the eko.
@@ -68,8 +70,13 @@ def construct_eko_cards(
     if "nfref" not in theory:
         theory["nfref"] = NFREF_DEFAULT
 
+    # if is eko_photon then mu0 = q_gamma
+    if not is_eko_photon:
+        mu0 = theory["Q0"]
+    else:
+        mu0 = q_gamma
+
     # Set nf_0 according to the fitting scale unless set explicitly
-    mu0 = theory["Q0"]
     if "nf0" not in theory:
         if mu0 < theory["mc"] * thresholds["c"]:
             theory["nf0"] = 3
@@ -87,18 +94,31 @@ def construct_eko_cards(
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
-    # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default
-    q2_grid = utils.generate_q2grid(
-        mu0,
-        q_fin,
-        q_points,
-        {
-            theory["mc"]: thresholds["c"],
-            theory["mb"]: thresholds["b"],
-            theory["mt"]: thresholds["t"],
-        },
-        theory["nf0"],
-    )
+    if not is_eko_photon:
+        # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default
+        q2_grid = utils.generate_q2grid(
+            mu0,
+            q_fin,
+            q_points,
+            {
+                theory["mc"]: thresholds["c"],
+                theory["mb"]: thresholds["b"],
+                theory["mt"]: thresholds["t"],
+            },
+            theory["nf0"],
+        )
+    else:
+        q_fin = theory["Q0"]
+        # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default
+        q2_grid = [q_fin]
+        if q_fin < theory["mc"] * thresholds["c"]:
+            nf_fin = 3
+        elif q_fin < theory["mb"] * thresholds["b"]:
+            nf_fin = 4
+        elif q_fin < theory["mt"] * thresholds["t"]:
+            nf_fin = 5
+        else:
+            nf_fin = 6
 
     # construct operator card
     op_card = default_op_card
@@ -110,27 +130,31 @@ def construct_eko_cards(
         origin=(mu0**2, theory["nf0"]),
     )
 
-    # Create the eko operator q2grid
-    # This is a grid which contains information on (q, nf)
-    # in VFNS values at the matching scales need to be doubled so that they are considered in both sides
-
-    ep = 1e-4
-    mugrid = []
-    for q2 in q2_grid:
-        q = float(np.sqrt(q2))
-        if nf_default(q2 + ep, atlas) != nf_default(q2 - ep, atlas):
-            nf_l = int(nf_default(q2 - ep, atlas))
-            nf_u = int(nf_default(q2 + ep, atlas))
-            mugrid.append((q, nf_l))
-            mugrid.append((q, nf_u))
-        else:
-            mugrid.append((q, int(nf_default(q2, atlas))))
-
-    op_card.update({"mu0": theory["Q0"], "mugrid": mugrid})
+    if not is_eko_photon:
+        # Create the eko operator q2grid
+        # This is a grid which contains information on (q, nf)
+        # in VFNS values at the matching scales need to be doubled so that they are considered in both sides
+        ep = 1e-4
+        mugrid = []
+        for q2 in q2_grid:
+            q = float(np.sqrt(q2))
+            if nf_default(q2 + ep, atlas) != nf_default(q2 - ep, atlas):
+                nf_l = int(nf_default(q2 - ep, atlas))
+                nf_u = int(nf_default(q2 + ep, atlas))
+                mugrid.append((q, nf_l))
+                mugrid.append((q, nf_u))
+            else:
+                mugrid.append((q, int(nf_default(q2, atlas))))
+
+        op_card.update({"mu0": theory["Q0"], "mugrid": mugrid})
+    else:
+        op_card.update({"mu0": mu0, "mugrid": [(q_fin, nf_fin)]})
 
     op_card["xgrid"] = x_grid
     # Specific defaults for evolven3fit evolution
     if theory["ModEv"] == "TRN":
+        if op_card_dict["configs"]["ev_op_iterations"]:
+            _logger.warning("Provided ev_op_iterations for a TRN theory. It will be unused")
         op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN)
     if theory["ModEv"] == "EXA":
         op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA)

n3fit/src/n3fit/scripts/evolven3fit_new.py

@@ -49,6 +49,47 @@ def construct_eko_parser(subparsers):
     )
     return parser
 
+def construct_eko_photon_parser(subparsers):
+    parser = subparsers.add_parser(
+        "produce_eko_photon",
+        help=(
+            """Produce the eko_photon for the specified theory.
+            The q_gamma will be the provided one.
+            The x_grid starts at x_grid_ini and ends at 1.0 and contains the 
+            provided number of points. The eko will be dumped in the provided path."""
+        ),
+    )
+    parser.add_argument(
+        "theoryID", type=int, help="ID of the theory used to produce the eko"
+    )
+    parser.add_argument(
+        "dump",
+        type=pathlib.Path,
+        help="Path where the EKO is dumped",
+    )
+    parser.add_argument(
+        "-i",
+        "--x-grid-ini",
+        default=None,
+        type=float,
+        help="Starting point of the x-grid",
+    )
+    parser.add_argument(
+        "-p",
+        "--x-grid-points",
+        default=None,
+        type=int,
+        help="Number of points of the x-grid",
+    )
+    parser.add_argument(
+        "-g",
+        "--q-gamma",
+        default=100,
+        type=float,
+        help="Scale at which the photon is generated",
+    )
+    return parser
+
 
 def construct_evolven3fit_parser(subparsers):
     parser = subparsers.add_parser(
@@ -107,13 +148,22 @@ def main():
         default=1,
         help="Number of cores to be used",
     )
+    parser.add_argument(
+        "-e",
+        "--ev-op-iterations",
+        type=int,
+        default=30,
+        help="ev_op_iterations for the EXA theory",
+    )
     subparsers = parser.add_subparsers(title="actions", dest="actions")
     eko_parser = construct_eko_parser(subparsers)
+    eko_photon_parser = construct_eko_photon_parser(subparsers)
     evolven3fit_parser = construct_evolven3fit_parser(subparsers)
 
     args = parser.parse_args()
     op_card_info = { "configs" : {
-        "n_integration_cores": args.n_cores,}
+        "n_integration_cores": args.n_cores,
+        "ev_op_iterations": args.ev_op_iterations,}
     }
     theory_card_info = {}
     if args.actions == "evolve":
@@ -127,7 +177,7 @@ def main():
             args.load,
             args.force,
         )
-    elif args.actions == "produce_eko":
+    else:
         stdout_log = logging.StreamHandler(sys.stdout)
         stdout_log.setLevel(evolve.LOGGING_SETTINGS["level"])
         stdout_log.setFormatter(evolve.LOGGING_SETTINGS["formatter"])
@@ -148,9 +198,14 @@ def main():
             )
         else:
             x_grid = np.geomspace(args.x_grid_ini, 1.0, args.x_grid_points)
-        tcard, opcard = eko_utils.construct_eko_cards(
-            args.theoryID, args.q_fin, args.q_points, x_grid, op_card_info, theory_card_info
-        )
+        if args.actions == "produce_eko":
+            tcard, opcard = eko_utils.construct_eko_cards(
+                args.theoryID, args.q_fin, args.q_points, x_grid, op_card_info, theory_card_info
+            )
+        elif args.actions == "produce_eko_photon":
+            tcard, opcard = eko_utils.construct_eko_cards(
+                args.theoryID, args.q_fin, args.q_points, x_grid, op_card_info, theory_card_info, args.q_gamma, is_eko_photon = True
+            )
         runner.solve(tcard, opcard, args.dump)
 
 

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

@@ -77,8 +77,8 @@ def test_params():
     replica = 1
     test_theory = API.theoryid(theoryid=THEORY_QED)
     theory = test_theory.get_description()
-    for channel in ["F2", "FL"]:
-        tmp = "fastkernel/FIATLUX_DIS_" + channel + ".pineappl.lz4"
+    for kind in ["F2", "FL"]:
+        tmp = "fastkernel/FIATLUX_DIS_" + kind + ".pineappl.lz4"
         path_to_fktable = test_theory.path / tmp
         struct_func = sf.InterpStructureFunction(path_to_fktable, pdfs.members[replica])
         np.testing.assert_allclose(struct_func.q2_max, 1e8)
@@ -93,8 +93,8 @@ def test_interpolation_grid():
     pdfs = PDFset(PDF).load()
     test_theory = API.theoryid(theoryid=THEORY_QED)
     for replica in [1, 2, 3]:
-        for channel in ["F2", "FL"]:
-            tmp = "fastkernel/FIATLUX_DIS_" + channel + ".pineappl.lz4"
+        for kind in ["F2", "FL"]:
+            tmp = "fastkernel/FIATLUX_DIS_" + kind + ".pineappl.lz4"
             path_to_fktable = test_theory.path / tmp
             fktable = pineappl.fk_table.FkTable.read(path_to_fktable)
             x = np.unique(fktable.bin_left(1))