Add photon

conda-recipe/meta.yaml

@@ -51,9 +51,9 @@ requirements:
         - sphinxcontrib-bibtex
         - curio >=1.0
         - pineappl >=0.5.8
-        - eko >=0.12.0
+        - eko >=0.12.2
         - banana-hep >=0.6.8
-
+        - fiatlux
 
 test:
     requires:

doc/sphinx/source/n3fit/runcard_detailed.rst

@@ -419,4 +419,5 @@ It is however possible to disable them by setting to false the ``sum_rules`` fla
 
 
 It is also possible to impose just the valence or the momentum sum rules by using the
-``VSR`` or ``MSR`` flags, respectively (``True`` is equal to ``All``).
+``VSR`` or ``MSR`` flags, respectively (``True`` is equal to ``All``).
+

doc/sphinx/source/references.bib

@@ -604,3 +604,32 @@ @article{Carrazza:2016htc
     pages = "205",
     year = "2016"
 }
+
+@article{Manohar:2017eqh,
+    author = "Manohar, Aneesh V. and Nason, Paolo and Salam, Gavin P. and Zanderighi, Giulia",
+    title = "{The Photon Content of the Proton}",
+    eprint = "1708.01256",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "CERN-TH-2017-141",
+    doi = "10.1007/JHEP12(2017)046",
+    journal = "JHEP",
+    volume = "12",
+    pages = "046",
+    year = "2017"
+}
+
+@article{Manohar:2016nzj,
+    author = "Manohar, Aneesh and Nason, Paolo and Salam, Gavin P. and Zanderighi, Giulia",
+    title = "{How bright is the proton? A precise determination of the photon parton distribution function}",
+    eprint = "1607.04266",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "CERN-TH-2016-155",
+    doi = "10.1103/PhysRevLett.117.242002",
+    journal = "Phys. Rev. Lett.",
+    volume = "117",
+    number = "24",
+    pages = "242002",
+    year = "2016"
+}

doc/sphinx/source/tutorials/run-fit.md

@@ -256,3 +256,8 @@ Iterate the fit
 
 It may be desirable to iterate a fit to achieve a higher degree of convergence/stability in the fit.
 To read more about this, see [How to run an iterated fit](run-iterated-fit).
+
+QED fit
+-------
+
+In order to run a QED fit see [How to run a QED fit](run-qed-fit)

doc/sphinx/source/tutorials/run-qed-fit.rst

@@ -0,0 +1,30 @@
+.. _run-qed-fit:
+
+==========================
+How to run a QED fit
+==========================
+
+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:
+
+.. code-block:: yaml
+
+    fiatlux:
+      luxset: NNPDF40_nnlo_as_01180
+      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 the 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.
+
+Whenever the photon PDF is generated, it will remain constant during the fit and will enter in the `MSR`.      

n3fit/runcards/examples/Basic_runcard_qed.yml

