Implementation of HT and power corrections

Author: achiefa

doc/sphinx/source/vp/theorycov/power_corrections.rst

@@ -0,0 +1,238 @@
+.. _vptheorycov-pc:
+
+Power corrections
+=================
+
+Power corrections (also referred to as higher twist corrections for DIS-like
+processes) model contributions from non-perturbative effects that scale as
+inverse powers of the hard scale. They are implemented in the
+``theorycovariance`` module and can be included as a theory covariance matrix in
+a fit. Power corrections for jets and higher twists for DIS data have been
+determined in :cite:p:`Ball:2025xtj`, based on NNPDF4.0, where the reader can
+find further details on the methodology and phenomenological implications.
+
+The implementation is in
+`higher_twist_functions.py <https://github.com/NNPDF/nnpdf/tree/master/validphys2/src/validphys/theorycovariance/higher_twist_functions.py>`_
+and
+`construction.py <https://github.com/NNPDF/nnpdf/tree/master/validphys2/src/validphys/theorycovariance/construction.py>`_.
+
+
+Overview
+--------
+
+In NNPDF, power corrections modify theoretical predictions by introducing multiplicative
+shifts. For a generic observable :math:`O`, the corrected prediction is
+
+.. math:: O \to O \times (1 + \mathrm{PC}),
+
+where :math:`\mathrm{PC}` is the power correction. The shift to the prediction
+is therefore
+
+.. math:: \Delta O = O \times \mathrm{PC}.
+
+Different functional forms for the power correction are used depending on the
+process type:
+
+- **DIS** (neutral current and charged current): the correction depends on
+  Bjorken-:math:`x` and :math:`Q^2`, and scales as :math:`1/Q^2`.
+- **Single-inclusive jets**: the correction depends on rapidity and transverse
+  momentum :math:`p_T`, and scales as :math:`1/p_T`.
+- **Dijets**: the correction depends on a rapidity variable and the dijet
+  invariant mass :math:`m_{jj}`, and scales as :math:`1/m_{jj}`.
+
+
+Parametrisation
+---------------
+
+Power corrections are parametrised using a piecewise-linear interpolation
+between a set of nodes. The node positions (``nodes``) and the function values
+at each node (``yshift``) are specified in the runcard.
+
+The interpolation is constructed as a sum of triangular basis functions: each
+node :math:`i` is associated with a triangle that peaks at the node position
+with value ``yshift[i]`` and drops linearly to zero at the two neighbouring
+nodes. The resulting function is continuous and piecewise-linear.
+
+For DIS processes, the nodes are placed in Bjorken-:math:`x` and the power
+correction for a data point at :math:`(x, Q^2)` is
+
+.. math:: \mathrm{PC}(x, Q^2) = \frac{h(x)}{Q^2},
+
+where :math:`h(x)` is the piecewise-linear interpolation.
+
+For jet processes, the nodes are placed in rapidity and the correction at
+:math:`(y, p_T)` is
+
+.. math:: \mathrm{PC}(y, p_T) = \frac{h(y)}{p_T}.
+
+For dijets, the same functional form is used but the suppression scale is the
+dijet invariant mass :math:`m_{jj}`.
+
+
+Dataset routing
+---------------
+
+Each dataset is mapped to one or more power correction parameter keys via the
+function ``get_pc_type``. The mapping depends on the process type and dataset
+name:
+
+.. list-table::
+   :header-rows: 1
+   :widths: 30 30 40
+
+   * - Process type
+     - PC type key
+     - Datasets
+   * - DIS NC (proton :math:`F_2`)
+     - ``f2p``
+     - SLAC, BCDMS proton :math:`F_2`; NMC, HERA :math:`\sigma_{\mathrm{red}}`
+   * - DIS NC (deuteron :math:`F_2`)
+     - ``f2d``
+     - SLAC, BCDMS deuteron :math:`F_2`
+   * - DIS NC (NMC ratio :math:`F_2^d / F_2^p`)
+     - ``(f2p, f2d)``
+     - NMC ratio dataset
+   * - DIS CC
+     - ``dis_cc``
+     - CHORUS, NuTeV, HERA CC
+   * - Jets
+     - ``Hj``
+     - Single-inclusive jet datasets
+   * - Dijets (ATLAS)
+     - ``H2j_ATLAS``
+     - ATLAS dijet datasets (falls back to ``H2j`` if key absent)
+   * - Dijets (CMS)
+     - ``H2j_CMS``
+     - CMS dijet datasets (falls back to ``H2j`` if key absent)
+
+
+Special case: NMC ratio
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+The NMC ratio dataset :math:`F_2^d / F_2^p` receives contributions from both
+the proton and deuteron power corrections. The corrected ratio is
+
+.. math::
+
+   \frac{F_2^d}{F_2^p} \to \frac{F_2^d \,(1 + \mathrm{PC}_d)}{F_2^p \,(1 + \mathrm{PC}_p)},
+
+and the shift is
+
+.. math::
+
+   \Delta\!\left(\frac{F_2^d}{F_2^p}\right) = \frac{F_2^d}{F_2^p} \,
+   \frac{\mathrm{PC}_d - \mathrm{PC}_p}{1 + \mathrm{PC}_p}.
+
+
+Covariance matrix construction
+------------------------------
+
+The theory covariance matrix is constructed from the shifts :math:`\Delta O` by
+taking outer products. For each combination of power correction parameters, a
+shift vector is computed per dataset. The sub-matrix between datasets :math:`i`
+and :math:`j` is then
+
+.. math::
+
+   S_{ij} = \sum_k \Delta_i^{(k)} \otimes \Delta_j^{(k)},
+
+where :math:`k` runs over all parameter combinations (one non-zero ``yshift``
+entry at a time, with all others set to zero). This corresponds to the
+``covmat_power_corrections`` function in ``construction.py``.
+
+
+Runcard configuration
+---------------------
+
+Power corrections are included via the ``theorycovmatconfig`` section of the
+runcard. The key ``"power corrections"`` must be added to the
+``point_prescriptions`` list, alongside any scale variation prescriptions.
+
+The following keys are used:
+
+- ``pc_parameters``: a dictionary mapping PC type keys to their parametrisation
+  (``yshift`` and ``nodes`` arrays). The length of ``yshift`` must match the
+  length of ``nodes``.
+- ``pc_included_procs``: list of process types to which power corrections
+  are applied (e.g. ``["DIS NC", "DIS CC", "JETS", "DIJET"]``).
+- ``pc_excluded_datasets``: list of dataset names to exclude from power corrections
+  even if their process type is included.
+- ``pdf``: the PDF used for computing the theory predictions that enter the
+  multiplicative shifts.
+
+Example
+~~~~~~~
+
+.. code:: yaml
+
+    theorycovmatconfig:
+      point_prescriptions: ["9 point", "power corrections"]
+      pc_parameters:
+        f2p:
+          yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0]
+          nodes:  [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0]
+        f2d:
+          yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0]
+          nodes:  [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0]
+        dis_cc:
+          yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0]
+          nodes:  [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0]
+        Hj:
+          yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
+          nodes:  [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]
+        H2j_ATLAS:
+          yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
+          nodes:  [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]
+        H2j_CMS:
+          yshift: [2.0, 2.0, 2.0, 2.0, 2.0]
+          nodes:  [0.25, 0.75, 1.25, 1.75, 2.25]
+      pc_included_procs: ["JETS", "DIJET", "DIS NC", "DIS CC"]
+      pc_excluded_datasets:
+        - HERA_NC_318GEV_EAVG_CHARM-SIGMARED
+        - HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED
+      pdf: NNPDF40_nnlo_as_01180
+      use_thcovmat_in_fitting: true
+      use_thcovmat_in_sampling: true
+
+.. warning::
+   The lengths of ``yshift`` and ``nodes`` must be equal for each PC type.
+   A mismatch will raise an error at initialisation time.
+
+.. note::
+   Power corrections can be combined with scale variation prescriptions.
+   Both contributions are summed into a single theory covariance matrix.
+   See the tutorial on :ref:`including a theory covmat in a fit <thcov_tutorial>`.
+
+
+Module reference
+----------------
+
+``higher_twist_functions.py`` provides the following public functions:
+
+- ``get_pc_type(exp_name, process_type, experiment, pc_dict)``:
+  determines which PC type key(s) apply to a given dataset.
+- ``linear_bin_function(a, y_shift, bin_edges)``:
+  evaluates the piecewise-linear triangular interpolation at points ``a``.
+- ``dis_pc_func(delta_h, nodes, x, Q2)``:
+  computes the DIS power correction :math:`h(x)/Q^2`.
+- ``jets_pc_func(delta_h, nodes, pT, rap)``:
+  computes the jet power correction :math:`h(y)/p_T`.
+- ``mult_dis_pc(nodes, x, q2, dataset_sp, pdf)``:
+  returns a function that computes the multiplicative DIS shift given node values.
+- ``mult_dis_ratio_pc(p_nodes, d_nodes, x, q2, dataset_sp, pdf)``:
+  returns a function that computes the shift for the :math:`F_2^d/F_2^p` ratio.
+- ``mult_jet_pc(nodes, pT, rap, dataset_sp, pdf)``:
+  returns a function that computes the multiplicative jet shift given node values.
+- ``construct_pars_combs(parameters_dict)``:
+  builds the list of one-at-a-time parameter combinations used to construct
+  the covariance matrix.
+- ``compute_deltas_pc(dataset_sp, pdf, pc_dict)``:
+  computes the full set of shifts for a single dataset.
+
+``construction.py`` provides:
+
+- ``covmat_power_corrections(deltas1, deltas2)``:
+  computes the theory covariance sub-matrix between two datasets from their
+  shift dictionaries.
+- ``covs_pt_prescrip_pc(combine_by_type, point_prescription, pdf, pc_parameters, pc_included_procs, pc_excluded_datasets)``:
+  assembles the full power correction covariance matrix across all datasets.

