Merge pull request #1778 from NNPDF/produce_eko_photon_new

Add produce_eko_photon to evolven3fit_new

n3fit/src/evolven3fit_new/eko_utils.py

@@ -48,45 +48,21 @@ def construct_eko_cards(
     op_card_dict and theory_card_dict are optional updates that can be provided respectively to the
     operator card and to the theory card.
     """
-    if theory_card_dict is None:
-        theory_card_dict = {}
-    if op_card_dict is None:
-        op_card_dict = {}
-
-    # theory_card construction
-    theory = Loader().check_theoryID(theoryID).get_description()
-    theory.pop("FNS")
-    theory.update(theory_card_dict)
-
-    # Prepare the thresholds according to MaxNfPdf
-    thresholds = {"c": theory["kcThr"], "b": theory["kbThr"], "t": theory["ktThr"]}
-    if theory["MaxNfPdf"] < 5:
-        thresholds["b"] = np.inf
-    if theory["MaxNfPdf"] < 6:
-        thresholds["t"] = np.inf
+    theory, thresholds = load_theory(theoryID, theory_card_dict)
 
-    if "nfref" not in theory:
-        theory["nfref"] = NFREF_DEFAULT
+    # if is eko_photon then mu0 = q_gamma
+    mu0 = theory["Q0"]
 
     # 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
-        elif mu0 < theory["mb"] * thresholds["b"]:
-            theory["nf0"] = 4
-        elif mu0 < theory["mt"] * thresholds["t"]:
-            theory["nf0"] = 5
-        else:
-            theory["nf0"] = 6
-
-    # Setting the thresholds in the theory card to inf if necessary
-    theory.update({"kbThr": thresholds["b"], "ktThr": thresholds["t"]})
+        theory["nf0"] = find_nf(mu0, theory, thresholds)
 
     # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
+    # construct mugrid
+
     # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default
     q2_grid = utils.generate_q2grid(
         mu0,
@@ -100,8 +76,6 @@ def construct_eko_cards(
         theory["nf0"],
     )
 
-    # construct operator card
-    op_card = default_op_card
     masses = np.array([theory["mc"], theory["mb"], theory["mt"]]) ** 2
     thresholds_ratios = np.array([thresholds["c"], thresholds["b"], thresholds["t"]]) ** 2
 
@@ -113,7 +87,6 @@ def construct_eko_cards(
     # 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:
@@ -126,7 +99,87 @@ def construct_eko_cards(
         else:
             mugrid.append((q, int(nf_default(q2, atlas))))
 
-    op_card.update({"mu0": theory["Q0"], "mugrid": mugrid})
+    # construct operator card
+    op_card = build_opcard(op_card_dict, theory, x_grid, mu0, mugrid)
+
+    return theory_card, op_card
+
+
+def construct_eko_photon_cards(
+    theoryID,
+    q_fin,
+    x_grid,
+    q_gamma,
+    op_card_dict: Optional[Dict[str, Any]] = None,
+    theory_card_dict: Optional[Dict[str, Any]] = None,
+):
+    """
+    Return the theory and operator cards used to construct the eko_photon.
+    theoryID is the ID of the theory for which we are computing the theory and operator card.
+    q_fin is the final point of the q grid while q_points is the number of points of the grid.
+    x_grid is the x grid to be used.
+    op_card_dict and theory_card_dict are optional updates that can be provided respectively to the
+    operator card and to the theory card.
+    """
+    theory, thresholds = load_theory(theoryID, theory_card_dict)
+
+    # if is eko_photon then mu0 = q_gamma
+    mu0 = q_gamma
+
+    # Set nf_0 according to mu0 unless set explicitly
+    if "nf0" not in theory:
+        theory["nf0"] = find_nf(mu0, theory, thresholds)
+
+    # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
+    legacy_class = runcards.Legacy(theory, {})
+    theory_card = legacy_class.new_theory
+
+    q_fin = theory["Q0"]
+
+    nf_fin = find_nf(q_fin, theory, thresholds)
+
+    # construct mugrid
+    mugrid = [(q_fin, nf_fin)]
+
+    # construct operator card
+    op_card = build_opcard(op_card_dict, theory, x_grid, mu0, mugrid)
+
+    return theory_card, op_card
+
+
+def load_theory(theoryID, theory_card_dict):
+    """loads and returns the theory dictionary and the thresholds"""
+    if theory_card_dict is None:
+        theory_card_dict = {}
+    # theory_card construction
+    theory = Loader().check_theoryID(theoryID).get_description()
+    theory.pop("FNS")
+    theory.update(theory_card_dict)
+
+    if "nfref" not in theory:
+        theory["nfref"] = NFREF_DEFAULT
+
+    # Prepare the thresholds according to MaxNfPdf
+    thresholds = {"c": theory["kcThr"], "b": theory["kbThr"], "t": theory["ktThr"]}
+    if theory["MaxNfPdf"] < 5:
+        thresholds["b"] = np.inf
+    if theory["MaxNfPdf"] < 6:
+        thresholds["t"] = np.inf
+
+    # Setting the thresholds in the theory card to inf if necessary
+    theory.update({"kbThr": thresholds["b"], "ktThr": thresholds["t"]})
+
+    return theory, thresholds
+
+
+def build_opcard(op_card_dict, theory, x_grid, mu0, mugrid):
+    """builds the opcard"""
+    if op_card_dict is None:
+        op_card_dict = {}
+
+    op_card = default_op_card
+
+    op_card.update({"mu0": mu0, "mugrid": mugrid})
 
     op_card["xgrid"] = x_grid
     # Specific defaults for evolven3fit evolution
@@ -144,5 +197,30 @@ def construct_eko_cards(
         elif key in op_card_dict:
             _logger.warning("Entry %s is not a dictionary and will be ignored", key)
 
+    # if no -e was given, take ev_op_iterations from EVOLVEN3FIT_CONFIGS_DEFAULTS_{TRN,EXA}
+    if op_card['configs']['ev_op_iterations'] is None:
+        if theory["ModEv"] == "TRN":
+            op_card['configs']['ev_op_iterations'] = EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN[
+                "ev_op_iterations"
+            ]
+        if theory["ModEv"] == "EXA":
+            op_card['configs']['ev_op_iterations'] = EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA[
+                "ev_op_iterations"
+            ]
+
     op_card = runcards.OperatorCard.from_dict(op_card)
-    return theory_card, op_card
+
+    return op_card
+
+
+def find_nf(mu, theory, thresholds):
+    """compute nf for a given mu"""
+    if mu < theory["mc"] * thresholds["c"]:
+        nf = 3
+    elif mu < theory["mb"] * thresholds["b"]:
+        nf = 4
+    elif mu < theory["mt"] * thresholds["t"]:
+        nf = 5
+    else:
+        nf = 6
+    return nf

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=None,
+        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_photon_cards(
+                args.theoryID, args.q_fin, x_grid, args.q_gamma, op_card_info, theory_card_info,
+            )
         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))