@@ -0,0 +1,171 @@
+#
+# Configuration file for n3fit
+#
+############################################################
+description: Basic runcard qed
+
+############################################################
+# frac: training fraction
+# ewk: apply ewk k-factors
+# sys: systematics treatment (see systypes)
+dataset_inputs:
+# - {dataset: NMCPD_dw_ite, frac: 0.75}
+- {dataset: NMC, frac: 0.75}
+- {dataset: SLACP_dwsh, frac: 0.75}
+# - {dataset: SLACD_dw_ite, frac: 0.75}
+# - {dataset: BCDMSP_dwsh, frac: 0.75}
+# - {dataset: BCDMSD_dw_ite, frac: 0.75}
+# - {dataset: CHORUSNUPb_dw_ite, frac: 0.75}
+# - {dataset: CHORUSNBPb_dw_ite, frac: 0.75}
+# - {dataset: NTVNUDMNFe_dw_ite, frac: 0.75, cfac: [MAS]}
+# - {dataset: NTVNBDMNFe_dw_ite, frac: 0.75, cfac: [MAS]}
+# - {dataset: HERACOMBNCEM, frac: 0.75}
+# - {dataset: HERACOMBNCEP460, frac: 0.75}
+# - {dataset: HERACOMBNCEP575, frac: 0.75}
+# - {dataset: HERACOMBNCEP820, frac: 0.75}
+# - {dataset: HERACOMBNCEP920, frac: 0.75}
+# - {dataset: HERACOMBCCEM, frac: 0.75}
+# - {dataset: HERACOMBCCEP, frac: 0.75}
+# - {dataset: HERACOMB_SIGMARED_C, frac: 0.75}
+# - {dataset: HERACOMB_SIGMARED_B, frac: 0.75}
+# - {dataset: DYE886R_dw_ite, frac: 0.75, cfac: [QCD]}
+# - {dataset: DYE886P, frac: 0.75, cfac: [QCD]}
+# - {dataset: DYE605_dw_ite, frac: 0.75, cfac: [QCD]}
+# - {dataset: DYE906R_dw_ite, frac: 0.75, cfac: [ACC, QCD]}
+# - {dataset: CDFZRAP_NEW, frac: 0.75, cfac: [QCD]}
+# - {dataset: D0ZRAP_40, frac: 0.75, cfac: [QCD]}
+# - {dataset: D0WMASY, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASWZRAP36PB, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASZHIGHMASS49FB, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASLOMASSDY11EXT, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASWZRAP11CC, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASWZRAP11CF, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASDY2D8TEV, frac: 0.75, cfac: [QCDEWK]}
+# - {dataset: ATLAS_DY_2D_8TEV_LOWMASS, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_WZ_TOT_13TEV, frac: 0.75, cfac: [NRM, QCD]}
+# - {dataset: ATLAS_WP_JET_8TEV_PT, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_WM_JET_8TEV_PT, frac: 0.75, cfac: [QCD]}
+- {dataset: ATLASZPT8TEVMDIST, frac: 0.75, cfac: [QCD], sys: 10}
+# - {dataset: ATLASZPT8TEVYDIST, frac: 0.75, cfac: [QCD], sys: 10}
+# - {dataset: ATLASTTBARTOT7TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASTTBARTOT8TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_TTBARTOT_13TEV_FULLLUMI, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_TTB_DIFF_8TEV_LJ_TRAPNORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_TTB_DIFF_8TEV_LJ_TTRAPNORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_TOPDIFF_DILEPT_8TEV_TTRAPNORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_1JET_8TEV_R06_DEC, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_2JET_7TEV_R06, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLASPHT15_SF, frac: 0.75, cfac: [QCD, EWK]}
+# - {dataset: ATLAS_SINGLETOP_TCH_R_7TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_SINGLETOP_TCH_R_13TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_7TEV_T_RAP_NORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_7TEV_TBAR_RAP_NORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_8TEV_T_RAP_NORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: ATLAS_SINGLETOP_TCH_DIFF_8TEV_TBAR_RAP_NORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSWEASY840PB, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSWMASY47FB, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSDY2D11, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSWMU8TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSZDIFF12, frac: 0.75, cfac: [QCD, NRM], sys: 10}
+# - {dataset: CMS_2JET_7TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_1JET_8TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSTTBARTOT7TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSTTBARTOT8TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSTTBARTOT13TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSTOPDIFF8TEVTTRAPNORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMSTTBARTOT5TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_TTBAR_2D_DIFF_MTT_TRAP_NORM, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_TTB_DIFF_13TEV_2016_2L_TRAP, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_TTB_DIFF_13TEV_2016_LJ_TRAP, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_SINGLETOP_TCH_TOT_7TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_SINGLETOP_TCH_R_8TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: CMS_SINGLETOP_TCH_R_13TEV, frac: 0.75, cfac: [QCD]}
+# - {dataset: LHCBZ940PB, frac: 0.75, cfac: [QCD]}
+# - {dataset: LHCBZEE2FB_40, frac: 0.75, cfac: [QCD]}
+# - {dataset: LHCBWZMU7TEV, frac: 0.75, cfac: [NRM, QCD]}
+# - {dataset: LHCBWZMU8TEV, frac: 0.75, cfac: [NRM, QCD]}
+# - {dataset: LHCB_Z_13TEV_DIMUON, frac: 0.75, cfac: [QCD]}
+# - {dataset: LHCB_Z_13TEV_DIELECTRON, frac: 0.75, cfac: [QCD]}
+
+############################################################
+datacuts:
+  t0pdfset: NNPDF40_nnlo_as_01180  # PDF set to generate t0 covmat
+  q2min: 3.49                        # Q2 minimum
+  w2min: 12.5                        # W2 minimum
+
+############################################################
+theory:
+  theoryid: 522     # database id
+
+############################################################
+trvlseed: 1551864071
+nnseed: 676150632
+mcseed: 619859729
+save: false
+genrep: true      # true = generate MC replicas, false = use real data
+
+parameters: # This defines the parameter dictionary that is passed to the Model Trainer
+  nodes_per_layer: [25, 20, 8]
+  activation_per_layer: [tanh, tanh, linear]
+  initializer: glorot_normal
+  optimizer:
+    clipnorm: 6.073e-6
+    learning_rate: 2.621e-3
+    optimizer_name: Nadam
+  epochs: 900
+  positivity:
+    initial: 184.8
+    multiplier:
+  integrability:
+    initial: 10
+    multiplier:
+  stopping_patience: 0.1
+  layer_type: dense
+  dropout: 0.0
+  threshold_chi2: 3.5
+
+fitting:
+  fitbasis: EVOL  # EVOL (7), EVOLQED (8), etc.
+  basis:
+  - {fl: sng, trainable: false, smallx: [1.089, 1.117], largex: [1.462, 3.008]}
+  - {fl: g, trainable: false, smallx: [0.7542, 1.105], largex: [2.826, 5.407]}
+  - {fl: v, trainable: false, smallx: [0.4715, 0.7253], largex: [1.564, 3.48]}
+  - {fl: v3, trainable: false, smallx: [0.1372, 0.4205], largex: [1.755, 3.451]}
+  - {fl: v8, trainable: false, smallx: [0.5641, 0.7702], largex: [1.513, 3.433]}
+  - {fl: t3, trainable: false, smallx: [-0.4942, 0.9992], largex: [1.751, 3.383]}
+  - {fl: t8, trainable: false, smallx: [0.532, 0.8572], largex: [1.541, 3.349]}
+  - {fl: t15, trainable: false, smallx: [1.052, 1.14], largex: [1.487, 3.09]}
+
+############################################################
+positivity:
+  posdatasets:
+  - {dataset: POSF2U, maxlambda: 1e6}        # Positivity Lagrange Multiplier
+  # - {dataset: POSF2DW, maxlambda: 1e6}
+  # - {dataset: POSF2S, maxlambda: 1e6}
+  - {dataset: POSFLL, maxlambda: 1e6}
+  # - {dataset: POSDYU, maxlambda: 1e10}
+  # - {dataset: POSDYD, maxlambda: 1e10}
+  # - {dataset: POSDYS, maxlambda: 1e10}
+  # - {dataset: POSF2C, maxlambda: 1e6}
+  # - {dataset: POSXUQ, maxlambda: 1e6}        # Positivity of MSbar PDFs
+  # - {dataset: POSXUB, maxlambda: 1e6}
+  # - {dataset: POSXDQ, maxlambda: 1e6}
+  # - {dataset: POSXDB, maxlambda: 1e6}
+  # - {dataset: POSXSQ, maxlambda: 1e6}
+  # - {dataset: POSXSB, maxlambda: 1e6}
+  # - {dataset: POSXGL, maxlambda: 1e6}
+
+############################################################
+integrability:
+  integdatasets:
+  # - {dataset: INTEGXT8, maxlambda: 1e2}
+  - {dataset: INTEGXT3, maxlambda: 1e2}
+
+############################################################
+debug: True
+maxcores: 8
+
+fiatlux:
+  luxset: NNPDF40_nnlo_as_01180
+  additional_errors: true # should be set to true only for the last iteration
+  luxseed: 1234567890

