diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index b9194d68c..c51cae83b 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -37,11 +37,7 @@ class BenchmarkVFNS(ApfelBenchmark): vfns_theory = { "FNS": "ZM-VFNS", - "ModEv": [ - "EXA", - "EXP", - "TRN", - ], + "ModEv": "EXA", "kcThr": 1.0, "kbThr": 1.0, "ktThr": 1.0, @@ -132,12 +128,10 @@ class BenchmarkFFNS(ApfelBenchmark): ffns_theory = { "FNS": "FFNS", - "ModEv": [ - "EXA", - "EXP", - "TRN", - ], + "ModEv": "EXA", "NfFF": 4, + "nf0": 4, + "nfref": 4, "kcThr": 0.0, "kbThr": np.inf, "ktThr": np.inf, @@ -159,7 +153,6 @@ def benchmark_plain_pol(self, pto): th = self.ffns_theory.copy() th.update({"PTO": [pto]}) - th["ModEv"] = ["EXA"] # TODO for the time one is sufficient op = operators.apfel_config.copy() op["polarized"] = [True] self.run(cartesian_product(th), operators.build(op), ["ToyLH"]) @@ -209,6 +202,8 @@ class BenchmarkFFNS_qed(ApfelBenchmark): # "TRN", ], "NfFF": 5, + "nf0": 5, + "nfref": 5, "kcThr": 0.0, "kbThr": 0.0, "ktThr": np.inf, @@ -291,6 +286,8 @@ class BenchmarkVFNS_qed(ApfelBenchmark): "kbThr": 1.0, "ktThr": 1.0, "Q0": 1.25, + "nf0": 3, + "nfref": 5, "alphas": 0.118000, "alphaqed": 0.007496, } @@ -319,7 +316,7 @@ def benchmark_plain(self, pto, qed): # obj = BenchmarkFFNS() # obj.benchmark_plain_pol(2) - # obj.benchmark_plain(2) + # obj.benchmark_plain(0) # obj.benchmark_sv(2, "exponentiated") # obj.benchmark_kthr(2) # obj.benchmark_msbar(2) diff --git a/benchmarks/sandbox.py b/benchmarks/sandbox.py index 0a9a04e48..9fd576950 100644 --- a/benchmarks/sandbox.py +++ b/benchmarks/sandbox.py @@ -44,20 +44,20 @@ class Sandbox(Runner): sandbox = True # select here the external program between LHA, LHAPDF, apfel, pegasus - external = "apfel" + external = "void" # external = "pegasus" # select to plot operators - plot_operator = True + # plot_operator = True - rotate_to_evolution_basis = True + # rotate_to_evolution_basis = True @staticmethod def generate_operators(): ops = { "ev_op_iterations": [1], # "ev_op_max_order": [20], - "mugrid": [[10]], + "mugrid": [[10.0]], # "debug_skip_singlet": [True], } return ops @@ -118,12 +118,12 @@ def lha(self): "FNS": "FFNS", "NfFF": 4, "ModEv": "EXA", - "XIF": np.sqrt(2), - "Q0": np.sqrt(2), + "XIF": float(np.sqrt(2)), + "Q0": float(np.sqrt(2)), "kcThr": 0.0, "kbThr": np.inf, "ktThr": np.inf, - "Qref": np.sqrt(2.0), + "Qref": float(np.sqrt(2.0)), "alphas": 0.35, } self.skip_pdfs = lambda _theory: [ @@ -139,11 +139,11 @@ def lha(self): ] self.run( [theory_updates], - [{"mugrid": [100], "debug_skip_singlet": True}], + [{"mugrid": [100.0], "debug_skip_singlet": True}], ["ToyLH"], ) if __name__ == "__main__": sand = Sandbox() - sand.doit() + sand.lha() diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index a6c5fa9ae..b79cc48b2 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -7,13 +7,15 @@ use std::ffi::c_void; pub mod bib; pub mod mellin; -/// Wrapper to pass arguments back to Python +/// Wrapper to pass arguments back to Python. struct RawCmplx { re: Vec, im: Vec, } -/// Map tensors to c-ordered list +/// Map tensor with shape (o,d,d) to c-ordered list. +/// +/// This is needed for the QCD singlet. fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: usize) -> RawCmplx { let mut target = RawCmplx { re: Vec::::new(), @@ -30,14 +32,63 @@ fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: us target } -/// QCD intergration kernel inside quad. +/// Map tensor with shape (o,o',d,d) to c-ordered list. +/// +/// This is needed for the QED singlet and valence. +fn unravel_qed( + res: Vec; DIM]; DIM]>>, + order_qcd: usize, + order_qed: usize, +) -> RawCmplx { + let mut target = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + for obj_ in res.iter().take(order_qcd + 1) { + for obj in obj_.iter().take(order_qed + 1) { + for col in obj.iter().take(DIM) { + for el in col.iter().take(DIM) { + target.re.push(el.re); + target.im.push(el.im); + } + } + } + } + target +} + +/// Map tensor with shape (o,o',d) to c-ordered list. +/// +/// This is needed for the QED non-singlet. +fn unravel_qed_ns(res: Vec>>, order_qcd: usize, order_qed: usize) -> RawCmplx { + let mut target = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + for col in res.iter().take(order_qcd + 1) { + for el in col.iter().take(order_qed + 1) { + target.re.push(el.re); + target.im.push(el.im); + } + } + target +} + +/// Intergration kernel inside quad. /// /// # Safety /// This is the connection from Python, so we don't know what is on the other side. #[no_mangle] -pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { - let args = *(rargs as *mut QuadQCDargs); - let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0); +pub unsafe extern "C" fn rust_quad_ker(u: f64, rargs: *mut c_void) -> f64 { + let args = *(rargs as *mut QuadArgs); + + let is_singlet = (100 == args.mode0) + || (21 == args.mode0) + || (90 == args.mode0) + || (22 == args.mode0) + || (101 == args.mode0); + + let is_qed_valence = (10200 == args.mode0) || (10204 == args.mode0); // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); let jac = path.jac() * path.prefactor(); @@ -70,14 +121,41 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { ); } } else if is_singlet { - let gamma_singlet_qcd = match args.is_polarized { - true => ekore::anomalous_dimensions::polarized::spacelike::gamma_singlet_qcd, - false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd, - }; - raw = unravel( - gamma_singlet_qcd(args.order_qcd, &mut c, args.nf), - args.order_qcd, - ); + if args.order_qed > 0 { + let gamma_singlet_qed = + ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qed; + raw = unravel_qed( + gamma_singlet_qed(args.order_qcd, args.order_qed, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } else { + let gamma_singlet_qcd = match args.is_polarized { + true => ekore::anomalous_dimensions::polarized::spacelike::gamma_singlet_qcd, + false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd, + }; + raw = unravel( + gamma_singlet_qcd(args.order_qcd, &mut c, args.nf), + args.order_qcd, + ); + } + } else if args.order_qed > 0 { + if is_qed_valence { + let gamma_valence_qed = + ekore::anomalous_dimensions::unpolarized::spacelike::gamma_valence_qed; + raw = unravel_qed( + gamma_valence_qed(args.order_qcd, args.order_qed, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } else { + let gamma_ns_qed = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qed; + raw = unravel_qed_ns( + gamma_ns_qed(args.order_qcd, args.order_qed, args.mode0, &mut c, args.nf), + args.order_qcd, + args.order_qed, + ); + } } else { // we can not do 1D let gamma_ns_qcd = match args.is_polarized { @@ -100,6 +178,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { jac.re, jac.im, args.order_qcd, + args.order_qed, is_singlet, args.mode0, args.mode1, @@ -109,7 +188,6 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.areas, args.areas_x, args.areas_y, - args.L, args.method_num, args.as1, args.as0, @@ -118,50 +196,68 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.sv_mode_num, args.is_threshold, args.Lsv, + // additional QED params + args.as_list, + args.as_list_len, + args.mu2_from, + args.mu2_to, + args.a_half, + args.a_half_x, + args.a_half_y, + args.alphaem_running, ) } /// Python callback signature -type PyQuadKerQCDT = unsafe extern "C" fn( - *const f64, - *const f64, - f64, - f64, - f64, - f64, - usize, - bool, - u16, - u16, - u8, - bool, - f64, - *const f64, - u8, - u8, - f64, - u8, - f64, - f64, - u8, - u8, - u8, - bool, - f64, +type PyQuadKerT = unsafe extern "C" fn( + *const f64, // re_gamma + *const f64, // im_gamma + f64, // re_n + f64, // im_n + f64, // re_jac + f64, // im_jac + usize, // order_qcd + usize, // order_qed + bool, // is_singlet + u16, // mode0 + u16, // mode1 + u8, // nf + bool, // is_log + f64, // logx + *const f64, // areas + u8, // areas_x + u8, // areas_y + u8, // method_num + f64, // as1 + f64, // as0 + u8, // ev_op_iterations + u8, // ev_op_max_order_qcd + u8, // sv_mode_num + bool, // is_threshold + f64, // lsv + *const f64, // as_list + u8, // as_list_len + f64, // mu2_from + f64, // mu2_to + *const f64, // a_half + u8, // a_half_x + u8, // a_half_y + bool, // alphaem_running ) -> f64; /// Additional integration parameters #[allow(non_snake_case)] #[repr(C)] #[derive(Clone, Copy)] -pub struct QuadQCDargs { +pub struct QuadArgs { pub order_qcd: usize, + pub order_qed: usize, pub mode0: u16, pub mode1: u16, pub is_polarized: bool, pub is_time_like: bool, pub nf: u8, - pub py: PyQuadKerQCDT, + pub py: PyQuadKerT, pub is_log: bool, pub logx: f64, pub areas: *const f64, @@ -177,6 +273,15 @@ pub struct QuadQCDargs { pub is_threshold: bool, pub is_ome: bool, pub Lsv: f64, + // additional param required for QED + pub as_list: *const f64, + pub as_list_len: u8, + pub mu2_from: f64, + pub mu2_to: f64, + pub a_half: *const f64, + pub a_half_x: u8, + pub a_half_y: u8, + pub alphaem_running: bool, } /// Empty placeholder function for python callback. @@ -191,6 +296,7 @@ pub unsafe extern "C" fn my_py( _re_jac: f64, _im_jac: f64, _order_qcd: usize, + _order_qed: usize, _is_singlet: bool, _mode0: u16, _mode1: u16, @@ -200,7 +306,6 @@ pub unsafe extern "C" fn my_py( _areas: *const f64, _areas_x: u8, _areas_y: u8, - _l: f64, _method_num: u8, _as1: f64, _as0: f64, @@ -209,6 +314,14 @@ pub unsafe extern "C" fn my_py( _sv_mode_num: u8, _is_threshold: bool, _lsv: f64, + _as_list: *const f64, + _as_list_len: u8, + _mu2_from: f64, + _mu2_to: f64, + _a_half: *const f64, + _a_half_x: u8, + _a_half_y: u8, + _alphaem_running: bool, ) -> f64 { 0. } @@ -216,14 +329,15 @@ pub unsafe extern "C" fn my_py( /// Return empty additional arguments. /// /// This is required to make the arguments part of the API, otherwise it won't be added to the compiled -/// package (since it does not appear in the signature of `rust_quad_ker_qcd`). +/// package (since it does not appear in the signature of `rust_quad_ker`). /// /// # Safety /// This is the connection from and back to Python, so we don't know what is on the other side. #[no_mangle] -pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { - QuadQCDargs { +pub unsafe extern "C" fn empty_args() -> QuadArgs { + QuadArgs { order_qcd: 0, + order_qed: 0, mode0: 0, mode1: 0, is_polarized: false, @@ -245,5 +359,13 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { is_threshold: false, is_ome: false, Lsv: 0., + as_list: [].as_ptr(), + as_list_len: 0, + mu2_from: 0., + mu2_to: 0., + a_half: [].as_ptr(), + a_half_x: 0, + a_half_y: 0, + alphaem_running: false, } } diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index 0035b7a39..13d85e8b5 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -107,3 +107,52 @@ @article{Bierenbaum2009zt pages = "401--406", year = "2009" } +@article{Gluck1995yr, + author = "Gluck, M. and Reya, E. and Stratmann, M. and Vogelsang, W.", + title = "{Next-to-leading order radiative parton model analysis of polarized deep inelastic lepton - nucleon scattering}", + eprint = "hep-ph/9508347", + archivePrefix = "arXiv", + reportNumber = "DO-TH-95-13, RAL-TR-95-042", + doi = "10.1103/PhysRevD.53.4775", + journal = "Phys. Rev. D", + volume = "53", + pages = "4775--4786", + year = "1996" +} +@phdthesis{Carrazza2015dea, + author = "Carrazza, Stefano", + title = "{Parton distribution functions with QED corrections}", + eprint = "1509.00209", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + doi = "10.13130/carrazza-stefano_phd2015-07-06", + school = "Milan U.", + year = "2015" +} +@article{deFlorian2015ujt, + author = "de Florian, Daniel and Sborlini, Germ\'an F. R. and Rodrigo, Germ\'an", + title = "{QED corrections to the Altarelli\textendash{}Parisi splitting functions}", + eprint = "1512.00612", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "ICAS-03-15, IFIC-15-81", + doi = "10.1140/epjc/s10052-016-4131-8", + journal = "Eur. Phys. J. C", + volume = "76", + number = "5", + pages = "282", + year = "2016" +} +@article{deFlorian2016gvk, + author = "de Florian, Daniel and Sborlini, Germ\'an F. R. and Rodrigo, Germ\'an", + title = "{Two-loop QED corrections to the Altarelli-Parisi splitting functions}", + eprint = "1606.02887", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "ICAS-09-16, IFIC-15-88", + doi = "10.1007/JHEP10(2016)056", + journal = "JHEP", + volume = "10", + pages = "056", + year = "2016" +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs index dec73bf04..08de553fe 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs @@ -1,10 +1,15 @@ //! The unpolarized, space-like anomalous dimensions at various couplings power. -use crate::constants::{PID_NSM, PID_NSP, PID_NSV}; +use crate::constants::{ + ED2, EU2, PID_NSM, PID_NSM_D, PID_NSM_U, PID_NSP, PID_NSP_D, PID_NSP_U, PID_NSV, +}; use crate::harmonics::cache::Cache; use num::complex::Complex; use num::Zero; +pub mod aem1; +pub mod aem2; pub mod as1; +pub mod as1aem1; pub mod as2; pub mod as3; @@ -55,3 +60,285 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl } gamma_S } + +/// Compute the tower of the |QCD| x |QED| non-singlet anomalous dimensions. +pub fn gamma_ns_qed( + order_qcd: usize, + order_qed: usize, + mode: u16, + c: &mut Cache, + nf: u8, +) -> Vec>> { + let col = vec![Complex::::zero(); order_qed + 1]; + let mut gamma_ns = vec![col; order_qcd + 1]; + + // QCD corrections + let qcd_mode = match mode { + PID_NSP_U | PID_NSP_D => PID_NSP, + PID_NSM_U | PID_NSM_D => PID_NSM, + _ => mode, + }; + let gamma_qcd = gamma_ns_qcd(order_qcd, qcd_mode, c, nf); + for j in 0..order_qcd { + gamma_ns[1 + j][0] = gamma_qcd[j]; + } + // QED corrections + gamma_ns[0][1] = match mode { + PID_NSP_U | PID_NSM_U => EU2 * aem1::gamma_ns(c, nf), + PID_NSP_D | PID_NSM_D => ED2 * aem1::gamma_ns(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; + if order_qed >= 2 { + gamma_ns[0][2] = match mode { + PID_NSP_U => EU2 * aem2::gamma_nspu(c, nf), + PID_NSP_D => ED2 * aem2::gamma_nspd(c, nf), + PID_NSM_U => EU2 * aem2::gamma_nsmu(c, nf), + PID_NSM_D => ED2 * aem2::gamma_nsmd(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; + } + // QCDxQED corrections + gamma_ns[1][1] = match mode { + PID_NSP_U => EU2 * as1aem1::gamma_nsp(c, nf), + PID_NSP_D => ED2 * as1aem1::gamma_nsp(c, nf), + PID_NSM_U => EU2 * as1aem1::gamma_nsm(c, nf), + PID_NSM_D => ED2 * as1aem1::gamma_nsm(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; + gamma_ns +} + +/// Compute the tower of the |QCD| x |QED| singlet anomalous dimensions matrices. +pub fn gamma_singlet_qed( + order_qcd: usize, + order_qed: usize, + c: &mut Cache, + nf: u8, +) -> Vec; 4]; 4]>> { + let col = vec![ + [[ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ]; 4]; + order_qed + 1 + ]; + let mut gamma_s = vec![col; order_qcd + 1]; + + // QCD corrections + let gamma_qcd_s = gamma_singlet_qcd(order_qcd, c, nf); + let gamma_qcd_nsp = gamma_ns_qcd(order_qcd, PID_NSP, c, nf); + for j in 0..order_qcd { + let s = gamma_qcd_s[j]; + gamma_s[1 + j][0] = [ + [ + s[1][1], + Complex::::zero(), + s[1][0], + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + s[0][1], + Complex::::zero(), + s[0][0], + Complex::::zero(), + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + gamma_qcd_nsp[j], + ], + ]; + } + // QED corrections + gamma_s[0][1] = aem1::gamma_singlet(c, nf); + if order_qed >= 2 { + gamma_s[0][2] = aem2::gamma_singlet(c, nf); + } + // QCDxQED corrections + gamma_s[1][1] = as1aem1::gamma_singlet(c, nf); + gamma_s +} + +/// Compute the tower of the |QCD| x |QED| valence anomalous dimensions matrices. +pub fn gamma_valence_qed( + order_qcd: usize, + order_qed: usize, + c: &mut Cache, + nf: u8, +) -> Vec; 2]; 2]>> { + let col = vec![[[Complex::::zero(), Complex::::zero(),]; 2]; order_qed + 1]; + let mut gamma_v = vec![col; order_qcd + 1]; + + // QCD corrections + let gamma_qcd_nsv = gamma_ns_qcd(order_qcd, PID_NSV, c, nf); + let gamma_qcd_nsm = gamma_ns_qcd(order_qcd, PID_NSM, c, nf); + for j in 0..order_qcd { + gamma_v[1 + j][0] = [ + [gamma_qcd_nsv[j], Complex::::zero()], + [Complex::::zero(), gamma_qcd_nsm[j]], + ]; + } + // QED corrections + gamma_v[0][1] = aem1::gamma_valence(c, nf); + if order_qed >= 2 { + gamma_v[0][2] = aem2::gamma_valence(c, nf); + } + // QCDxQED corrections + gamma_v[1][1] = as1aem1::gamma_valence(c, nf); + gamma_v +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{ + assert_approx_eq_cmplx, assert_approx_eq_cmplx_1d, assert_approx_eq_cmplx_2d, cmplx, + }; + use num::complex::Complex; + + #[test] + fn test_gamma_ns_qcd() { + const NF: u8 = 3; + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + assert_approx_eq_cmplx!( + f64, + gamma_ns_qcd(3, PID_NSP, &mut c, NF)[0], + cmplx!(0., 0.), + epsilon = 1e-14 + ); + + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(2, PID_NSM, &mut c, NF), + [cmplx!(0., 0.); 2], + 2, + epsilon = 2e-6 + ); + + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(3, PID_NSM, &mut c, NF), + [cmplx!(0., 0.); 3], + 3, + epsilon = 2e-4 + ); + + assert_approx_eq_cmplx_1d!( + f64, + gamma_ns_qcd(3, PID_NSV, &mut c, NF), + [cmplx!(0., 0.); 3], + 3, + epsilon = 8e-4 + ); + } + + #[test] + fn test_gamma_ns_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + + // ns- + for pid in [PID_NSM_U, PID_NSM_D] { + assert_approx_eq_cmplx_2d!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF), + [[cmplx!(0., 0.); 2]; 2], + 2, + epsilon = 1e-5 + ); + } + + // ns+ + for pid in [PID_NSP_U, PID_NSP_U] { + // as^0 a^0 must be trivial + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF)[0][0], + cmplx!(0., 0.) + ); + assert_approx_eq_cmplx!( + f64, + gamma_ns_qed(1, 1, pid, &mut c, NF)[0][1], + cmplx!(0., 0.), + epsilon = 1e-5 + ); + } + } + + #[test] + fn test_gamma_valence_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + let g = gamma_valence_qed(3, 2, &mut c, NF); + // as^0 a^0 must be trivial + assert_approx_eq_cmplx_2d!(f64, g[0][0], [[cmplx!(0., 0.); 2]; 2], 2); + // reference from Python side + assert_approx_eq_cmplx_2d!( + f64, + g[3][0], + [ + [cmplx!(459.646893789751, 0.), cmplx!(0., 0.)], + [cmplx!(0., 0.), cmplx!(437.60340375, 0.)] + ], + 2, + epsilon = 1e-5 + ); + } + + #[test] + fn test_gamma_singlet_qed() { + const NF: u8 = 3; + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let g = gamma_singlet_qed(3, 2, &mut c, NF); + // as^0 a^0 must be trivial + assert_approx_eq_cmplx_2d!(f64, g[0][0], [[cmplx!(0., 0.); 4]; 4], 4); + + // reference from Python side + assert_approx_eq_cmplx_2d!( + f64, + g[3][0], + [ + [ + cmplx!(3.857918949669738, 0.), + cmplx!(0., 0.), + cmplx!(-290.42193908689745, 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(-3.859554320251334, 0.), + cmplx!(0., 0.), + cmplx!(290.4252052962147, 0.), + cmplx!(0., 0.) + ], + [ + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(0., 0.), + cmplx!(448.0752570151872, 0.) + ] + ], + 4, + epsilon = 1e-5 + ); + } +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs new file mode 100644 index 000000000..2cbb04b24 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs @@ -0,0 +1,132 @@ +//! |LO| |QED|. +use num::complex::Complex; +use num::Zero; + +use crate::cmplx; +use crate::constants::{ChargeCombinations, CF, ED2, EU2, NC, TR}; +use crate::harmonics::cache::Cache; + +use crate::anomalous_dimensions::unpolarized::spacelike::as1; + +/// Compute the photon-quark anomalous dimension. +/// +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. +pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_gq(c, nf) / CF +} + +/// Compute the quark-photon anomalous dimension. +/// +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. +/// However, we are adding the $N_C$ and the $2n_f$ factors from $\theta$ inside the +/// definition of $\gamma_{q \gamma}^{(0)}(N)$. +pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_qg(c, nf) / TR * (NC as f64) +} + +/// Compute the photon-photon anomalous dimension. +/// +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. +pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { + let cc = ChargeCombinations { nf }; + cmplx!( + 4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2), + 0. + ) +} + +/// Compute the non-singlet anomalous dimension. +/// +/// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea]. +pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { + as1::gamma_ns(c, nf) / CF +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = ChargeCombinations { nf }; + + let gamma_ph_q = gamma_phq(c, nf); + let gamma_q_ph = gamma_qph(c, nf); + let gamma_nonsinglet = gamma_ns(c, nf); + + [ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + Complex::::zero(), + gamma_phph(c, nf), + cc.e2avg() * gamma_ph_q, + cc.vue2m() * gamma_ph_q, + ], + [ + Complex::::zero(), + cc.e2avg() * gamma_q_ph, + cc.e2avg() * gamma_nonsinglet, + cc.vue2m() * gamma_nonsinglet, + ], + [ + Complex::::zero(), + cc.vde2m() * gamma_q_ph, + cc.vde2m() * gamma_nonsinglet, + cc.e2delta() * gamma_nonsinglet, + ], + ] +} + +/// Compute the valence anomalous dimension matrix. +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = ChargeCombinations { nf }; + let g = gamma_ns(c, nf); + [ + [cc.e2avg() * g, cc.vue2m() * g], + [cc.vde2m() * g, cc.e2delta() * g], + ] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + const NF: u8 = 5; + + #[test] + fn number_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + let me = gamma_ns(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); + } + + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let me = gamma_ns(&mut c, NF) + gamma_phq(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); + } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + assert_approx_eq_cmplx!( + f64, + EU2 * gamma_qph(&mut c, cc.nu()) + + ED2 * gamma_qph(&mut c, cc.nd()) + + gamma_phph(&mut c, nf), + cmplx!(0., 0.), + epsilon = 2e-6 + ); + } + } +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs new file mode 100644 index 000000000..1a39ef6a6 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs @@ -0,0 +1,226 @@ +//! |NLO| |QED|. +use num::complex::Complex; +use num::Zero; + +use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC}; +use crate::harmonics::cache::{Cache, K}; + +use crate::anomalous_dimensions::unpolarized::spacelike::as1aem1; + +/// Compute the photon-photon anomalous dimension. +/// +/// Implements Eq. (68) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phph(c: &mut Cache, nf: u8) -> Complex { + let cc = ChargeCombinations { nf }; + (NC as f64) + * ((cc.nu() as f64) * EU2.powi(2) + (cc.nd() as f64) * ED2.powi(2)) + * (as1aem1::gamma_gph(c, nf) / CF / CA + 4.) +} + +/// Compute the up-photon anomalous dimension. +/// +/// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_uph(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_qph(c, nf) / CF +} + +/// Compute the down-photon anomalous dimension. +/// +/// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_dph(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_qph(c, nf) / CF +} + +/// Compute the photon-quark anomalous dimension. +/// +/// Implements singlet part of Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +fn gamma_phq(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let gamma = (-16.0 * (-16.0 - 27.0 * N - 13.0 * N.powu(2) - 8.0 * N.powu(3))) + / (9.0 * (-1.0 + N) * N * (1.0 + N).powu(2)) + - 16.0 * (2.0 + 3.0 * N + 2.0 * N.powu(2) + N.powu(3)) + / (3.0 * (-1.0 + N) * N * (1.0 + N).powu(2)) + * S1; + let cc = ChargeCombinations { nf }; + let eSigma2 = (NC as f64) * cc.e2avg() * (nf as f64); + eSigma2 * gamma +} + +/// Compute the photon-up anomalous dimension. +/// +/// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf) +} + +/// Compute the photon-down anomalous dimension. +/// +/// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_phd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf) +} + +/// Compute the non-singlet anomalous dimension. +/// +/// Implements singlet part of Eq. (57) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +fn gamma_nsq(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let gamma = 2.0 * (-12.0 + 20.0 * N + 47.0 * N.powu(2) + 6.0 * N.powu(3) + 3.0 * N.powu(4)) + / (9.0 * N.powu(2) * (1.0 + N).powu(2)) + - 80.0 / 9.0 * S1 + + 16.0 / 3.0 * S2; + let cc = ChargeCombinations { nf }; + let eSigma2 = (NC as f64) * cc.e2avg() * (nf as f64); + eSigma2 * gamma +} + +/// Compute the singlet-like non-singlet up anomalous dimension. +/// +/// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_nspu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the singlet-like non-singlet down anomalous dimension. +/// +/// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_nspd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the valence-like non-singlet up anomalous dimension. +/// +/// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u. +pub fn gamma_nsmu(c: &mut Cache, nf: u8) -> Complex { + EU2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the valence-like non-singlet down anomalous dimension. +/// +/// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d. +pub fn gamma_nsmd(c: &mut Cache, nf: u8) -> Complex { + ED2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf) +} + +/// Compute the pure-singlet anomalous dimension. +/// +/// Implements Eq. (59) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk]. +pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let result = -4.0 * (2.0 + N * (5.0 + N)) * (4.0 + N * (4.0 + N * (7.0 + 5.0 * N))) + / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3) * (2.0 + N).powu(2)); + 2. * (nf as f64) * CA * result +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = ChargeCombinations { nf }; + const E2M: f64 = EU2 - ED2; + + let gamma_phu = gamma_phu(c, nf); + let gamma_phd = gamma_phd(c, nf); + let gamma_uph = gamma_uph(c, nf); + let gamma_dph = gamma_dph(c, nf); + let gamma_nspu = gamma_nspu(c, nf); + let gamma_nspd = gamma_nspd(c, nf); + let gamma_ps = gamma_ps(c, nf); + + [ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [ + Complex::::zero(), + gamma_phph(c, nf), + cc.vu() * EU2 * gamma_phu + cc.vd() * ED2 * gamma_phd, + cc.vu() * (EU2 * gamma_phu - ED2 * gamma_phd), + ], + [ + Complex::::zero(), + cc.vu() * EU2 * gamma_uph + cc.vd() * ED2 * gamma_dph, + cc.vu() * EU2 * gamma_nspu + cc.vd() * ED2 * gamma_nspd + cc.e2avg().powi(2) * gamma_ps, + cc.vu() * (EU2 * gamma_nspu - ED2 * gamma_nspd + E2M * cc.e2avg() * gamma_ps), + ], + [ + Complex::::zero(), + cc.vd() * (EU2 * gamma_uph - ED2 * gamma_dph), + cc.vd() * (EU2 * gamma_nspu - ED2 * gamma_nspd + E2M * cc.e2avg() * gamma_ps), + cc.vd() * EU2 * gamma_nspu + + cc.vu() * ED2 * gamma_nspd + + cc.vu() * cc.vd() * E2M.powi(2) * gamma_ps, + ], + ] +} + +/// Compute the valence anomalous dimension matrix. +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = ChargeCombinations { nf }; + let gamma_nsmu = gamma_nsmu(c, nf); + let gamma_nsmd = gamma_nsmd(c, nf); + [ + [ + cc.vu() * EU2 * gamma_nsmu + cc.vd() * ED2 * gamma_nsmd, + cc.vu() * (EU2 * gamma_nsmu - ED2 * gamma_nsmd), + ], + [ + cc.vd() * (EU2 * gamma_nsmu - ED2 * gamma_nsmd), + cc.vd() * EU2 * gamma_nsmu + cc.vu() * ED2 * gamma_nsmd, + ], + ] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + + #[test] + fn number_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + for nf in 2..7 { + assert_approx_eq_cmplx!(f64, gamma_nsmu(&mut c, nf), Complex::zero(), epsilon = 1e-5); + assert_approx_eq_cmplx!(f64, gamma_nsmd(&mut c, nf), Complex::zero(), epsilon = 1e-5); + } + } + + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + let ps = EU2 * gamma_ps(&mut c, cc.nu()) + ED2 * gamma_ps(&mut c, cc.nd()); + let me_u = gamma_nspu(&mut c, nf) + ps + gamma_phu(&mut c, nf); + assert_approx_eq_cmplx!(f64, me_u, Complex::zero(), epsilon = 1e-5); + let me_d = gamma_nspd(&mut c, nf) + ps + gamma_phd(&mut c, nf); + assert_approx_eq_cmplx!(f64, me_d, Complex::zero(), epsilon = 1e-5); + } + } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + assert_approx_eq_cmplx!( + f64, + EU2 * gamma_uph(&mut c, cc.nu()) + + ED2 * gamma_dph(&mut c, cc.nd()) + + gamma_phph(&mut c, nf), + Complex::zero(), + epsilon = 1e-12 + ); + } + } +} diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs new file mode 100644 index 000000000..ac4596edf --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs @@ -0,0 +1,330 @@ +//! |LO| |QED| x |LO| |QCD|. +use crate::cmplx; +use num::complex::Complex; + +use crate::constants::{ChargeCombinations, CA, CF, ED2, EU2, NC, TR, ZETA2, ZETA3}; +use crate::harmonics::cache::{recursive_harmonic_sum, Cache, K}; +use std::f64::consts::PI; + +/// Compute the photon-quark anomalous dimension. +/// +/// Implements Eq. (36) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + + #[rustfmt::skip] + let tmp_const = + 2.0 + * ( + -4.0 + - 12.0 * N + - N.powu(2) + + 28.0 * N.powu(3) + + 43.0 * N.powu(4) + + 30.0 * N.powu(5) + + 12.0 * N.powu(6) + ) / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3)); + + #[rustfmt::skip] + let tmp_S1 = -4.0 + * (10.0 + 27.0 * N + 25.0 * N.powu(2) + 13.0 * N.powu(3) + 5.0 * N.powu(4)) + / ((-1.0 + N) * N * (1.0 + N).powu(3)); + + let tmp_S12 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N)); + let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N)); + + CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) +} + +/// Compute the quark-photon anomalous dimension. +/// +/// Implements Eq. (26) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + + #[rustfmt::skip] + let tmp_const = -2.0 + * (4.0 + + 8.0 * N + + 25.0 * N.powu(2) + + 51.0 * N.powu(3) + + 36.0 * N.powu(4) + + 15.0 * N.powu(5) + + 5.0 * N.powu(6)) + / (N.powu(3) * (1.0 + N).powu(3) * (2.0 + N)); + + let tmp_S1 = 8.0 / N.powu(2); + let tmp_S12 = -4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N)); + let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N)); + + 2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2) +} + +/// Compute the gluon-photon anomalous dimension. +/// +/// Implements Eq. (27) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + CF * CA + * (8.0 * (-4.0 + N * (-4.0 + N * (-5.0 + N * (-10.0 + N + 2.0 * N.powu(2) * (2.0 + N)))))) + / (N.powu(3) * (1.0 + N).powu(3) * (-2.0 + N + N.powu(2))) +} + +/// Compute the photon-gluon anomalous dimension. +/// +/// Implements Eq. (30) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex { + TR / CF / CA * (NC as f64) * gamma_gph(c, nf) +} + +/// Compute the quark-gluon singlet anomalous dimension. +/// +/// Implements Eq. (29) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + TR / CF / CA * (NC as f64) * gamma_qph(c, nf) +} + +/// Compute the gluon-quark singlet anomalous dimension. +/// +/// Implements Eq. (35) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { + gamma_phq(c, nf) +} + +/// Compute the photon-photon singlet anomalous dimension. +/// +/// Implements Eq. (28) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex { + let cc = ChargeCombinations { nf }; + cmplx!( + 4.0 * CF * CA * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2), + 0. + ) +} + +/// Compute the gluon-gluon singlet anomalous dimension. +/// +/// Implements Eq. (31) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex { + cmplx!(4.0 * TR * (NC as f64), 0.) +} + +/// Shift for $g_3(N)$ by 2. +/// +/// This is $g_3(N+2) - g_3(N)$. +fn g3_shift(c: &mut Cache) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + (6. * (n + 1.) * (2. * n + 1.) * S1 + n * (-PI.powi(2) * (n + 1.).powu(2) - 6. * n)) + / (6. * n.powu(2) * (n + 1.).powu(3)) +} + +/// Compute the singlet-like non-singlet anomalous dimension. +/// +/// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + let S1h = c.get(K::S1h); + let S2h = c.get(K::S2h); + let S3h = c.get(K::S3h); + let S1p1h = recursive_harmonic_sum(c.get(K::S1mh), (N - 1.) / 2., 1, 1); + let S2p1h = recursive_harmonic_sum(c.get(K::S2mh), (N - 1.) / 2., 1, 2); + let S3p1h = recursive_harmonic_sum(c.get(K::S3mh), (N - 1.) / 2., 1, 3); + + let g3N = c.get(K::G3); + let g3Np2 = g3N + g3_shift(c); + + #[rustfmt::skip] + let result = 32.0 * ZETA2 * S1h - 32.0 * ZETA2 * S1p1h + + 8.0 / (N + N.powu(2)) * S2h + - 4.0 * S3h + (24.0 + 16.0 / (N + N.powu(2))) * S2 + - 32.0 * S3 - 8.0 / (N + N.powu(2)) * S2p1h + + S1 * (16.0 * (3.0 / N.powu(2) - 3.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) - 16.0 * S2h + - 32.0 * S2 + 16.0 * S2p1h ) + + (-8.0 + N * (-32.0 + N * ( -8.0 - 3.0 * N * (3.0 + N) * (3.0 + N.powu(2)) - 48.0 * (1.0 + N).powu(2) * ZETA2))) + / (N.powu(3) * (1.0 + N).powu(3)) + + 32.0 * (g3N + g3Np2) + 4.0 * S3p1h - 16.0 * ZETA3; + + CF * result +} + +/// Compute the valence-like non-singlet anomalous dimension. +/// +/// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt]. +pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + let S1h = c.get(K::S1h); + let S2h = c.get(K::S2h); + let S3h = c.get(K::S3h); + let S1p1h = recursive_harmonic_sum(c.get(K::S1mh), (N - 1.) / 2., 1, 1); + let S2p1h = recursive_harmonic_sum(c.get(K::S2mh), (N - 1.) / 2., 1, 2); + let S3p1h = recursive_harmonic_sum(c.get(K::S3mh), (N - 1.) / 2., 1, 3); + let g3N = c.get(K::G3); + let g3Np2 = g3N + g3_shift(c); + + #[rustfmt::skip] + let result = + -32.0 * ZETA2 * S1h + - 8.0 / (N + N.powu(2)) * S2h + + (24.0 + 16.0 / (N + N.powu(2))) * S2 + + 8.0 / (N + N.powu(2)) * S2p1h + + S1 + * ( + 16.0 * (-1.0 / N.powu(2) + 1.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) + + 16.0 * S2h + - 32.0 * S2 + - 16.0 * S2p1h + ) + + ( + 72.0 + + N + * ( + 96.0 + - 3.0 * N * (8.0 + 3.0 * N * (3.0 + N) * (3.0 + N.powu(2))) + + 48.0 * N * (1.0 + N).powu(2) * ZETA2 + ) + ) + / (3.0 * N.powu(3) * (1.0 + N).powu(3)) + - 32.0 * (g3N + g3Np2) + + 32.0 * ZETA2 * S1p1h + + 4.0 * S3h + - 32.0 * S3 + - 4.0 * S3p1h + - 16.0 * ZETA3; + + CF * result +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 4]; 4] { + let cc = ChargeCombinations { nf }; + let e2_tot = nf as f64 * cc.e2avg(); + let gq = gamma_gq(c, nf); + let phq = gamma_phq(c, nf); + let qg = gamma_qg(c, nf); + let qph = gamma_qph(c, nf); + let nsp = gamma_nsp(c, nf); + [ + [ + e2_tot * gamma_gg(c, nf), + e2_tot * gamma_gph(c, nf), + cc.e2avg() * gq, + cc.vue2m() * gq, + ], + [ + e2_tot * gamma_phg(c, nf), + gamma_phph(c, nf), + cc.e2avg() * phq, + cc.vue2m() * phq, + ], + [ + cc.e2avg() * qg, + cc.e2avg() * qph, + cc.e2avg() * nsp, + cc.vue2m() * nsp, + ], + [ + cc.vde2m() * qg, + cc.vde2m() * qph, + cc.vde2m() * nsp, + cc.e2delta() * nsp, + ], + ] +} + +/// Compute the valence anomalous dimension matrix. +pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let cc = ChargeCombinations { nf }; + let g = gamma_nsm(c, nf); + [ + [cc.e2avg() * g, cc.vue2m() * g], + [cc.vde2m() * g, cc.e2delta() * g], + ] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::Zero; + const NF: u8 = 5; + + #[test] + fn number_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + let me = gamma_nsm(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-4); + } + + #[test] + fn gluon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + assert_approx_eq_cmplx!( + f64, + EU2 * gamma_qg(&mut c, cc.nu()) + + ED2 * gamma_qg(&mut c, cc.nd()) + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_phg(&mut c, nf) + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_gg(&mut c, nf), + cmplx!(0., 0.), + epsilon = 1e-14 + ); + } + } + + #[test] + fn photon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + + for nf in 2..7 { + let cc = ChargeCombinations { nf }; + assert_approx_eq_cmplx!( + f64, + EU2 * gamma_qph(&mut c, cc.nu()) + + ED2 * gamma_qph(&mut c, cc.nd()) + + gamma_phph(&mut c, nf) + + (cc.nu() as f64 * EU2 + cc.nd() as f64 * ED2) * gamma_gph(&mut c, nf), + cmplx!(0., 0.), + epsilon = 1e-14 + ); + } + } + + #[test] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let me = gamma_nsp(&mut c, NF) + gamma_gq(&mut c, NF) + gamma_phq(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-4); + } + + #[test] + fn test_g3_shift() { + for N in [cmplx!(1.23, 0.), cmplx!(5., 0.), cmplx!(2., 1.)] { + let mut c = Cache::new(N); + let mut c2 = Cache::new(N + 2.); + assert_approx_eq_cmplx!( + f64, + g3_shift(&mut c), + c2.get(K::G3) - c.get(K::G3), + epsilon = 1e-6 + ); + } + } +} diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index 3e8a1eb8f..89ca13ccf 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2024-09-10T11:22:16.645053). +//! List of References (autogenerated on 2025-01-13T15:14:10.528903). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -107,3 +107,51 @@ pub fn Buza1996wv() {} /// /// DOI: [10.1016/j.physletb.2009.01.057](https:dx.doi.org/10.1016/j.physletb.2009.01.057) pub fn Bierenbaum2009zt() {} + +#[allow(non_snake_case)] +/// Next-to-leading order radiative parton model analysis of polarized deep inelastic lepton - nucleon scattering +/// +/// Gluck, M. and Reya, E. and Stratmann, M. and Vogelsang, W. +/// +/// Published in: Phys. Rev. D 53 (1996), 4775--4786 +/// +/// e-Print: [hep-ph/9508347](https://arxiv.org/abs/hep-ph/9508347) +/// +/// DOI: [10.1103/PhysRevD.53.4775](https:dx.doi.org/10.1103/PhysRevD.53.4775) +pub fn Gluck1995yr() {} + +#[allow(non_snake_case)] +/// Parton distribution functions with QED corrections +/// +/// Carrazza, Stefano +/// +/// Published as PhD thesis at Milan U. (2015) +/// +/// e-Print: [1509.00209](https://arxiv.org/abs/1509.00209) +/// +/// DOI: [10.13130/carrazza-stefano_phd2015-07-06](https:dx.doi.org/10.13130/carrazza-stefano_phd2015-07-06) +pub fn Carrazza2015dea() {} + +#[allow(non_snake_case)] +/// QED corrections to the Altarelli-Parisi splitting functions +/// +/// de Florian, Daniel and Sborlini, Germán F. R. and Rodrigo, Germán +/// +/// Published in: Eur. Phys. J. C 76 (2016), 282 +/// +/// e-Print: [1512.00612](https://arxiv.org/abs/1512.00612) +/// +/// DOI: [10.1140/epjc/s10052-016-4131-8](https:dx.doi.org/10.1140/epjc/s10052-016-4131-8) +pub fn deFlorian2015ujt() {} + +#[allow(non_snake_case)] +/// Two-loop QED corrections to the Altarelli-Parisi splitting functions +/// +/// de Florian, Daniel and Sborlini, Germán F. R. and Rodrigo, Germán +/// +/// Published in: JHEP 10 (2016), 056 +/// +/// e-Print: [1606.02887](https://arxiv.org/abs/1606.02887) +/// +/// DOI: [10.1007/JHEP10(2016)056](https:dx.doi.org/10.1007/JHEP10(2016)056) +pub fn deFlorian2016gvk() {} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index 8cb90eb1f..98023e016 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -20,6 +20,16 @@ pub const CA: f64 = NC as f64; /// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$. pub const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64); +/// Up quark electric charge square. +/// +/// Defaults to $e_u^2 = 4./9$. +pub const EU2: f64 = 4. / 9.; + +/// Down quark electric charge square. +/// +/// Defaults to $e_d^2 = 1./9$. +pub const ED2: f64 = 1. / 9.; + /// Riemann zeta function at z = 2. /// /// $\zeta(2) = \pi^2 / 6$. @@ -41,3 +51,62 @@ pub const PID_NSM: u16 = 10201; /// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; + +/// singlet-like non-singlet up-sector |PID|. +pub const PID_NSP_U: u16 = 10102; + +/// singlet-like non-singlet down-sector |PID|. +pub const PID_NSP_D: u16 = 10103; + +/// valence-like non-singlet up-sector |PID|. +pub const PID_NSM_U: u16 = 10202; + +/// valence-like non-singlet down-sector |PID|. +pub const PID_NSM_D: u16 = 10203; + +/// |QED| electric charge combinations. +pub struct ChargeCombinations { + pub nf: u8, +} + +impl ChargeCombinations { + /// Number of up-like flavors. + pub fn nu(&self) -> u8 { + self.nf / 2 + } + + /// Number of down-like flavors. + pub fn nd(&self) -> u8 { + self.nf - self.nu() + } + + /// Charge average. + pub fn e2avg(&self) -> f64 { + (self.nu() as f64 * EU2 + self.nd() as f64 * ED2) / (self.nf as f64) + } + + /// Relative up contribution. + pub fn vu(&self) -> f64 { + self.nu() as f64 / (self.nf as f64) + } + + /// Relative down contribution. + pub fn vd(&self) -> f64 { + self.nd() as f64 / (self.nf as f64) + } + + /// Relative up contribution to charge difference. + pub fn vue2m(&self) -> f64 { + self.vu() * (EU2 - ED2) + } + + /// Relative down contribution to charge difference. + pub fn vde2m(&self) -> f64 { + self.vd() * (EU2 - ED2) + } + + /// Asymmetric charge combination. + pub fn e2delta(&self) -> f64 { + self.vde2m() - self.vue2m() + self.e2avg() + } +} diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs index b5c812165..a9a8a818d 100644 --- a/crates/ekore/src/util.rs +++ b/crates/ekore/src/util.rs @@ -22,3 +22,27 @@ macro_rules! assert_approx_eq_cmplx { float_cmp::assert_approx_eq!($size, $ref.im, $target.im $(, $set = $val)*); }; } + +/// Shorthand complex list comparators. +#[cfg(test)] +#[macro_export] +macro_rules! assert_approx_eq_cmplx_1d { + ($size:ty, $ref:expr, $target:expr, $d:expr $(, $set:ident = $val:expr)*) => { + for j in 0..$d { + assert_approx_eq_cmplx!($size, $ref[j], $target[j] $(, $set = $val)*); + } + } +} + +/// Shorthand complex matrix comparators. +#[cfg(test)] +#[macro_export] +macro_rules! assert_approx_eq_cmplx_2d { + ($size:ty, $ref:expr, $target:expr, $d:expr $(, $set:ident = $val:expr)*) => { + for j in 0..$d { + for k in 0..$d { + assert_approx_eq_cmplx!($size, $ref[j][k], $target[j][k] $(, $set = $val)*); + } + } + } +} diff --git a/crates/make_bib.py b/crates/make_bib.py index e758c41ad..ef1f06472 100644 --- a/crates/make_bib.py +++ b/crates/make_bib.py @@ -25,9 +25,15 @@ PHD = """Published as PhD thesis at {school} ({year})""" -def clean_nl(t: str) -> str: - """Shrink whitespace to one.""" - return re.sub("\n\\s+", " ", t) +def sanitize(t: str) -> str: + """Reformat for HTML.""" + # Shrink whitespace to one. + t = re.sub("\n\\s+", " ", t) + # Recover some TeX + t = t.replace(r"\textendash{}", "-").replace(r"\'a", "á") + # Remove curly brackets + t = re.sub(r"^\{(.+)\}$", r"\1", t) + return t # collect output @@ -39,8 +45,8 @@ def clean_nl(t: str) -> str: bib_database = bibtexparser.parse_string(bib_str) # iterate on all elements for el in bib_database.entries: - title = re.sub(r"^\{(.+)\}$", r"\1", clean_nl(el.fields_dict["title"].value)) - author = el.fields_dict["author"].value + title = sanitize(el.fields_dict["title"].value) + author = sanitize(el.fields_dict["author"].value) if el.entry_type == "phdthesis": publication = PHD.format( school=el.fields_dict["school"].value, diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index a4951e687..576c34ae6 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,8 +1,8 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index bd1b19d6..f543f7bc 100644 +index bd1b19d6..24f99075 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py -@@ -3,16 +3,16 @@ r"""Contains the central operator classes. +@@ -3,120 +3,28 @@ r"""Contains the central operator classes. See :doc:`Operator overview `. """ @@ -17,15 +17,28 @@ index bd1b19d6..f543f7bc 100644 import numba as nb import numpy as np -from scipy import integrate +- +-import ekore.anomalous_dimensions.polarized.space_like as ad_ps +-import ekore.anomalous_dimensions.unpolarized.space_like as ad_us +-import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut +from scipy import LowLevelCallable, integrate - import ekore.anomalous_dimensions.polarized.space_like as ad_ps - import ekore.anomalous_dimensions.unpolarized.space_like as ad_us -@@ -32,91 +32,10 @@ from ..matchings import Segment, lepton_number + from .. import basis_rotation as br + from .. import interpolation, mellin + from .. import scale_variations as sv + from ..io.types import EvolutionMethod, OperatorLabel + from ..kernels import ev_method +-from ..kernels import non_singlet as ns +-from ..kernels import non_singlet_qed as qed_ns +-from ..kernels import singlet as s +-from ..kernels import singlet_qed as qed_s +-from ..kernels import valence_qed as qed_v +-from ..matchings import Segment, lepton_number ++from ..matchings import Segment from ..member import OpMember - from ..scale_variations import expanded as sv_expanded - from ..scale_variations import exponentiated as sv_exponentiated -+from .quad_ker import cb_quad_ker_qcd +-from ..scale_variations import expanded as sv_expanded +-from ..scale_variations import exponentiated as sv_exponentiated ++from .quad_ker import cb_quad_ker_qcd, cb_quad_ker_qed logger = logging.getLogger(__name__) @@ -114,7 +127,7 @@ index bd1b19d6..f543f7bc 100644 spec = [ ("is_singlet", nb.boolean), ("is_QEDsinglet", nb.boolean), -@@ -188,421 +107,6 @@ class QuadKerBase: +@@ -188,422 +96,6 @@ class QuadKerBase: return self.path.prefactor * pj * self.path.jac @@ -533,10 +546,11 @@ index bd1b19d6..f543f7bc 100644 - ) - return ker - - +- OpMembers = Dict[OperatorLabel, OpMember] """Map of all operators.""" -@@ -792,49 +296,6 @@ class Operator(sv.ScaleVariationModeMixin): + +@@ -792,49 +284,6 @@ class Operator(sv.ScaleVariationModeMixin): """Return the evolution method.""" return ev_method(EvolutionMethod(self.config["method"])) @@ -586,7 +600,7 @@ index bd1b19d6..f543f7bc 100644 def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( -@@ -857,10 +318,14 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -857,10 +306,29 @@ class Operator(sv.ScaleVariationModeMixin): else: self.op_members[n] = zero.copy() @@ -596,16 +610,31 @@ index bd1b19d6..f543f7bc 100644 - ): + def update_cfg(self, cfg): + """Adjust integration config.""" ++ # prepare couplings for c + cfg.as1 = self.as_list[1] + cfg.as0 = self.as_list[0] -+ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) ++ as_list_len = self.as_list.shape[0] ++ for j, a in enumerate(self.as_list.tolist()): ++ cfg.as_list[j] = a ++ cfg.as_list_len = as_list_len ++ a_half_x = self.a_half_list.shape[0] ++ a_half_y = self.a_half_list.shape[1] ++ for j, a in enumerate(self.a_half_list.flatten().tolist()): ++ cfg.a_half[j] = a ++ cfg.a_half_x = a_half_x ++ cfg.a_half_y = a_half_y ++ ++ if self.order[1] == 0: ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qcd.address) ++ else: ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_qed.address) + cfg.method_num = self.ev_method + + def run_op_integration(self, log_grid): """Run the integration for each grid point. Parameters -@@ -875,18 +339,53 @@ class Operator(sv.ScaleVariationModeMixin): +@@ -875,18 +343,62 @@ class Operator(sv.ScaleVariationModeMixin): """ column = [] k, logx = log_grid @@ -613,8 +642,9 @@ index bd1b19d6..f543f7bc 100644 + labels = self.labels start_time = time.perf_counter() + # start preparing C arguments -+ cfg = ekors.lib.empty_qcd_args() ++ cfg = ekors.lib.empty_args() + cfg.order_qcd = self.order[0] ++ cfg.order_qed = self.order[1] + cfg.is_polarized = self.config["polarized"] + cfg.is_time_like = self.config["time_like"] + cfg.nf = self.nf @@ -625,6 +655,18 @@ index bd1b19d6..f543f7bc 100644 + cfg.ev_op_max_order_qcd = self.config["ev_op_max_order"][0] + cfg.sv_mode_num = self.sv_mode + cfg.is_threshold = self.is_threshold ++ cfg.mu2_from = self.q2_from ++ cfg.mu2_to = self.q2_to ++ cfg.alphaem_running = self.alphaem_running ++ a_half_ffi = ekors.ffi.new("double[]", 2 * self.config["ev_op_iterations"]) ++ cfg.a_half = a_half_ffi ++ as_list_ffi = ekors.ffi.new("double[]", self.config["ev_op_iterations"] + 1) ++ cfg.as_list = as_list_ffi ++ max_areas_shape = self.int_disp.max_areas_shape() ++ max_area_len = max_areas_shape[0] * max_areas_shape[1] ++ areas_ffi = ekors.ffi.new("double[]", max_area_len) ++ cfg.areas = areas_ffi ++ + self.update_cfg(cfg) + # iterate basis functions @@ -638,12 +680,8 @@ index bd1b19d6..f543f7bc 100644 temp_dict = {} + # prepare areas for C + curareas = bf.areas_representation -+ areas_len = curareas.shape[0] * curareas.shape[1] -+ # force the variable in scope -+ areas_ffi = ekors.ffi.new( -+ f"double[{areas_len}]", curareas.flatten().tolist() -+ ) -+ cfg.areas = areas_ffi ++ for j, x in enumerate(curareas.flatten().tolist()): ++ cfg.areas[j] = x + cfg.areas_x = curareas.shape[0] + cfg.areas_y = curareas.shape[1] # iterate sectors @@ -653,7 +691,7 @@ index bd1b19d6..f543f7bc 100644 + cfg.mode1 = label[1] + # construct the low level object + func = LowLevelCallable( -+ ekors.lib.rust_quad_ker_qcd, ekors.ffi.addressof(cfg) ++ ekors.lib.rust_quad_ker, ekors.ffi.addressof(cfg) + ) res = integrate.quad( - self.quad_ker( diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 087f6bbb7..ba111e3bc 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -9,7 +9,11 @@ from .. import scale_variations as sv from ..io.types import InversionMethod from ..kernels import non_singlet as ns +from ..kernels import non_singlet_qed as qed_ns from ..kernels import singlet as s +from ..kernels import singlet_qed as qed_s +from ..kernels import valence_qed as qed_v +from ..matchings import lepton_number from ..scale_variations import expanded as sv_expanded from ..scale_variations import exponentiated as sv_exponentiated @@ -47,6 +51,7 @@ def select_singlet_element(ker, mode0, mode1): nb.types.double, # re_jac nb.types.double, # im_jac nb.types.uintc, # order_qcd + nb.types.uintc, # order_qed nb.types.bool_, # is_singlet nb.types.uintc, # mode0 nb.types.uintc, # mode1 @@ -56,7 +61,6 @@ def select_singlet_element(ker, mode0, mode1): nb.types.CPointer(nb.types.double), # areas_raw nb.types.uintc, # areas_x nb.types.uintc, # areas_y - nb.types.double, # L nb.types.uintc, # method_num nb.types.double, # as1 nb.types.double, # as0 @@ -65,6 +69,14 @@ def select_singlet_element(ker, mode0, mode1): nb.types.uintc, # sv_mode_num nb.types.bool_, # is_threshold nb.types.double, # Lsv + nb.types.CPointer(nb.types.double), # as_list + nb.types.uintc, # as_list_len + nb.types.double, # mu2_from + nb.types.double, # mu2_to + nb.types.CPointer(nb.types.double), # a_half + nb.types.uintc, # a_half_x + nb.types.uintc, # a_half_y + nb.types.bool_, # alphaem_running ) @@ -81,6 +93,7 @@ def cb_quad_ker_qcd( re_jac, im_jac, order_qcd, + _order_qed, is_singlet, mode0, mode1, @@ -90,7 +103,6 @@ def cb_quad_ker_qcd( areas_raw, areas_x, areas_y, - _L, ev_method, as1, as0, @@ -99,6 +111,14 @@ def cb_quad_ker_qcd( sv_mode, is_threshold, Lsv, + _as_list, + _as_list_len, + _mu2_from, + _mu2_to, + _a_half, + _a_half_x, + _a_half_y, + _alphaem_running, ): """C Callback inside integration kernel.""" # recover complex variables @@ -225,6 +245,7 @@ def cb_quad_ker_ome( re_jac, im_jac, order_qcd, + _order_qed, is_singlet, mode0, mode1, @@ -234,7 +255,6 @@ def cb_quad_ker_ome( areas_raw, areas_x, areas_y, - L, backward_method, as1, _as0, @@ -243,6 +263,14 @@ def cb_quad_ker_ome( sv_mode, _is_threshold, Lsv, + _as_list, + _as_list_len, + _mu2_from, + _mu2_to, + _a_half, + _a_half_x, + _a_half_y, + _alphaem_running, ): """C Callback inside integration kernel.""" # recover complex variables @@ -280,210 +308,224 @@ def cb_quad_ker_ome( return np.real(res) -# from ..kernels import singlet_qed as qed_s -# from ..kernels import non_singlet_qed as qed_ns -# from ..kernels import valence_qed as qed_v - -# @nb.njit(cache=True) -# def select_QEDsinglet_element(ker, mode0, mode1): -# """Select element of the QEDsinglet matrix. - -# Parameters -# ---------- -# ker : numpy.ndarray -# QEDsinglet integration kernel -# mode0 : int -# id for first sector element -# mode1 : int -# id for second sector element -# Returns -# ------- -# ker : complex -# QEDsinglet integration kernel element -# """ -# if mode0 == 21: -# index1 = 0 -# elif mode0 == 22: -# index1 = 1 -# elif mode0 == 100: -# index1 = 2 -# else: -# index1 = 3 -# if mode1 == 21: -# index2 = 0 -# elif mode1 == 22: -# index2 = 1 -# elif mode1 == 100: -# index2 = 2 -# else: -# index2 = 3 -# return ker[index1, index2] - - -# @nb.njit(cache=True) -# def select_QEDvalence_element(ker, mode0, mode1): -# """ -# Select element of the QEDvalence matrix. - -# Parameters -# ---------- -# ker : numpy.ndarray -# QEDvalence integration kernel -# mode0 : int -# id for first sector element -# mode1 : int -# id for second sector element -# Returns -# ------- -# ker : complex -# QEDvalence integration kernel element -# """ -# index1 = 0 if mode0 == 10200 else 1 -# index2 = 0 if mode1 == 10200 else 1 -# return ker[index1, index2] - - -# @nb.njit(cache=True) -# def quad_ker_qed( -# ker_base, -# order, -# mode0, -# mode1, -# method, -# as_list, -# mu2_from, -# mu2_to, -# a_half, -# alphaem_running, -# nf, -# L, -# ev_op_iterations, -# ev_op_max_order, -# sv_mode, -# is_threshold, -# ): -# """Raw evolution kernel inside quad. - -# Parameters -# ---------- -# ker_base : QuadKerBase -# quad argument -# order : int -# perturbation order -# mode0: int -# pid for first sector element -# mode1 : int -# pid for second sector element -# method : str -# method -# as1 : float -# target coupling value -# as0 : float -# initial coupling value -# mu2_from : float -# initial value of mu2 -# mu2_from : float -# final value of mu2 -# aem_list : list -# list of electromagnetic coupling values -# alphaem_running : bool -# whether alphaem is running or not -# nf : int -# number of active flavors -# L : float -# logarithm of the squared ratio of factorization and renormalization scale -# ev_op_iterations : int -# number of evolution steps -# ev_op_max_order : int -# perturbative expansion order of U -# sv_mode: int, `enum.IntEnum` -# scale variation mode, see `eko.scale_variations.Modes` -# is_threshold : boolean -# is this an itermediate threshold operator? - -# Returns -# ------- -# float -# evaluated integration kernel -# """ -# # compute the actual evolution kernel for QEDxQCD -# if ker_base.is_QEDsinglet: -# gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf) -# # scale var exponentiated is directly applied on gamma -# if sv_mode == sv.Modes.exponentiated: -# gamma_s = sv.exponentiated.gamma_variation_qed( -# gamma_s, order, nf, L, alphaem_running -# ) -# ker = qed_s.dispatcher( -# order, -# method, -# gamma_s, -# as_list, -# a_half, -# nf, -# ev_op_iterations, -# ev_op_max_order, -# ) -# # scale var expanded is applied on the kernel -# # TODO : in this way a_half[-1][1] is the aem value computed in -# # the middle point of the last step. Instead we want aem computed in mu2_final. -# # However the distance between the two is very small and affects only the running aem -# if sv_mode == sv.Modes.expanded and not is_threshold: -# ker = np.ascontiguousarray( -# sv.expanded.singlet_variation_qed( -# gamma_s, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L -# ) -# ) @ np.ascontiguousarray(ker) -# ker = select_QEDsinglet_element(ker, mode0, mode1) -# elif ker_base.is_QEDvalence: -# gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) -# # scale var exponentiated is directly applied on gamma -# if sv_mode == sv.Modes.exponentiated: -# gamma_v = sv.exponentiated.gamma_variation_qed( -# gamma_v, order, nf, L, alphaem_running -# ) -# ker = qed_v.dispatcher( -# order, -# method, -# gamma_v, -# as_list, -# a_half, -# nf, -# ev_op_iterations, -# ev_op_max_order, -# ) -# # scale var expanded is applied on the kernel -# if sv_mode == sv.Modes.expanded and not is_threshold: -# ker = np.ascontiguousarray( -# sv.expanded.valence_variation_qed( -# gamma_v, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L -# ) -# ) @ np.ascontiguousarray(ker) -# ker = select_QEDvalence_element(ker, mode0, mode1) -# else: -# gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) -# # scale var exponentiated is directly applied on gamma -# if sv_mode == sv.Modes.exponentiated: -# gamma_ns = sv.exponentiated.gamma_variation_qed( -# gamma_ns, order, nf, L, alphaem_running -# ) -# ker = qed_ns.dispatcher( -# order, -# method, -# gamma_ns, -# as_list, -# a_half[:, 1], -# alphaem_running, -# nf, -# ev_op_iterations, -# mu2_from, -# mu2_to, -# ) -# if sv_mode == sv.Modes.expanded and not is_threshold: -# ker = ( -# sv.expanded.non_singlet_variation_qed( -# gamma_ns, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L -# ) -# * ker -# ) -# return ker +@nb.njit(cache=True) +def select_QEDsinglet_element(ker, mode0, mode1): + """Select element of the QEDsinglet matrix. + + Parameters + ---------- + ker : numpy.ndarray + QEDsinglet integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + Returns + ------- + ker : complex + QEDsinglet integration kernel element + """ + if mode0 == 21: + index1 = 0 + elif mode0 == 22: + index1 = 1 + elif mode0 == 100: + index1 = 2 + else: + index1 = 3 + if mode1 == 21: + index2 = 0 + elif mode1 == 22: + index2 = 1 + elif mode1 == 100: + index2 = 2 + else: + index2 = 3 + return ker[index1, index2] + + +@nb.njit(cache=True) +def select_QEDvalence_element(ker, mode0, mode1): + """Select element of the QEDvalence matrix. + + Parameters + ---------- + ker : numpy.ndarray + QEDvalence integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + Returns + ------- + ker : complex + QEDvalence integration kernel element + """ + index1 = 0 if mode0 == 10200 else 1 + index2 = 0 if mode1 == 10200 else 1 + return ker[index1, index2] + + +@nb.cfunc( + CB_SIGNATURE, + cache=True, + nopython=True, +) +def cb_quad_ker_qed( + re_gamma_raw, + im_gamma_raw, + re_n, + im_n, + re_jac, + im_jac, + order_qcd, + order_qed, + is_singlet, + mode0, + mode1, + nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, + ev_method, + _as1, + _as0, + ev_op_iterations, + ev_op_max_order_qcd, + sv_mode, + is_threshold, + Lsv, + as_list_raw, + as_list_len, + mu2_from, + mu2_to, + a_half_raw, + a_half_x, + a_half_y, + alphaem_running, +): + """C Callback inside integration kernel.""" + # recover complex variables + n = re_n + im_n * 1j + jac = re_jac + im_jac * 1j + # compute basis functions + areas = nb.carray(areas_raw, (areas_x, areas_y)) + pj = interpolation.evaluate_grid(n, is_log, logx, areas) + order = (order_qcd, order_qed) + ev_op_max_order = (ev_op_max_order_qcd, order_qed) + is_valence = (mode0 == 10200) or (mode0 == 10204) + + as_list = nb.carray(as_list_raw, as_list_len) + a_half = nb.carray(a_half_raw, (a_half_x, a_half_y)) + + if is_singlet: + # reconstruct singlet matrices + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1, 4, 4)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1, 4, 4)) + gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + + # scale var exponentiated is directly applied on gamma + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation_qed( + gamma_singlet, order, nf, lepton_number(mu2_to), Lsv, alphaem_running + ) + + ker = qed_s.dispatcher( + order, + ev_method, + gamma_singlet, + as_list, + a_half, + nf, + ev_op_iterations, + ev_op_max_order, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.singlet_variation_qed( + gamma_singlet, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + Lsv, + ) + ) @ np.ascontiguousarray(ker) + ker = select_QEDsinglet_element(ker, mode0, mode1) + + elif is_valence: + # reconstruct valence matrices + re_gamma_valence = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1, 2, 2)) + im_gamma_valence = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1, 2, 2)) + gamma_valence = re_gamma_valence + im_gamma_valence * 1j + + if sv_mode == sv.Modes.exponentiated: + gamma_valence = sv.exponentiated.gamma_variation_qed( + gamma_valence, order, nf, lepton_number(mu2_to), Lsv, alphaem_running + ) + ker = qed_v.dispatcher( + order, + ev_method, + gamma_valence, + as_list, + a_half, + nf, + ev_op_iterations, + ev_op_max_order, + ) + # scale var expanded is applied on the kernel + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.valence_variation_qed( + gamma_valence, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + Lsv, + ) + ) @ np.ascontiguousarray(ker) + ker = select_QEDvalence_element(ker, mode0, mode1) + + else: + # construct non-singlet matrices + re_gamma_ns = nb.carray(re_gamma_raw, (order_qcd + 1, order_qed + 1)) + im_gamma_ns = nb.carray(im_gamma_raw, (order_qcd + 1, order_qed + 1)) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv_exponentiated.gamma_variation_qed( + gamma_ns, order, nf, lepton_number(mu2_to), Lsv, alphaem_running + ) + # construct eko + ker = qed_ns.dispatcher( + order, + ev_method, + gamma_ns, + as_list, + a_half[:, 1], + alphaem_running, + nf, + ev_op_iterations, + mu2_from, + mu2_to, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = ( + sv_expanded.non_singlet_variation_qed( + gamma_ns, + as_list[-1], + a_half[-1][1], + alphaem_running, + order, + nf, + Lsv, + ) + * ker + ) + # recombine everything + res = ker * pj * jac + return np.real(res) diff --git a/src/eko/interpolation.py b/src/eko/interpolation.py index 2948f028d..526bd5c76 100644 --- a/src/eko/interpolation.py +++ b/src/eko/interpolation.py @@ -569,6 +569,11 @@ def __iter__(self): # return iter(self.basis) yield from self.basis + def max_areas_shape(self) -> tuple[int, int]: + """Maximum dimensions of the polynomial areas.""" + # 2. dim: xmin, xmax, 1+degree coefficients + return (max(len(bf.areas) for bf in self), 2 + 1 + self.polynomial_degree) + def __getitem__(self, item): return self.basis[item] diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py index 1d25d5845..0e8352dda 100644 --- a/src/eko/kernels/non_singlet.py +++ b/src/eko/kernels/non_singlet.py @@ -375,7 +375,7 @@ def dispatcher(order, method, gamma_ns, a1, a0, nf, ev_op_iterations): # pylint return eko_ordered_truncated( gamma_ns, a1, a0, betalist, order, ev_op_iterations ) - if method is EvoMethods.TRUNCATED: + if method == EvoMethods.TRUNCATED: return eko_truncated(gamma_ns, a1, a0, betalist, order, ev_op_iterations) # NLO diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index 96f02a972..2ce7bce21 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -613,7 +613,7 @@ def dispatcher( # pylint: disable=too-many-return-statements # Common method for NLO and NNLO if method in [EvoMethods.ITERATE_EXACT, EvoMethods.ITERATE_EXPANDED]: return eko_iterate(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) - if method is EvoMethods.PERTURBATIVE_EXACT: + if method == EvoMethods.PERTURBATIVE_EXACT: return eko_perturbative( gamma_singlet, a1, @@ -624,7 +624,7 @@ def dispatcher( # pylint: disable=too-many-return-statements ev_op_max_order, True, ) - if method is EvoMethods.PERTURBATIVE_EXPANDED: + if method == EvoMethods.PERTURBATIVE_EXPANDED: return eko_perturbative( gamma_singlet, a1, @@ -638,13 +638,13 @@ def dispatcher( # pylint: disable=too-many-return-statements if method in [EvoMethods.TRUNCATED, EvoMethods.ORDERED_TRUNCATED]: return eko_truncated(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) # These methods are scattered for nlo and nnlo - if method is EvoMethods.DECOMPOSE_EXACT: + if method == EvoMethods.DECOMPOSE_EXACT: if order[0] == 2: return nlo_decompose_exact(gamma_singlet, a1, a0, betalist) if order[0] == 3: return nnlo_decompose_exact(gamma_singlet, a1, a0, betalist) return n3lo_decompose_exact(gamma_singlet, a1, a0, nf) - if method is EvoMethods.DECOMPOSE_EXPANDED: + if method == EvoMethods.DECOMPOSE_EXPANDED: if order[0] == 2: return nlo_decompose_expanded(gamma_singlet, a1, a0, betalist) if order[0] == 3: diff --git a/src/eko/kernels/singlet_qed.py b/src/eko/kernels/singlet_qed.py index d63a2f1ac..0a5716df8 100644 --- a/src/eko/kernels/singlet_qed.py +++ b/src/eko/kernels/singlet_qed.py @@ -97,7 +97,7 @@ def dispatcher( e_s : numpy.ndarray singlet EKO """ - if method is EvoMethods.ITERATE_EXACT: + if method == EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_singlet, as_list, a_half, nf, order, ev_op_iterations, 4 ) diff --git a/src/eko/kernels/valence_qed.py b/src/eko/kernels/valence_qed.py index 8b83e1917..b5da158b4 100644 --- a/src/eko/kernels/valence_qed.py +++ b/src/eko/kernels/valence_qed.py @@ -45,7 +45,7 @@ def dispatcher( e_v : numpy.ndarray singlet EKO """ - if method is EvoMethods.ITERATE_EXACT: + if method == EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_valence, as_list, a_half, nf, order, ev_op_iterations, 2 ) diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index e2ea83091..0451d136c 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -147,6 +147,20 @@ def run_external(self, theory, ocard, pdf): DGLAP result """ # pylint:disable=import-error,import-outside-toplevel + if self.external.lower() == "void": + xgrid = ocard["interpolation_xgrid"] + mugrid = ocard["mugrid"] + labels = br.flavor_basis_pids + if self.rotate_to_evolution_basis: + labels = br.unified_evol_basis if theory["QED"] > 0 else br.evol_basis + + void = { + "target_xgrid": xgrid, + "values": { + mu**2: {pid: [0.0] * len(xgrid) for pid in labels} for mu in mugrid + }, + } + return void if self.external.lower() == "lha": from .external import LHA_utils @@ -211,7 +225,7 @@ def log(self, theory, _, pdf, me, ext): # LHA NNLO VFNS needs a special treatment # Valence contains only u and d rotate_to_evolution = None - labels = None + labels = br.flavor_basis_pids qed = theory["QED"] > 0 if self.rotate_to_evolution_basis: if not qed: diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py index db0f01c38..54be58b16 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py @@ -254,7 +254,7 @@ def gamma_nsmu(N, nf, cache): nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) tmp = ( - 2 + 2.0 * (-12.0 + 20.0 * N + 47.0 * N**2 + 6.0 * N**3 + 3.0 * N**4) / (9.0 * N**2 * (1.0 + N) ** 2) - 80.0 / 9.0 * S1 diff --git a/tests/eko/test_interpolation.py b/tests/eko/test_interpolation.py index b3d032415..9bb5d3de7 100644 --- a/tests/eko/test_interpolation.py +++ b/tests/eko/test_interpolation.py @@ -133,7 +133,7 @@ def test_get_interpolation(self): inter_x = interpolation.InterpolatorDispatcher(xg, 1, False) i = inter_x.get_interpolation(xg.raw) np.testing.assert_array_almost_equal(i, np.eye(len(xg))) - # .75 is exactly inbetween + # .75 is exactly in between i = inter_x.get_interpolation([0.75]) np.testing.assert_array_almost_equal(i, [[0.5, 0.5]]) @@ -163,6 +163,23 @@ def test_math(self): inter_N = interpolation.InterpolatorDispatcher(xgrid, poly_deg, mode_N=True) check_correspondence_interpolators(inter_x, inter_N) + def test_max_areas_shape(self): + for ng in [3, 6, 7, 8, 10, 11, 14]: + for pd in range(1, ng // 2): + xgrid = interpolation.XGrid(np.linspace(0.1, 1, ng), False) + interp = interpolation.InterpolatorDispatcher(xgrid, pd) + max_shape = interp.max_areas_shape() + # all areas must respect the maximum + for bf in interp: + s = bf.areas_to_const().shape + assert s[0] <= max_shape[0], f"{ng=},{pd=},k={bf.poly_number}" + assert s[1] <= max_shape[1], f"{ng=},{pd=},k={bf.poly_number}" + # the end must be the maximum, because we are leaning upwards + end = ng - pd - 1 + ms = interp[end].areas_to_const().shape + assert ms[0] == max_shape[0], f"{ng=},{pd=},{end=}" + assert ms[1] == max_shape[1], f"{ng=},{pd=},{end=}" + class TestBasisFunction: def test_init(self):