Add final-state QED corrections

.github/workflows/release.yml

@@ -0,0 +1,38 @@
+name: Release
+
+on:
+  push:
+    tags:
+      - 'v[0-9]+*'
+  workflow_dispatch:
+
+env:
+  GH_TOKEN: ${{ github.token }}
+
+jobs:
+  create-release:
+    runs-on: ubuntu-latest
+    if: "startsWith(github.ref, 'refs/tags/')"
+    steps:
+      - uses: actions/checkout@v4
+      - name: Create new release
+        run: |
+          gh release create v${GITHUB_REF_NAME#v} --generate-notes
+
+  publish-crates:
+    runs-on: ubuntu-latest
+    if: "startsWith(github.ref, 'refs/tags/')"
+    steps:
+      - uses: actions/checkout@v4
+      - name: Install Rust toolchain
+        uses: dtolnay/rust-toolchain@stable
+      - name: Publish all crates
+        env:
+          CARGO_REGISTRY_TOKEN: ${{ secrets.CARGO_REGISTRY_TOKEN }}
+        run: |
+          cargo publish -p navis-core
+          sleep 30
+          cargo publish -p navis-unpol
+          cargo publish -p navis-pol
+          sleep 30
+          cargo publish -p navis-cli

Cargo.toml

@@ -11,6 +11,7 @@ members = [
 version = "0.1.0"
 edition = "2021"
 rust-version = "1.91.0"
+repository = "https://github.com/QCDLab/navis"
 license = "MIT"
 
 [workspace.dependencies]

README.md

@@ -3,19 +3,22 @@
 <p align="justify">
   <b>ηaνis</b> is a code that computes both polarized and unpolarized single-inclusive
   hadron production in proton-proton collision up to next-to-leading order (NLO) in
-  QCD and dump the results in the form of
+  QCD⊗QED and dump the results in the form of
   <a href="https://github.com/NNPDF/pineappl">PineAPPL</a>
   fast interpolation grids. It is based on the following papers:
-  <a href="https://arxiv.org/abs/hep-ph/0210442">[arXiv:0210442]</a> and
-  <a href="https://arxiv.org/abs/hep-ph/0211007">[arXiv:0211007]</a>.
+  <a href="https://arxiv.org/abs/hep-ph/0210442">[arXiv:0210442]</a>,
+  <a href="https://arxiv.org/abs/hep-ph/0211007">[arXiv:0211007]</a>,
+  <a href="https://bib-pubdb1.desy.de/record/517576/files/">[DESY94-042]</a>,
+  and <a href="https://inspirehep.net/literature/342979">[iNSPIRE:342979]</a>.
 </p>
 
 ## Installation
 
 To compile the code, you only need [cargo](https://doc.rust-lang.org/cargo/) and then
-run the following command:
+run the following commands:
 
 ```sh
+git clone https://github.com/QCDLab/navis.git
 cargo install --path navis-cli
 ```
 
@@ -28,7 +31,7 @@ navis-unpol ./runcards/example.yml
 
 ## Citation
 
-If you use **ηaνis**, please cite the following:
+If you use **ηaνis**, please cite the papers above and the following:
 
 <table>
   <tr>

maintainer/make-release.sh

@@ -0,0 +1,73 @@
+#!/bin/bash
+
+set -euo pipefail
+
+crates=(
+    navis-core
+    navis-unpol
+    navis-pol
+    navis-cli
+)
+
+main=main
+root_dir=$(cd "$(dirname "$0")/.." && pwd)
+cd "${root_dir}"
+
+this_branch=$(git rev-parse --abbrev-ref HEAD)
+
+if [[ $# != 1 ]]; then
+    echo "No version number given."
+    exit 1
+fi
+
+version=$1
+
+prerelease=$(echo ${version} | perl -pe 's/^(0|[1-9]\d*)\.(0|[1-9]\d*)\.(0|[1-9]\d*)(?:-((?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*)(?:\.(?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*))*))?(?:\+([0-9a-zA-Z-]+(?:\.[0-9a-zA-Z-]+)*))?$/\4/')
+
+if [[ $(echo ${version} | grep -oP '^(0|[1-9]\d*)\.(0|[1-9]\d*)\.(0|[1-9]\d*)(?:-((?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*)(?:\.(?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*))*))?(?:\+([0-9a-zA-Z-]+(?:\.[0-9a-zA-Z-]+)*))?$') != ${version} ]]; then
+    echo "Version string incorrect."
+    exit 1
+fi
+
+if [[ $(git tag -l v${version}) ]]; then
+    echo "Version already exists."
+    exit 1
+fi
+
+# in branches that are not main we only allow prereleases
+if [[ ${this_branch} != ${main} ]] && [[ ${prerelease} == "" ]]; then
+    echo "Ordinary releases are only allowed in the '${main}' branch."
+    echo "If you really want to make a release from '${this_branch}', consider making a prerelease."
+    exit 1
+fi
+
+echo ">>> Updating version strings ..."
+
+# modify the version in the main workspace
+sed -i \
+    -e "/^\[workspace\.package\]/,/^\[/ s:^version = \".*\":version = \"${version}\":" \
+    Cargo.toml
+
+# the internal crates are only depended upon via path in
+# `[workspace.dependencies]`; `cargo publish` refuses path dependencies that
+# don't also carry a version, so pin one here for whichever depends on them
+for crate in navis-core navis-pol navis-unpol; do
+    sed -i \
+        -e "s:${crate} = { path = \"\./${crate}\" }:${crate} = { path = \"./${crate}\", version = \"${version}\" }:" \
+        -e "s:${crate} = { path = \"\./${crate}\", version = \"[^\"]*\" }:${crate} = { path = \"./${crate}\", version = \"${version}\" }:" \
+        Cargo.toml
+done
+git add Cargo.toml
+
+echo ">>> Updating Cargo.lock ..."
+
+for crate in "${crates[@]}"; do
+    cargo pkgid path+file://$(cd "${crate}" && pwd)
+done | xargs printf " -p %s" | xargs cargo update
+git add Cargo.lock
+
+echo ">>> Committing and pushing changes ..."
+
+git commit -m "Release v${version}"
+git tag -a v${version} -m v${version}
+git push --follow-tags

navis-cli/Cargo.toml

@@ -3,7 +3,8 @@ name = "navis-cli"
 version.workspace = true
 edition.workspace = true
 rust-version.workspace = true
-description = "Command-line drivers for unpolarized and polarized single-inclusive hadron production"
+description = "CLI interface"
+repository.workspace = true
 license.workspace = true
 
 [[bin]]

navis-cli/src/bin/navis-pol.rs

@@ -43,6 +43,7 @@ fn main() -> anyhow::Result<()> {
         isigm: CrossSectionType::from_isigm(card.isigm)?,
         hc2: params.hc2,
         pi: params.pi,
+        q2pho: card.q2pho,
     };
     let targets = Targets {
         ih1: card.ih1,
@@ -64,7 +65,7 @@ fn main() -> anyhow::Result<()> {
     let edges = card.bin_edges();
     let grid_config = card
         .generate_grids
-        .then(|| GridConfig::new(order, convolutions, &edges));
+        .then(|| GridConfig::new(order, convolutions, &edges, false));
 
     let n_iter = card.iter_max as usize;
     let n_eval = card.ncalls_vegas as usize;

navis-cli/src/bin/navis-unpol.rs

@@ -43,6 +43,7 @@ fn main() -> anyhow::Result<()> {
         isigm: CrossSectionType::from_isigm(card.isigm)?,
         hc2: params.hc2,
         pi: params.pi,
+        q2pho: card.q2pho,
     };
     let targets = Targets {
         ih1: card.ih1,
@@ -64,7 +65,7 @@ fn main() -> anyhow::Result<()> {
     let edges = card.bin_edges();
     let grid_config = card
         .generate_grids
-        .then(|| GridConfig::new(order, convolutions, &edges));
+        .then(|| GridConfig::new(order, convolutions, &edges, true));
 
     let n_iter = card.iter_max as usize;
     let n_eval = card.ncalls_vegas as usize;

navis-core/Cargo.toml

@@ -3,7 +3,8 @@ name = "navis-core"
 version.workspace = true
 edition.workspace = true
 rust-version.workspace = true
-description = "Shared kinematics, VEGAS integrator, PineAPPL grid setup, and runcard parsing for navis"
+description = "navis core API"
+repository.workspace = true
 license.workspace = true
 
 [dependencies]

navis-core/src/constants.rs

@@ -4,8 +4,10 @@
 pub const CF: f64 = 4.0 / 3.0;
 /// Number of colors `C_A = N_C = 3`.
 pub const CA: f64 = 3.0;
-/// Conversion factor `(hbar*c)^2` in GeV^2 pbarn, i.e. `HC2` in the Fortran code.
+/// Conversion factor `(hbar*c)^2` in GeV^2 pbarn.
 pub const HC2: f64 = 0.389_429_57e9;
+/// QED fine-structure constant `AEM`.
+pub const AEM: f64 = 1.0 / 137.036;
 
 /// Physics parameters derived from the runcard, computed once per run.
 #[derive(Debug, Clone, Copy)]
@@ -15,12 +17,12 @@ pub struct Params {
     pub nf: f64,
     pub pi: f64,
     pub hc2: f64,
+    pub aem: f64,
     /// `VC = C_A^2 - 1`.
     pub vc: f64,
     /// `GTR = NF / 2`.
     pub gtr: f64,
-    /// The 16-entry channel normalization prefactor `CC(J0)`, see [`prefactors`].
-    pub cc: [f64; 16],
+    pub cc: [f64; 21],
 }
 
 impl Params {
@@ -32,25 +34,26 @@ impl Params {
             nf,
             pi: std::f64::consts::PI,
             hc2: HC2,
+            aem: AEM,
             vc: CA * CA - 1.0,
             gtr: nf / 2.0,
             cc: prefactors(CA),
         }
     }
 }
 
-/// Per-channel color prefactor `CC(J0)`, `J0 = 1..16`, computed in the `DO 2
-/// J0=1,16`.
-///
-/// Channels 15/16 get `8*VC^2`; channels 8,9,10,13,14 get `8*VC*NC`; all
-/// others get `8*NC^2`.
+/// Per-channel color prefactor `CC(J0)`. Channels 15/16 get `8*VC^2`;
+/// channels 8,9,10,13,14 get `8*VC*NC`; channels 17-21 (photon) get `1`;
+/// all others get `8*NC^2`.
 #[must_use]
-pub fn prefactors(nc: f64) -> [f64; 16] {
+pub fn prefactors(nc: f64) -> [f64; 21] {
     let vc = nc * nc - 1.0;
-    let mut cc = [0.0; 16];
+    let mut cc = [0.0; 21];
     for (j0, entry) in cc.iter_mut().enumerate() {
         let channel = j0 + 1;
-        *entry = if channel == 16 || channel == 15 {
+        *entry = if channel >= 17 {
+            1.0
+        } else if channel == 16 || channel == 15 {
             8.0 * vc * vc
         } else if matches!(channel, 14 | 13 | 10 | 9 | 8) {
             8.0 * vc * nc

navis-core/src/grid_setup.rs

@@ -75,7 +75,12 @@ impl GridConfig {
     /// 1`) across *all* pT bins in the run: the Q2 interpolation range is
     /// derived from the overall pT range.
     #[must_use]
-    pub fn new(order: PertOrder, convolutions: Vec<Conv>, global_edges: &[f64]) -> Self {
+    pub fn new(
+        order: PertOrder,
+        convolutions: Vec<Conv>,
+        global_edges: &[f64],
+        include_qed: bool,
+    ) -> Self {
         let bin_min = global_edges[0];
         let bin_max = global_edges[global_edges.len() - 1];
 
@@ -119,7 +124,7 @@ impl GridConfig {
         ];
 
         Self {
-            channels: sihp_pdg_channels(),
+            channels: sihp_pdg_channels(include_qed),
             orders: order.orders(),
             convolutions,
             interps,

navis-core/src/kinematics.rs

@@ -58,14 +58,14 @@ pub struct BinInputs {
 #[derive(Debug, Clone, Copy)]
 pub struct RunParams {
     pub sqs: f64,
-    /// `S = SQS**2`.
     pub s: f64,
     pub scfac: f64,
     pub scmu: f64,
     pub scfrag: f64,
     pub isigm: CrossSectionType,
     pub hc2: f64,
     pub pi: f64,
+    pub q2pho: f64,
 }
 
 /// The result of mapping one VEGAS point through the 5D phase space.
@@ -113,8 +113,7 @@ impl PhaseSpace {
     }
 }
 
-/// Map a 5D VEGAS point `xx` (`XX(1..5)` in Fortran, 0-indexed here) to a
-/// [`PhaseSpace`] point.
+/// Map a 5D VEGAS point `xx` to a [`PhaseSpace`] point.
 #[must_use]
 pub fn map_phase_space(
     xx: &[f64; 5],