n3fit/src/evolven3fit_new/eko_utils.py

@@ -10,14 +10,22 @@
 
 _logger = logging.getLogger(__name__)
 
-EVOLVEN3FIT_CONFIGS_DEFAULTS = {
+EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN = {
     "ev_op_iterations": 1,
     "ev_op_max_order": (1, 0),
     "evolution_method": "truncated",
     "inversion_method": "expanded",
     "n_integration_cores": 1,
 }
 
+EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
+    "ev_op_iterations": 10,
+    "ev_op_max_order": (1, 0),
+    "evolution_method": "iterate-exact",
+    "inversion_method": "exact",
+    "n_integration_cores": 1,
+}
+
 NFREF_DEFAULT = 5
 NF0_DEFAULT = 4
 
@@ -68,7 +76,10 @@ def construct_eko_cards(
     )
     op_card["rotations"]["xgrid"] = x_grid
     # Specific defaults for evolven3fit evolution
-    op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS)
+    if theory["ModEv"] == "TRN":
+        op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN)
+    if theory["ModEv"] == "EXA":
+        op_card["configs"].update(EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA)
     # User can still change the configs via op_card_dict
 
     # Note that every entry that is not a dictionary should not be

n3fit/src/evolven3fit_new/evolve.py

@@ -6,8 +6,7 @@
 import eko
 from eko import basis_rotation
 from eko import runner