n3fit/runcards/examples/Basic_runcard_pc_covmat.yml

@@ -0,0 +1,113 @@
+#
+# Configuration file for n3fit
+#
+######################################################################################
+description: NNPDF4.1 with TCM higher twists and jet power corrections
+
+######################################################################################
+dataset_inputs:
+- {dataset: NMC_NC_NOTFIXED_EM-F2, variant: legacy_dw}
+- {dataset: NMC_NC_NOTFIXED_P_EM-SIGMARED, variant: legacy}
+- {dataset: SLAC_NC_NOTFIXED_P_EM-F2, variant: legacy_dw}
+- {dataset: SLAC_NC_NOTFIXED_D_EM-F2, variant: legacy_dw}
+- {dataset: BCDMS_NC_NOTFIXED_P_EM-F2, variant: legacy_dw}
+- {dataset: BCDMS_NC_NOTFIXED_D_EM-F2, variant: legacy_dw}
+- {dataset: CHORUS_CC_NOTFIXED_PB_NU-SIGMARED, variant: legacy_dw}
+- {dataset: CHORUS_CC_NOTFIXED_PB_NB-SIGMARED, variant: legacy_dw}
+- {dataset: NUTEV_CC_NOTFIXED_FE_NU-SIGMARED, cfac: [MAS], variant: legacy_dw}
+- {dataset: NUTEV_CC_NOTFIXED_FE_NB-SIGMARED, cfac: [MAS], variant: legacy_dw}
+- {dataset: HERA_NC_318GEV_EM-SIGMARED}
+- {dataset: ATLAS_1JET_8TEV_R06_PTY, variant: decorrelated}
+- {dataset: ATLAS_2JET_7TEV_R06_M12Y}
+- {dataset: CMS_1JET_8TEV_PTY}
+- {dataset: CMS_2JET_7TEV_M12-Y}
+- {dataset: CMS_2JET_13TEV_M12-YSTAR-YB-R08}
+
+################################################################################
+diagonal_frac: 0.75
+
+datacuts:
+  t0pdfset: 260202-jk-nnpdf41-mhou
+  q2min: 2.5
+  w2min: 3.24
+
+theory:
+  theoryid: 41_000_000
+
+theorycovmatconfig:
+  point_prescriptions: ["power corrections"]
+  pc_parameters:
+    f2p:        {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
+    f2d:        {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
+    dis_cc:     {yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
+    Hj:         {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]}
+    H2j_ystar:  {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]}
+    H2j_yb:   {yshift: [2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25]}
+  pc_included_procs: ["DIS NC", "DIS CC", "JETS", "DIJET"]
+  pc_excluded_datasets: []
+  pdf: 260202-jk-nnpdf41-mhou
+  use_thcovmat_in_fitting: true
+  use_thcovmat_in_sampling: true
+resample_negative_pseudodata: false
+
+trvlseed: 130582403
+nnseed: 953262798
+mcseed: 1437981271
+genrep: true
+parameters: # This defines the parameter dictionary that is passed to the Model Trainer
+  nodes_per_layer: [25, 20, 9]
+  activation_per_layer: [tanh, tanh, linear]
+  initializer: glorot_normal
+  optimizer:
+    clipnorm: 6.073e-6
+    learning_rate: 2.621e-3
+    optimizer_name: Nadam
+  epochs: 27000
+  positivity:
+    initial: 184.8
+    multiplier:
+  integrability:
+    initial: 10
+    multiplier:
+  stopping_patience: 0.1
+  layer_type: dense
+  dropout: 0.0
+  threshold_chi2: 3.5
+  feature_scaling_points: 5
+
+fitting:
+  fitbasis: CCBAR_ASYMM  # EVOL (7), EVOLQED (8), etc.
+  savepseudodata: true
+  basis:
+  - {fl: sng, trainable: false, smallx: [1.058, 1.155]}
+  - {fl: g, trainable: false, smallx: [0.9017, 1.084]}
+  - {fl: v, trainable: false, smallx: [0.481, 0.6499]}
+  - {fl: v3, trainable: false, smallx: [0.08225, 0.502]}
+  - {fl: v8, trainable: false, smallx: [0.5823, 0.7928]}
+  - {fl: t3, trainable: false, smallx: [-0.3987, 0.9689]}
+  - {fl: t8, trainable: false, smallx: [0.6077, 0.9459]}
+  - {fl: t15, trainable: false, smallx: [1.023, 1.147]}
+  - {fl: v15, trainable: false, smallx: [0.5005, 0.7189]}
+
+################################################################################
+positivity:
+  posdatasets:
+  # Positivity of MSbar PDFs
+  - {dataset: NNPDF_POS_100GEV_XUQ, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XUB, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XDQ, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XDB, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XSQ, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XSB, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XCQ, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XCB, maxlambda: 1e6}
+  - {dataset: NNPDF_POS_100GEV_XGL, maxlambda: 1e6}
+
+integrability:
+  integdatasets:
+  - {dataset: NNPDF_INTEG_3GEV_XT8, maxlambda: 1e2}
+  - {dataset: NNPDF_INTEG_3GEV_XT3, maxlambda: 1e2}
+
+################################################################################
+debug: false
+maxcores: 16

validphys2/src/validphys/checks.py

@@ -361,3 +361,15 @@ def check_darwin_single_process(NPROC):
     """
     if platform.system() == "Darwin" and NPROC != 1:
         raise CheckError("NPROC must be set to 1 on OSX, because multithreading is not supported.")
+
+
+@make_argcheck
+def check_pc_parameters(pc_parameters):
+    """Check that the parameters for the PC method are set correctly"""
+    for par in pc_parameters.values():
+        # Check that the length of shifts is the same as the length of nodes
+        if len(par['yshift']) != len(par['nodes']):
+            raise ValueError(
+                f"The length of nodes does not match that of the list in {par['ht']}."
+                f"Check the runcard. Got {len(par['yshift'])} != {len(par['nodes'])}"
+            )