-from ekobox import genpdf, info_file
-from ekomark import apply
+from ekobox import genpdf, info_file, apply
 from reportengine.compat import yaml
 from validphys.loader import Loader
 
@@ -69,7 +68,7 @@ def evolve_fit(
     stdout_log = logging.StreamHandler(sys.stdout)
     for log in [log_file, stdout_log]:
         log.setLevel(LOGGING_SETTINGS["level"])
-        log.setFormatter(logging.Formatter(LOGGING_SETTINGS["formatter"]))
+        log.setFormatter(LOGGING_SETTINGS["formatter"])
     for logger in (_logger, *[logging.getLogger("eko")]):
         logger.handlers = []
         logger.setLevel(LOGGING_SETTINGS["level"])
@@ -85,6 +84,7 @@ def evolve_fit(
     theory, op = eko_utils.construct_eko_cards(
         theoryID, q_fin, q_points, x_grid, op_card_dict, theory_card_dict
     )
+    qed = theory.order[1] > 0
     if eko_path is not None:
         eko_path = pathlib.Path(eko_path)
         _logger.info(f"Loading eko from : {eko_path}")
@@ -105,10 +105,12 @@ def evolve_fit(
     info["XMin"] = float(x_grid[0])
     info["XMax"] = float(x_grid[-1])
     with eko.EKO.read(eko_path) as eko_op:
-        info["AlphaS_Qs"] = eko_op.mu2grid.tolist()
+        if eko.__version__ >= "0.13":
+            raise ModuleNotFoundError("Please, fix evolven3fit np.sqrt(Q) hack")
+        info["AlphaS_Qs"] = np.sqrt(info["AlphaS_Qs"]).tolist()
         dump_info_file(usr_path, info)
         for replica in initial_PDFs_dict.keys():
-            evolved_block = evolve_exportgrid(initial_PDFs_dict[replica], eko_op, x_grid)
+            evolved_block = evolve_exportgrid(initial_PDFs_dict[replica], eko_op, x_grid, qed)
             dump_evolved_replica(
                 evolved_block, usr_path, int(replica.removeprefix("replica_"))
             )
@@ -142,7 +144,7 @@ def load_fit(usr_path):
     return pdf_dict
 
 
-def evolve_exportgrid(exportgrid, eko, x_grid):
+def evolve_exportgrid(exportgrid, eko, x_grid, qed):
     """
     Evolves the provided exportgrid for the desired replica with the eko and returns the evolved block
 
@@ -154,6 +156,8 @@ def evolve_exportgrid(exportgrid, eko, x_grid):
             eko operator for evolution
         xgrid: list
             xgrid to be used as the targetgrid
+        qed: bool
+            whether qed is activated or not
     Returns
     -------
         : np.array
@@ -163,7 +167,7 @@ def evolve_exportgrid(exportgrid, eko, x_grid):
     pdf_grid = np.array(exportgrid["pdfgrid"]).transpose()
     pdf_to_evolve = utils.LhapdfLike(pdf_grid, exportgrid["q20"], x_grid)
     # evolve pdf
-    evolved_pdf = apply.apply_pdf(eko, pdf_to_evolve)
+    evolved_pdf = apply.apply_pdf(eko, pdf_to_evolve, qed=qed)
     # generate block to dump
     targetgrid = eko.rotations.targetgrid.tolist()