diff --git a/crates/dekoder/src/eko.rs b/crates/dekoder/src/eko.rs index 1df9d7811..50112e352 100644 --- a/crates/dekoder/src/eko.rs +++ b/crates/dekoder/src/eko.rs @@ -92,7 +92,7 @@ impl EKO { pub fn write(&self, dst: PathBuf) -> Result<()> { self.assert_working_dir()?; // create writer - let dst_file = File::create(&dst)?; + let dst_file = File::create(dst)?; let dst_file = BufWriter::with_capacity(TAR_WRITER_CAPACITY, dst_file); let mut ar = tar::Builder::new(dst_file); // do it! @@ -101,7 +101,7 @@ impl EKO { /// Extract tar file from `src` to `dst`. pub fn extract(src: PathBuf, dst: PathBuf) -> Result { - let mut ar = tar::Archive::new(File::open(&src)?); + let mut ar = tar::Archive::new(File::open(src)?); ar.unpack(&dst)?; Self::load_opened(dst) } diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 67bf8956a..4e6347494 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -1,11 +1,35 @@ //! Interface to the eko Python package. use ekore::harmonics::cache::Cache; +use num::Complex; use std::ffi::c_void; pub mod bib; pub mod mellin; +/// Wrapper to pass arguments back to Python +struct RawCmplx { + re: Vec, + im: Vec, +} + +/// Map tensors to c-ordered list +fn unravel(res: Vec<[[Complex; DIM]; DIM]>, order_qcd: usize) -> RawCmplx { + let mut target = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + for obj in res.iter().take(order_qcd) { + 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 +} + /// QCD intergration kernel inside quad. /// /// # Safety @@ -14,27 +38,48 @@ pub mod mellin; 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); - // prepare gamma + // prepare Mellin stuff let path = mellin::TalbotPath::new(u, args.logx, is_singlet); let jac = path.jac() * path.prefactor(); let mut c = Cache::new(path.n()); - let mut re = Vec::::new(); - let mut im = Vec::::new(); - if is_singlet { - let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd( + let mut raw = RawCmplx { + re: Vec::::new(), + im: Vec::::new(), + }; + + if args.is_ome { + if is_singlet { + raw = unravel( + ekore::operator_matrix_elements::unpolarized::spacelike::A_singlet( + args.order_qcd, + &mut c, + args.nf, + args.L, + ), + args.order_qcd, + ); + } else { + raw = unravel( + ekore::operator_matrix_elements::unpolarized::spacelike::A_non_singlet( + args.order_qcd, + &mut c, + args.nf, + args.L, + ), + args.order_qcd, + ); + } + } else if is_singlet { + raw = unravel( + ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd( + args.order_qcd, + &mut c, + args.nf, + ), args.order_qcd, - &mut c, - args.nf, ); - for gamma_s in res.iter().take(args.order_qcd) { - for col in gamma_s.iter().take(2) { - for el in col.iter().take(2) { - re.push(el.re); - im.push(el.im); - } - } - } } else { + // we can not do 1D let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd( args.order_qcd, args.mode0, @@ -42,16 +87,17 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.nf, ); for el in res.iter().take(args.order_qcd) { - re.push(el.re); - im.push(el.im); + raw.re.push(el.re); + raw.im.push(el.im); } } + // pass on (args.py)( - re.as_ptr(), - im.as_ptr(), - c.n.re, - c.n.im, + raw.re.as_ptr(), + raw.im.as_ptr(), + c.n().re, + c.n().im, jac.re, jac.im, args.order_qcd, @@ -72,6 +118,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { args.ev_op_max_order_qcd, args.sv_mode_num, args.is_threshold, + args.Lsv, ) } @@ -101,6 +148,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn( u8, u8, bool, + f64, ) -> f64; /// Additional integration parameters @@ -128,6 +176,8 @@ pub struct QuadQCDargs { pub ev_op_max_order_qcd: u8, pub sv_mode_num: u8, pub is_threshold: bool, + pub is_ome: bool, + pub Lsv: f64, } /// Empty placeholder function for python callback. @@ -159,6 +209,7 @@ pub unsafe extern "C" fn my_py( _ev_op_max_order_qcd: u8, _sv_mode_num: u8, _is_threshold: bool, + _lsv: f64, ) -> f64 { 0. } @@ -166,7 +217,7 @@ 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 `quad_ker_qcd`). +/// package (since it does not appear in the signature of `rust_quad_ker_qcd`). /// /// # Safety /// This is the connection from and back to Python, so we don't know what is on the other side. @@ -193,5 +244,7 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs { ev_op_max_order_qcd: 0, sv_mode_num: 0, is_threshold: false, + is_ome: false, + Lsv: 0., } } diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib index a8a049ac1..0035b7a39 100644 --- a/crates/ekore/refs.bib +++ b/crates/ekore/refs.bib @@ -94,3 +94,16 @@ @article{Buza1996wv pages = "301--320", year = "1998" } +@article{Bierenbaum2009zt, + author = "Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian", + title = "{The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production}", + eprint = "0901.0669", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "DESY-08-187, SFB-CPP-08-107, IFIC-08-68", + doi = "10.1016/j.physletb.2009.01.057", + journal = "Phys. Lett. B", + volume = "672", + pages = "401--406", + year = "2009" +} diff --git a/crates/ekore/src/anomalous_dimensions.rs b/crates/ekore/src/anomalous_dimensions.rs index 6a395d982..55a4dad43 100644 --- a/crates/ekore/src/anomalous_dimensions.rs +++ b/crates/ekore/src/anomalous_dimensions.rs @@ -1,3 +1,3 @@ -//! The anomalous dimensions for DGLAP evolution. +//! The anomalous dimensions for |DGLAP| evolution. pub mod unpolarized; diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs index 71e285ea9..6d25e1c51 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs @@ -1,4 +1,4 @@ -//! LO QCD +//! |LO| |QCD|. use num::complex::Complex; @@ -9,7 +9,7 @@ use crate::harmonics::cache::{Cache, K}; /// /// Implements Eq. (3.4) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0)); CF * gamma @@ -19,7 +19,7 @@ pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex { /// /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0)); 2.0 * TR * 2.0 * (nf as f64) * gamma } @@ -28,7 +28,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0)); 2.0 * CF * gamma } @@ -37,7 +37,7 @@ pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex { /// /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw). pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0); CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64) @@ -54,63 +54,64 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { #[cfg(test)] mod tests { - use crate::cmplx; - use crate::{anomalous_dimensions::unpolarized::spacelike::as1::*, harmonics::cache::Cache}; - use float_cmp::assert_approx_eq; + use super::*; + use crate::harmonics::cache::Cache; + 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.]; + const N: Complex = cmplx!(1., 0.); let mut c = Cache::new(N); let me = gamma_ns(&mut c, NF); - assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); - assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); } #[test] fn quark_momentum_conservation() { - const N: Complex = cmplx![2., 0.]; + const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF); - assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); - assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); } #[test] fn gluon_momentum_conservation() { - const N: Complex = cmplx![2., 0.]; + const N: Complex = cmplx!(2., 0.); let mut c = Cache::new(N); let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF); - assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12); - assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); } #[test] fn gamma_qg_() { - const N: Complex = cmplx![1., 0.]; + const N: Complex = cmplx!(1., 0.); let mut c = Cache::new(N); let me = gamma_qg(&mut c, NF); - assert_approx_eq!(f64, me.re, -20. / 3.0, ulps = 32); - assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12); + assert_approx_eq_cmplx!(f64, me, cmplx!(-20. / 3., 0.), ulps = 32, epsilon = 1e-12); } #[test] fn gamma_gq_() { - const N: Complex = cmplx![0., 1.]; + const N: Complex = cmplx!(0., 1.); let mut c = Cache::new(N); let me = gamma_gq(&mut c, NF); - assert_approx_eq!(f64, me.re, 4. / 3.0, ulps = 32); - assert_approx_eq!(f64, me.im, -4. / 3.0, ulps = 32); + assert_approx_eq_cmplx!(f64, me, cmplx!(4. / 3.0, -4. / 3.0), ulps = 32); } #[test] fn gamma_gg_() { - const N: Complex = cmplx![0., 1.]; + const N: Complex = cmplx!(0., 1.); let mut c = Cache::new(N); let me = gamma_gg(&mut c, NF); - assert_approx_eq!(f64, me.re, 5.195725159621, ulps = 32, epsilon = 1e-11); - assert_approx_eq!(f64, me.im, 10.52008856962, ulps = 32, epsilon = 1e-11); + assert_approx_eq_cmplx!( + f64, + me, + cmplx!(5.195725159621, 10.52008856962), + ulps = 32, + epsilon = 1e-11 + ); } } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs index b2b16459c..f1d0266d5 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs @@ -1,6 +1,6 @@ -//! NLO QCD +//! |NLO| |QCD|. -use ::num::complex::Complex; +use num::complex::Complex; use std::f64::consts::LN_2; use crate::constants::{CA, CF, TR, ZETA2, ZETA3}; @@ -10,7 +10,7 @@ use crate::harmonics::cache::{Cache, K}; /// /// Implements Eq. (3.6) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let Sp1m = c.get(K::S1mh); @@ -50,7 +50,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.5) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let Sp1p = c.get(K::S1h); @@ -90,7 +90,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.6) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let gqqps1_nfcf = (-4. * (2. + N * (5. + N)) * (4. + N * (4. + N * (7. + 5. * N)))) / ((-1. + N) * N.powu(3) * (1. + N).powu(3) * (2. + N).powu(2)); 2.0 * TR * (nf as f64) * CF * gqqps1_nfcf @@ -100,7 +100,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.7) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let Sp2p = c.get(K::S2h); @@ -128,7 +128,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.8) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let Sp2p = c.get(K::S2h); @@ -164,7 +164,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.9) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let Sp1p = c.get(K::S1h); let Sp2p = c.get(K::S2h); @@ -216,11 +216,12 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { #[cfg(test)] mod tests { - use crate::cmplx; - use crate::{anomalous_dimensions::unpolarized::spacelike::as2::*, harmonics::cache::Cache}; - use float_cmp::assert_approx_eq; + use super::*; + use crate::harmonics::cache::Cache; + use crate::{assert_approx_eq_cmplx, cmplx}; use num::complex::Complex; use num::traits::Pow; + use num::Zero; use std::f64::consts::PI; const NF: u8 = 5; @@ -228,29 +229,42 @@ mod tests { #[test] fn physical_constraints() { // number conservation - let mut c = Cache::new(cmplx![1., 0.]); - assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.0, epsilon = 2e-6); + let mut c = Cache::new(cmplx!(1., 0.)); + assert_approx_eq_cmplx!(f64, gamma_nsm(&mut c, NF), Complex::zero(), epsilon = 2e-6); // momentum conservation - let mut c = Cache::new(cmplx![2., 0.]); + let mut c = Cache::new(cmplx!(2., 0.)); let gS1 = gamma_singlet(&mut c, NF); // gluon momentum conservation - assert_approx_eq!(f64, (gS1[0][1] + gS1[1][1]).re, 0.0, epsilon = 4e-5); + assert_approx_eq_cmplx!( + f64, + (gS1[0][1] + gS1[1][1]), + Complex::zero(), + epsilon = 4e-5 + ); // quark momentum conservation - assert_approx_eq!(f64, (gS1[0][0] + gS1[1][0]).re, 0.0, epsilon = 2e-6); + assert_approx_eq_cmplx!( + f64, + (gS1[0][0] + gS1[1][0]), + Complex::zero(), + epsilon = 2e-6 + ); } #[test] fn N2() { // reference values are obtained from MMa - let mut c = Cache::new(cmplx![2., 0.]); + let mut c = Cache::new(cmplx!(2., 0.)); // ns+ - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gamma_nsp(&mut c, NF).re, - (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0, + gamma_nsp(&mut c, NF), + cmplx!( + (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0, + 0. + ), epsilon = 2e-6 ); @@ -259,80 +273,102 @@ mod tests { + (373.0 / 9.0 - 34.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA - 64.0 * (NF as f64) / 27.0) * CF; - assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, check, epsilon = 2e-6); + assert_approx_eq_cmplx!( + f64, + gamma_nsm(&mut c, NF), + cmplx!(check, 0.), + epsilon = 2e-6 + ); // singlet sector let gS1 = gamma_singlet(&mut c, NF); // ps - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gamma_ps(&mut c, NF).re, - -40.0 * CF * (NF as f64) / 27.0 + gamma_ps(&mut c, NF), + cmplx!(-40.0 * CF * (NF as f64) / 27.0, 0.) ); // qg - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gS1[0][1].re, - (-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0 + gS1[0][1], + cmplx!((-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0, 0.) ); // gq - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gS1[1][0].re, - (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0, + gS1[1][0], + cmplx!( + (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0, + 0. + ), epsilon = 1e-13 ); } #[test] fn N3() { - let mut c = Cache::new(cmplx![3., 0.]); + let mut c = Cache::new(cmplx!(3., 0.)); // ns+ let check = ((-34487.0 / 432.0 + 86.0 * PI.pow(2) / 9.0 - 16.0 * ZETA3) * CF + (459.0 / 8.0 - 43.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA - 415.0 * (NF as f64) / 108.0) * CF; - assert_approx_eq!(f64, gamma_nsp(&mut c, NF).re, check, epsilon = 2e-6); + assert_approx_eq_cmplx!( + f64, + gamma_nsp(&mut c, NF), + cmplx!(check, 0.), + epsilon = 2e-6 + ); // singlet sector let gS1 = gamma_singlet(&mut c, NF); // ps - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gamma_ps(&mut c, NF).re, - -1391.0 * CF * (NF as f64) / 5400.0 + gamma_ps(&mut c, NF), + cmplx!(-1391.0 * CF * (NF as f64) / 5400.0, 0.) ); // gq - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gS1[1][0].re, - (973.0 / 432.0 * CF - + (2801.0 / 5400.0 - 7.0 * PI.pow(2) / 9.0) * CA - + 61.0 / 54.0 * (NF as f64)) - * CF + gS1[1][0], + cmplx!( + (973.0 / 432.0 * CF + + (2801.0 / 5400.0 - 7.0 * PI.pow(2) / 9.0) * CA + + 61.0 / 54.0 * (NF as f64)) + * CF, + 0. + ) ); // gg - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gS1[1][1].re, - (-79909.0 / 3375.0 + 194.0 * PI.pow(2) / 45.0 - 8.0 * ZETA3) * CA.pow(2) - - 967.0 / 270.0 * CA * (NF as f64) - + 541.0 / 216.0 * CF * (NF as f64), + gS1[1][1], + cmplx!( + (-79909.0 / 3375.0 + 194.0 * PI.pow(2) / 45.0 - 8.0 * ZETA3) * CA.pow(2) + - 967.0 / 270.0 * CA * (NF as f64) + + 541.0 / 216.0 * CF * (NF as f64), + 0. + ), epsilon = 3e-5 ); } #[test] fn N4() { - let mut c = Cache::new(cmplx![4., 0.]); + let mut c = Cache::new(cmplx!(4., 0.)); // singlet sector let gS1 = gamma_singlet(&mut c, NF); // qg - assert_approx_eq!( + assert_approx_eq_cmplx!( f64, - gS1[0][1].re, - (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * (NF as f64), + gS1[0][1], + cmplx!( + (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * (NF as f64), + 0. + ), epsilon = 1e-14 - ) + ); } } diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs index 14354f736..66d61fc76 100644 --- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs +++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs @@ -1,4 +1,4 @@ -//! NNLO QCD +//! |NNLO| |QCD|. use ::num::complex::Complex; use num::traits::Pow; @@ -11,7 +11,7 @@ use crate::harmonics::cache::{Cache, K}; /// /// Implements Eq. (3.8) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -68,7 +68,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.7) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -125,7 +125,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.9) of [\[Moch:2004pa\]][crate::bib::Moch2004pa]. pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -135,7 +135,7 @@ pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex { let B12 = -(S1 + 1.0 / (N + 1.0) + 1.0 / (N + 2.0)) / (N + 2.0); let B1M = if N.im.abs() < 1.0e-5 && (N - 1.0).re.abs() < 1.0e-5 { - cmplx![-ZETA2, 0.] + cmplx!(-ZETA2, 0.) } else { -(S1 - 1.0 / N) / (N - 1.0) }; @@ -162,7 +162,7 @@ pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.10) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -218,7 +218,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.11) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -275,7 +275,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.12) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -348,7 +348,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { /// /// Implements Eq. (3.13) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw]. pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { - let N = c.n; + let N = c.n(); let S1 = c.get(K::S1); let S2 = c.get(K::S2); let S3 = c.get(K::S3); @@ -425,9 +425,9 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { #[cfg(test)] mod tests { - use crate::cmplx; - use crate::{anomalous_dimensions::unpolarized::spacelike::as3::*, harmonics::cache::Cache}; - use float_cmp::assert_approx_eq; + use super::*; + use crate::harmonics::cache::Cache; + use crate::{assert_approx_eq_cmplx, cmplx}; use num::complex::Complex; const NF: u8 = 5; @@ -435,21 +435,46 @@ mod tests { #[test] fn physical_constraints() { // number conservation - let mut c = Cache::new(cmplx![1., 0.]); - assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, -0.000960586, epsilon = 3e-7); - assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.000594225, epsilon = 6e-7); - - let mut c = Cache::new(cmplx![2., 0.]); + let mut c = Cache::new(cmplx!(1., 0.)); + assert_approx_eq_cmplx!( + f64, + gamma_nsv(&mut c, NF), + cmplx!(-0.000960586, 0.), + epsilon = 3e-7 + ); + assert_approx_eq_cmplx!( + f64, + gamma_nsm(&mut c, NF), + cmplx!(0.000594225, 0.), + epsilon = 6e-7 + ); + + let mut c = Cache::new(cmplx!(2., 0.)); let gS2 = gamma_singlet(&mut c, NF); // gluon momentum conservation - assert_approx_eq!(f64, (gS2[0][1] + gS2[1][1]).re, -0.00388726, epsilon = 2e-6); + assert_approx_eq_cmplx!( + f64, + (gS2[0][1] + gS2[1][1]), + cmplx!(-0.00388726, 0.), + epsilon = 2e-6 + ); // quark momentum conservation - assert_approx_eq!(f64, (gS2[0][0] + gS2[1][0]).re, 0.00169375, epsilon = 2e-6); + assert_approx_eq_cmplx!( + f64, + (gS2[0][0] + gS2[1][0]), + cmplx!(0.00169375, 0.), + epsilon = 2e-6 + ); } #[test] fn N2() { - let mut c = Cache::new(cmplx![2., 0.]); - assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, 188.325593, epsilon = 3e-7); + let mut c = Cache::new(cmplx!(2., 0.)); + assert_approx_eq_cmplx!( + f64, + gamma_nsv(&mut c, NF), + cmplx!(188.325593, 0.), + epsilon = 3e-7 + ); } } diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs index 01b7025f8..3e8a1eb8f 100644 --- a/crates/ekore/src/bib.rs +++ b/crates/ekore/src/bib.rs @@ -1,4 +1,4 @@ -//! List of References (autogenerated on 2024-07-18T17:42:47.650964). +//! List of References (autogenerated on 2024-09-10T11:22:16.645053). #[allow(non_snake_case)] /// The Three loop splitting functions in QCD: The Nonsinglet case @@ -95,3 +95,15 @@ pub fn Ball2015tna() {} /// /// DOI: [10.1007/BF01245820](https:dx.doi.org/10.1007/BF01245820) pub fn Buza1996wv() {} + +#[allow(non_snake_case)] +/// The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production +/// +/// Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian +/// +/// Published in: Phys. Lett. B 672 (2009), 401--406 +/// +/// e-Print: [0901.0669](https://arxiv.org/abs/0901.0669) +/// +/// DOI: [10.1016/j.physletb.2009.01.057](https:dx.doi.org/10.1016/j.physletb.2009.01.057) +pub fn Bierenbaum2009zt() {} diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs index c824bdf5a..8cb90eb1f 100644 --- a/crates/ekore/src/constants.rs +++ b/crates/ekore/src/constants.rs @@ -33,11 +33,11 @@ pub const ZETA3: f64 = 1.2020569031595942; /// $\zeta(4) = \pi^4 / 90$. pub const ZETA4: f64 = 1.082323233711138; -/// singlet-like non-singlet PID +/// singlet-like non-singlet |PID|. pub const PID_NSP: u16 = 10101; -/// valence-like non-singlet PID +/// valence-like non-singlet |PID|. pub const PID_NSM: u16 = 10201; -/// non-singlet all-valence PID +/// non-singlet all-valence |PID|. pub const PID_NSV: u16 = 10200; diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs index eb2a9090b..da0f4800b 100644 --- a/crates/ekore/src/harmonics/cache.rs +++ b/crates/ekore/src/harmonics/cache.rs @@ -30,12 +30,28 @@ pub enum K { S3mh, /// $g_3(N)$ G3, + /// $S_{-1}(N)$ even moments + Sm1e, + /// $S_{-1}(N)$ odd moments + Sm1o, + /// $S_{-2}(N)$ even moments + Sm2e, + /// $S_{-2}(N)$ odd moments + Sm2o, + /// $S_{-3}(N)$ even moments + Sm3e, + /// $S_{-3}(N)$ odd moments + Sm3o, + /// $S_{-2,1}(N)$ even moments + Sm21e, + /// $S_{-2,1}(N)$ odd moments + Sm21o, } /// Hold all cached values. pub struct Cache { /// Mellin N - pub n: Complex, + n: Complex, /// Mapping m: HashMap>, } @@ -49,6 +65,11 @@ impl Cache { } } + /// Get Mellin N. + pub fn n(&self) -> Complex { + Complex::new(self.n.re, self.n.im) + } + /// Retrieve an element. pub fn get(&mut self, k: K) -> Complex { let val = self.m.get(&k); @@ -69,6 +90,14 @@ impl Cache { K::S2mh => w2::S2((self.n - 1.) / 2.), K::S3mh => w3::S3((self.n - 1.) / 2.), K::G3 => g_functions::g3(self.n, self.get(K::S1)), + K::Sm1e => w1::Sm1e(self.get(K::S1), self.get(K::S1h)), + K::Sm1o => w1::Sm1o(self.get(K::S1), self.get(K::S1mh)), + K::Sm2e => w2::Sm2e(self.get(K::S2), self.get(K::S2h)), + K::Sm2o => w2::Sm2o(self.get(K::S2), self.get(K::S2mh)), + K::Sm3e => w3::Sm3e(self.get(K::S3), self.get(K::S3h)), + K::Sm3o => w3::Sm3o(self.get(K::S3), self.get(K::S3mh)), + K::Sm21e => w3::Sm21e(self.n, self.get(K::S1), self.get(K::Sm1e)), + K::Sm21o => w3::Sm21o(self.n, self.get(K::S1), self.get(K::Sm1o)), }; // insert self.m.insert(k, val); @@ -94,15 +123,25 @@ pub fn recursive_harmonic_sum( #[cfg(test)] mod tests { - use crate::harmonics::cache::recursive_harmonic_sum; + use super::*; use crate::harmonics::{w1, w2, w3, w4}; use crate::{assert_approx_eq_cmplx, cmplx}; use num::complex::Complex; + #[test] + fn n() { + let n = cmplx!(1., 0.); + let c = Cache::new(n); + let mut m = c.n(); + m += cmplx!(1., 0.); + assert_approx_eq_cmplx!(f64, c.n(), n); + assert_approx_eq_cmplx!(f64, m, cmplx!(2., 0.)); + } + #[test] fn test_recursive_harmonic_sum() { const SX: [fn(Complex) -> Complex; 4] = [w1::S1, w2::S2, w3::S3, w4::S4]; - const NS: [Complex; 2] = [cmplx![1.0, 0.0], cmplx![2.34, 3.45]]; + const NS: [Complex; 2] = [cmplx!(1.0, 0.0), cmplx!(2.34, 3.45)]; const ITERS: [usize; 2] = [1, 2]; for sit in SX.iter().enumerate() { for nit in NS.iter().enumerate() { diff --git a/crates/ekore/src/harmonics/g_functions.rs b/crates/ekore/src/harmonics/g_functions.rs index cfde9e9ea..9b5348476 100644 --- a/crates/ekore/src/harmonics/g_functions.rs +++ b/crates/ekore/src/harmonics/g_functions.rs @@ -32,19 +32,19 @@ mod tests { #[test] fn test_mellin_g3() { - const NS: [Complex; 3] = [cmplx![1.0, 0.0], cmplx![2.0, 0.0], cmplx![1.0, 1.0]]; + const NS: [Complex; 3] = [cmplx!(1.0, 0.0), cmplx!(2.0, 0.0), cmplx!(1.0, 1.0)]; // NIntegrate[x^({1, 2, 1 + I} - 1) PolyLog[2, x]/(1 + x), {x, 0, 1}] const REFVALS: [Complex; 3] = [ - cmplx![0.3888958462, 0.], - cmplx![0.2560382207, 0.], - cmplx![0.3049381491, -0.1589060625], + cmplx!(0.3888958462, 0.), + cmplx!(0.2560382207, 0.), + cmplx!(0.3049381491, -0.1589060625), ]; for it in NS.iter().enumerate() { let n = *it.1; let s1 = w1::S1(n); let refval = REFVALS[it.0]; let g3 = g3(n, s1); - assert_approx_eq_cmplx![f64, g3, refval, epsilon = 1e-6]; + assert_approx_eq_cmplx!(f64, g3, refval, epsilon = 1e-6); } } } diff --git a/crates/ekore/src/harmonics/polygamma.rs b/crates/ekore/src/harmonics/polygamma.rs index 7415db88e..c9c504adb 100644 --- a/crates/ekore/src/harmonics/polygamma.rs +++ b/crates/ekore/src/harmonics/polygamma.rs @@ -127,73 +127,73 @@ mod tests { #[test] fn fortran() { const ZS: [Complex; 9] = [ - cmplx![1.0, 0.], - cmplx![2.0, 0.], - cmplx![3.0, 0.], - cmplx![0., 1.], - cmplx![-1., 1.], - cmplx![-2., 1.], - cmplx![-1., 2.], - cmplx![-2., 2.], - cmplx![-3., 2.], + cmplx!(1.0, 0.), + cmplx!(2.0, 0.), + cmplx!(3.0, 0.), + cmplx!(0., 1.), + cmplx!(-1., 1.), + cmplx!(-2., 1.), + cmplx!(-1., 2.), + cmplx!(-2., 2.), + cmplx!(-3., 2.), ]; const KS: [usize; 5] = [0, 1, 2, 3, 4]; const FORTRAN_REF: [[Complex; 9]; 5] = [ [ - cmplx![-0.5772156649015332, 0.], - cmplx![0.42278433509846636, 0.], - cmplx![0.9227843350984672, 0.], - cmplx![0.09465032062247669, 2.0766740474685816], - cmplx![0.5946503206224767, 2.5766740474685816], - cmplx![0.9946503206224772, 2.7766740474685814], - cmplx![0.9145915153739776, 2.2208072826422303], - cmplx![1.1645915153739772, 2.47080728264223], - cmplx![1.395360746143208, 2.624653436488384], + cmplx!(-0.5772156649015332, 0.), + cmplx!(0.42278433509846636, 0.), + cmplx!(0.9227843350984672, 0.), + cmplx!(0.09465032062247669, 2.0766740474685816), + cmplx!(0.5946503206224767, 2.5766740474685816), + cmplx!(0.9946503206224772, 2.7766740474685814), + cmplx!(0.9145915153739776, 2.2208072826422303), + cmplx!(1.1645915153739772, 2.47080728264223), + cmplx!(1.395360746143208, 2.624653436488384), ], [ - cmplx![1.6449340668482264, 0.], - cmplx![0.6449340668482264, 0.], - cmplx![0.3949340668482264, 0.], - cmplx![-0.5369999033772361, -0.7942335427593189], - cmplx![-0.5369999033772362, -0.2942335427593189], - cmplx![-0.4169999033772362, -0.13423354275931887], - cmplx![-0.24506883785905695, -0.3178255501472297], - cmplx![-0.24506883785905695, -0.19282555014722969], - cmplx![-0.21548303904248892, -0.12181963298746643], + cmplx!(1.6449340668482264, 0.), + cmplx!(0.6449340668482264, 0.), + cmplx!(0.3949340668482264, 0.), + cmplx!(-0.5369999033772361, -0.7942335427593189), + cmplx!(-0.5369999033772362, -0.2942335427593189), + cmplx!(-0.4169999033772362, -0.13423354275931887), + cmplx!(-0.24506883785905695, -0.3178255501472297), + cmplx!(-0.24506883785905695, -0.19282555014722969), + cmplx!(-0.21548303904248892, -0.12181963298746643), ], [ - cmplx![-2.404113806319188, 0.], - cmplx![-0.40411380631918853, 0.], - cmplx![-0.15411380631918858, 0.], - cmplx![0.3685529315879351, -1.233347149654934], - cmplx![-0.13144706841206477, -0.7333471496549337], - cmplx![-0.09944706841206462, -0.5573471496549336], - cmplx![0.03902435405364951, -0.15743252404131272], - cmplx![-0.02347564594635048, -0.09493252404131272], - cmplx![-0.031668636387861625, -0.053057239562477945], + cmplx!(-2.404113806319188, 0.), + cmplx!(-0.40411380631918853, 0.), + cmplx!(-0.15411380631918858, 0.), + cmplx!(0.3685529315879351, -1.233347149654934), + cmplx!(-0.13144706841206477, -0.7333471496549337), + cmplx!(-0.09944706841206462, -0.5573471496549336), + cmplx!(0.03902435405364951, -0.15743252404131272), + cmplx!(-0.02347564594635048, -0.09493252404131272), + cmplx!(-0.031668636387861625, -0.053057239562477945), ], [ - cmplx![6.49393940226683, 0.], - cmplx![0.49393940226682925, 0.], - cmplx![0.11893940226682913, 0.], - cmplx![4.4771255510465044, -0.31728657866196064], - cmplx![2.9771255510464925, -0.3172865786619599], - cmplx![2.909925551046492, -0.08688657866195917], - cmplx![0.12301766661068443, -0.05523068481179527], - cmplx![0.02926766661068438, -0.055230684811795216], - cmplx![0.004268541930176011, -0.03002148345329936], + cmplx!(6.49393940226683, 0.), + cmplx!(0.49393940226682925, 0.), + cmplx!(0.11893940226682913, 0.), + cmplx!(4.4771255510465044, -0.31728657866196064), + cmplx!(2.9771255510464925, -0.3172865786619599), + cmplx!(2.909925551046492, -0.08688657866195917), + cmplx!(0.12301766661068443, -0.05523068481179527), + cmplx!(0.02926766661068438, -0.055230684811795216), + cmplx!(0.004268541930176011, -0.03002148345329936), ], [ - cmplx![-24.88626612344089, 0.], - cmplx![-0.8862661234408784, 0.], - cmplx![-0.13626612344087824, 0.], - cmplx![3.2795081690440493, 21.41938794863803], - cmplx![0.2795081690440445, 18.419387948637894], - cmplx![-0.012331830955960252, 18.734267948637896], - cmplx![0.14223316576854003, 0.10023607930398608], - cmplx![0.04848316576854002, 0.006486079303986134], - cmplx![0.009893695996688708, 0.014372034600746361], + cmplx!(-24.88626612344089, 0.), + cmplx!(-0.8862661234408784, 0.), + cmplx!(-0.13626612344087824, 0.), + cmplx!(3.2795081690440493, 21.41938794863803), + cmplx!(0.2795081690440445, 18.419387948637894), + cmplx!(-0.012331830955960252, 18.734267948637896), + cmplx!(0.14223316576854003, 0.10023607930398608), + cmplx!(0.04848316576854002, 0.006486079303986134), + cmplx!(0.009893695996688708, 0.014372034600746361), ], ]; for kit in KS.iter().enumerate() { diff --git a/crates/ekore/src/harmonics/w1.rs b/crates/ekore/src/harmonics/w1.rs index a1aa893d9..64431e6d0 100644 --- a/crates/ekore/src/harmonics/w1.rs +++ b/crates/ekore/src/harmonics/w1.rs @@ -10,3 +10,13 @@ use crate::harmonics::polygamma::cern_polygamma; pub fn S1(N: Complex) -> Complex { cern_polygamma(N + 1.0, 0) + 0.577_215_664_901_532_9 } + +/// Analytic continuation of harmonic sum $S_{-1}(N)$ for even moments. +pub fn Sm1e(hS1: Complex, hS1h: Complex) -> Complex { + hS1h - hS1 +} + +/// Analytic continuation of harmonic sum $S_{-1}(N)$ for odd moments. +pub fn Sm1o(hS1: Complex, hS1mh: Complex) -> Complex { + hS1mh - hS1 +} diff --git a/crates/ekore/src/harmonics/w2.rs b/crates/ekore/src/harmonics/w2.rs index d55a32679..a6d7cc89f 100644 --- a/crates/ekore/src/harmonics/w2.rs +++ b/crates/ekore/src/harmonics/w2.rs @@ -11,3 +11,13 @@ use crate::harmonics::polygamma::cern_polygamma; pub fn S2(N: Complex) -> Complex { -cern_polygamma(N + 1.0, 1) + ZETA2 } + +/// Analytic continuation of harmonic sum $S_{-2}(N)$ for even moments. +pub fn Sm2e(hS2: Complex, hS2h: Complex) -> Complex { + 1. / 2. * hS2h - hS2 +} + +/// Analytic continuation of harmonic sum $S_{-2}(N)$ for odd moments. +pub fn Sm2o(hS2: Complex, hS2mh: Complex) -> Complex { + 1. / 2. * hS2mh - hS2 +} diff --git a/crates/ekore/src/harmonics/w3.rs b/crates/ekore/src/harmonics/w3.rs index 7850c678d..bc7286dcb 100644 --- a/crates/ekore/src/harmonics/w3.rs +++ b/crates/ekore/src/harmonics/w3.rs @@ -1,7 +1,10 @@ //! Harmonic sums of weight 3. use num::complex::Complex; +use num::traits::Pow; +use std::f64::consts::LN_2; -use crate::constants::ZETA3; +use crate::constants::{ZETA2, ZETA3}; +use crate::harmonics::g_functions::g3; use crate::harmonics::polygamma::cern_polygamma; /// Compute the harmonic sum $S_3(N)$. @@ -11,3 +14,28 @@ use crate::harmonics::polygamma::cern_polygamma; pub fn S3(N: Complex) -> Complex { 0.5 * cern_polygamma(N + 1.0, 2) + ZETA3 } + +/// Analytic continuation of harmonic sum $S_{-3}(N)$ for even moments. +pub fn Sm3e(hS3: Complex, hS3h: Complex) -> Complex { + 1. / (2.).pow(2) * hS3h - hS3 +} + +/// Analytic continuation of harmonic sum $S_{-3}(N)$ for odd moments. +pub fn Sm3o(hS3: Complex, hS3mh: Complex) -> Complex { + 1. / (2.).pow(2) * hS3mh - hS3 +} + +/// Analytic continuation of harmonic sum $S_{-2,1}(N)$ for even moments. +pub fn Sm21e(N: Complex, hS1: Complex, hSm1: Complex) -> Complex { + Sm21(N, 1., hS1, hSm1) +} + +/// Analytic continuation of harmonic sum $S_{-2,1}(N)$ for odd moments. +pub fn Sm21o(N: Complex, hS1: Complex, hSm1: Complex) -> Complex { + Sm21(N, -1., hS1, hSm1) +} + +/// Analytic continuation of harmonic sum $S_{-2,1}(N)$ for odd moments. +fn Sm21(N: Complex, eta: f64, hS1: Complex, hSm1: Complex) -> Complex { + -eta * g3(N + 1., hS1 + 1. / (N + 1.)) + ZETA2 * hSm1 - 5. / 8. * ZETA3 + ZETA2 * LN_2 +} diff --git a/crates/ekore/src/lib.rs b/crates/ekore/src/lib.rs index af39644cc..db937583d 100644 --- a/crates/ekore/src/lib.rs +++ b/crates/ekore/src/lib.rs @@ -1,4 +1,4 @@ -//! Library for anomalous dimension in DGLAP and transition matrix elements. +//! Library for anomalous dimension in |DGLAP| and |OME|. // Let's stick to the original names which often come from FORTRAN, where such convention do not exists #![allow(non_snake_case)] @@ -7,4 +7,5 @@ pub mod anomalous_dimensions; pub mod bib; mod constants; pub mod harmonics; +pub mod operator_matrix_elements; pub mod util; diff --git a/crates/ekore/src/operator_matrix_elements.rs b/crates/ekore/src/operator_matrix_elements.rs new file mode 100644 index 000000000..98d3cc3b7 --- /dev/null +++ b/crates/ekore/src/operator_matrix_elements.rs @@ -0,0 +1,3 @@ +//! The |OME|. + +pub mod unpolarized; diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized.rs b/crates/ekore/src/operator_matrix_elements/unpolarized.rs new file mode 100644 index 000000000..a00885eec --- /dev/null +++ b/crates/ekore/src/operator_matrix_elements/unpolarized.rs @@ -0,0 +1,3 @@ +//! The unpolarized |OME| for space-like and time-like kinematics. + +pub mod spacelike; diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs new file mode 100644 index 000000000..3c57d33ef --- /dev/null +++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs @@ -0,0 +1,66 @@ +//! The unpolarized, space-like |OME| at various couplings power. +use crate::harmonics::cache::Cache; +use num::complex::Complex; +use num::Zero; +pub mod as1; +pub mod as2; + +/// Compute the tower of the singlet |OME|. +pub fn A_singlet( + matching_order_qcd: usize, + c: &mut Cache, + nf: u8, + L: f64, +) -> Vec<[[Complex; 3]; 3]> { + let mut A_s = vec![ + [ + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero() + ] + ]; + matching_order_qcd + ]; + if matching_order_qcd >= 1 { + A_s[0] = as1::A_singlet(c, nf, L); + } + if matching_order_qcd >= 2 { + // TODO recover MSbar mass + A_s[1] = as2::A_singlet(c, nf, L, false); + } + A_s +} + +/// Compute the tower of the non-singlet |OME|. +pub fn A_non_singlet( + matching_order_qcd: usize, + c: &mut Cache, + nf: u8, + L: f64, +) -> Vec<[[Complex; 2]; 2]> { + let mut A_ns = vec![ + [ + [Complex::::zero(), Complex::::zero()], + [Complex::::zero(), Complex::::zero()] + ]; + matching_order_qcd + ]; + if matching_order_qcd >= 1 { + A_ns[0] = as1::A_ns(c, nf, L); + } + if matching_order_qcd >= 2 { + A_ns[1] = as2::A_ns(c, nf, L); + } + A_ns +} diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs new file mode 100644 index 000000000..960efe43c --- /dev/null +++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs @@ -0,0 +1,142 @@ +//! |NLO| |QCD|. + +use num::complex::Complex; +use num::Zero; + +use crate::cmplx; +use crate::constants::CF; +use crate::harmonics::cache::{Cache, K}; + +/// Compute heavy-heavy |OME|. +/// +/// Implements Eq. (20a) of [\[Ball:2015tna\]](crate::bib::Ball2015tna). +pub fn A_hh(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S1m = c.get(K::S1) - 1. / N; + let S2m = c.get(K::S2) - 1. / N.powu(2); + let ahh_l = (2. + N - 3. * N.powu(2)) / (N * (1. + N)) + 4. * S1m; + let ahh = 2. + * (2. + 5. * N + N.powu(2) + - 6. * N.powu(3) + - 2. * N.powu(4) + - 2. * N * (-1. - 2. * N + N.powu(3)) * S1m) + / (N * (1. + N)).powu(2) + + 4. * (S1m.powu(2) + S2m); + + -CF * (ahh_l * L + ahh) +} + +/// Compute gluon-heavy |OME|. +/// +/// Implements Eq. (20b) of [\[Ball:2015tna\]](crate::bib::Ball2015tna). +pub fn A_gh(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let agh_l1 = (2. + N + N.powu(2)) / (N * (N.powu(2) - 1.)); + let agh_l0 = (-4. + 2. * N + N.powu(2) * (15. + N * (3. + N - N.powu(2)))) + / (N * (N.powu(2) - 1.)).powu(2); + 2. * CF * (agh_l0 + agh_l1 * L) +} + +/// Compute heavy-gluon |OME|. +/// +/// Implements Eq. (B.2) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let den = 1. / (N * (N + 1.) * (2. + N)); + let num = 2. * (2. + N + N.powu(2)); + num * den * L +} + +/// Compute gluon-gluon |OME|. +/// +/// Implements Eq. (B.6) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_gg(_c: &mut Cache, _nf: u8, L: f64) -> Complex { + cmplx!(-2.0 / 3.0 * L, 0.) +} + +/// Compute the |NLO| singlet |OME|. +pub fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex; 3]; 3] { + [ + [A_gg(c, nf, L), Complex::::zero(), A_gh(c, nf, L)], + [ + Complex::::zero(), + Complex::::zero(), + Complex::::zero(), + ], + [A_hg(c, nf, L), Complex::::zero(), A_hh(c, nf, L)], + ] +} + +/// Compute the |NLO| non-singlet |OME|. +pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex; 2]; 2] { + [ + [Complex::::zero(), Complex::::zero()], + [Complex::::zero(), A_hh(c, nf, L)], + ] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::harmonics::cache::Cache; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + const NF: u8 = 5; + + #[test] + fn test_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + const L: f64 = 100.; + let mut c = Cache::new(N); + let aS1 = A_singlet(&mut c, NF, L); + // heavy quark momentum conservation + assert_approx_eq_cmplx!( + f64, + (aS1[0][2] + aS1[1][2] + aS1[2][2]), + Complex::zero(), + epsilon = 1e-10 + ); + // gluon momentum conservation + assert_approx_eq_cmplx!(f64, (aS1[0][0] + aS1[1][0] + aS1[2][0]), Complex::zero()); + } + + #[test] + fn test_A1_intrinsic() { + const N: Complex = cmplx!(2., 0.); + const L: f64 = 3.0; + let mut c = Cache::new(N); + let aNS1i = A_ns(&mut c, NF, L); + let aS1i = A_singlet(&mut c, NF, L); + assert_eq!(aNS1i[1][1], aS1i[2][2]); + } + + #[test] + fn test_Blumlein_1() { + // Test against Blümlein OME implementation Bierenbaum:2009mv. + // Only even moments are available in that code. + // Note there is a minus sign in the definition of L. + const L: f64 = 10.; + let ref_val_gg = [ + -6.666666667, + -6.666666667, + -6.666666667, + -6.666666667, + -6.666666667, + ]; + let ref_val_Hg = [ + 6.666666667, + 3.666666667, + 2.61904761905, + 2.0555555556, + 1.696969697, + ]; + + for n in 0..4 { + let N = cmplx!(2. * (n as f64) + 2., 0.); + let mut c = Cache::new(N); + let aS1 = A_singlet(&mut c, NF, L); + assert_approx_eq_cmplx!(f64, aS1[0][0], cmplx!(ref_val_gg[n], 0.), epsilon = 1e-6); + assert_approx_eq_cmplx!(f64, aS1[2][0], cmplx!(ref_val_Hg[n], 0.), epsilon = 1e-6); + } + } +} diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs new file mode 100644 index 000000000..21b5e7056 --- /dev/null +++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs @@ -0,0 +1,480 @@ +//! |NNLO| |QCD|. + +use num::complex::Complex; +use num::traits::Pow; +use num::Zero; + +use crate::constants::{CA, CF, TR, ZETA2, ZETA3}; +use crate::harmonics::cache::{Cache, K}; + +use crate::operator_matrix_elements::unpolarized::spacelike::as1; + +/// |NNLO| light-light non-singlet |OME|. +/// +/// Implements Eq. (B.4) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + let S1m = c.get(K::S1) - 1. / N; + let S2m = c.get(K::S2) - 1. / N.powu(2); + + let a_qq_l0 = -224.0 / 27.0 * (S1 - 1.0 / N) - 8.0 / 3.0 * ZETA3 + + 40. / 9.0 * ZETA2 + + 73.0 / 18.0 + + 44.0 / 27.0 / N + - 268.0 / 27.0 / (N + 1.0) + + 8.0 / 3.0 * (-1.0 / N.powu(2) + 1.0 / (N + 1.0).powu(2)) + + 20.0 / 9.0 * (S2 - 1.0 / N.powu(2) - ZETA2 + S2 + 1.0 / (N + 1.0).powu(2) - ZETA2) + + 2.0 / 3.0 + * (-2.0 * (S3 - 1.0 / N.powu(3) - ZETA3) + - 2.0 * (S3 + 1.0 / (N + 1.0).powu(3) - ZETA3)); + + let a_qq_l1 = 2. * (-12. - 28. * N + 9. * N.powu(2) + 34. * N.powu(3) - 3. * N.powu(4)) + / (9. * (N * (N + 1.)).powu(2)) + + 80. / 9. * S1m + - 16. / 3. * S2m; + + let a_qq_l2 = -2. * ((2. + N - 3. * N.powu(2)) / (3. * N * (N + 1.)) + 4. / 3. * S1m); + + CF * TR * (a_qq_l2 * L.pow(2) + a_qq_l1 * L + a_qq_l0) +} + +/// |NNLO| heavy-light pure-singlet |OME|. +/// +/// Implements Eq. (B.1) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S2 = c.get(K::S2); + let F1M = 1.0 / (N - 1.0) * (ZETA2 - (S2 - 1.0 / N.powu(2))); + let F11 = 1.0 / (N + 1.0) * (ZETA2 - (S2 + 1.0 / (N + 1.0).powu(2))); + let F12 = 1.0 / (N + 2.0) * (ZETA2 - (S2 + 1.0 / (N + 1.0).powu(2) + 1.0 / (N + 2.0).powu(2))); + let F21 = -F11 / (N + 1.0); + + let a_hq_l0 = -(32.0 / 3.0 / (N - 1.0) + 8.0 * (1.0 / N - 1.0 / (N + 1.0)) + - 32.0 / 3.0 * 1.0 / (N + 2.0)) + * ZETA2 + - 448.0 / 27.0 / (N - 1.0) + - 4.0 / 3.0 / N + - 124.0 / 3.0 * 1.0 / (N + 1.0) + + 1600.0 / 27.0 / (N + 2.0) + - 4.0 / 3.0 * (-6.0 / N.powu(4) - 6.0 / (N + 1.0).powu(4)) + + 2.0 * 2.0 / N.powu(3) + + 10.0 * 2.0 / (N + 1.0).powu(3) + + 16.0 / 3.0 * 2.0 / (N + 2.0).powu(3) + - 16.0 * ZETA2 * (-1.0 / N.powu(2) - 1.0 / (N + 1.0).powu(2)) + + 56.0 / 3.0 / N.powu(2) + + 88.0 / 3.0 / (N + 1.0).powu(2) + + 448.0 / 9.0 / (N + 2.0).powu(2) + + 32.0 / 3.0 * F1M + + 8.0 * ((ZETA2 - S2) / N - F11) + - 32.0 / 3.0 * F12 + + 16.0 * (-(ZETA2 - S2) / N.powu(2) + F21); + + let a_hq_l1 = 8. * (2. + N * (5. + N)) * (4. + N * (4. + N * (7. + 5. * N))) + / ((N - 1.) * (N + 2.).powu(2) * (N * (N + 1.)).powu(3)); + + let a_hq_l2 = + -4. * (2. + N + N.powu(2)).powu(2) / ((N - 1.) * (N + 2.) * (N * (N + 1.)).powu(2)); + + CF * TR * (a_hq_l2 * L.pow(2) + a_hq_l1 * L + a_hq_l0) +} + +/// |NNLO| heavy-gluon |OME|. +/// +/// Implements Eq. (B.3) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +/// The expession for ``A_Hg_l0`` comes form [\[Bierenbaum:2009zt\]](crate::bib::Bierenbaum2009zt). +pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let Sm2e = c.get(K::Sm2e); + let S3 = c.get(K::S3); + let Sm21e = c.get(K::Sm21e); + let Sm3e = c.get(K::Sm3e); + let S1m = S1 - 1. / N; + let S2m = S2 - 1. / N.powu(2); + + #[rustfmt::skip] + let a_hg_l0 = + -( + 3084. + + 192. / N.powu(4) + + 1056. / N.powu(3) + + 2496. / N.powu(2) + + 2928. / N + + 2970. * N + + 1782. * N.powu(2) + + 6. * N.powu(3) + - 1194. * N.powu(4) + - 1152. * N.powu(5) + - 516. * N.powu(6) + - 120. * N.powu(7) + - 12. * N.powu(8) + ) + / ((N - 1.) * ((1. + N) * (2. + N)).powu(4)) + + ( + 764. + - 16. / N.powu(4) + - 80. / N.powu(3) + - 100. / N.powu(2) + + 3. * 72. / N + + 208. * N.powu(3) + + 3. * (288. * N + 176. * N.powu(2) + 16. * N.powu(4)) + ) + / (3. * (1. + N).powu(4) * (2. + N)) + + 12. * Sm3e * (2. + N + N.powu(2)) / (N * (1. + N) * (2. + N)) + - 24. * Sm2e * (4. + N - N.powu(2)) / ((1. + N) * (2. + N)).powu(2) + - S1 + * (48. / N + 432. + 564. * N + 324. * N.powu(2) + 138. * N.powu(3) + 48. * N.powu(4) + 6. * N.powu(5)) + / ((1. + N) * (2. + N)).powu(3) + + S1 + * (-160. - 32. / N.powu(2) - 80. / N + 8. * N * (N - 1.)) + / (3. * (1. + N).powu(2) * (2. + N)) + - 6. * S1.powu(2) * (11. + 8. * N + N.powu(2) + 2. / N) / ((1. + N) * (2. + N)).powu(2) + + 8. * S1.powu(2) * (2. / (3. * N) + 1.) / (N * (2. + N)) + - 2. + * S2 + * (63. + 48. / N.powu(2) + 54. / N + 39. * N + 63. * N.powu(2) + 21. * N.powu(3)) + / ((N - 1.) * (1. + N).powu(2) * (2. + N).powu(2)) + + 8. * S2 * (17. - 2. / N.powu(2) - 5. / N + N * (17. + N)) / (3. * (1. + N).powu(2) * (2. + N)) + + (1. + 2. / N + N) + / ((1. + N) * (2. + N)) + * (24. * Sm2e * S1 + 10. * S1.powu(3) / 9. + 46. * S1 * S2 / 3. + 176. * S3 / 9. - 24. * Sm21e); + + #[rustfmt::skip] + let mut a_hg_l1 = + 2. + * ( + 640. + 2192. * N + + 2072. * N.powu(2) + + 868. * N.powu(3) + + 518. * N.powu(4) + + 736. * N.powu(5) + + 806. * N.powu(6) + + 542. * N.powu(7) + + 228. * N.powu(8) + + 38. * N.powu(9) + ) + / (3. * (N * (N + 1.) * (N + 2.)).powu(3) * (N - 1.)); + + a_hg_l1 -= 2. + * (N * (N.powu(2) - 1.) + * (N + 2.) + * (4. * (36. + N * (88. + N * (33. + N * (8. + 9. * N)))) * S1m + + N * (N + 1.) + * (N + 2.) + * (2. + N + N.powu(2)) + * (10. * S1m.powu(2) + 18. * (2. * Sm2e + ZETA2) + 26. * S2m))) + / (3. * (N * (N + 1.) * (N + 2.)).powu(3) * (N - 1.)); + + a_hg_l1 += 12. * ZETA2 * (-2. + N + N.powu(3)) / (N * (N.powu(2) - 1.) * (N + 2.)); + + #[rustfmt::skip] + let a_hg_l2 = ( + 4. + * (2. + N + N.powu(2)) + * (2. * (-11. + N + N.powu(2)) * (1. + N + N.powu(2)) / (N - 1.)) + / (3. * (N * (N + 1.) * (N + 2.)).powu(2)) + ) + 20. * (2. + N + N.powu(2)) * S1 / (3. * N * (N + 1.) * (N + 2.)); + + a_hg_l2 * L.pow(2) + a_hg_l1 * L + a_hg_l0 +} + +/// |NNLO| gluon-quark |OME|. +/// +/// Implements Eq. (B.5) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S1m = S1 - 1. / N; + + let B2M = ((S1 - 1.0 / N).powu(2) + S2 - 1.0 / N.powu(2)) / (N - 1.0); + let B21 = ((S1 + 1.0 / (N + 1.0)).powu(2) + S2 + 1.0 / (N + 1.0).powu(2)) / (N + 1.0); + + let a_gq_l0 = 4.0 / 3.0 * (2.0 * B2M - 2.0 * (S1.powu(2) + S2) / N + B21) + + 8.0 / 9.0 + * (-10.0 * (S1 - 1.0 / N) / (N - 1.0) + 10.0 * S1 / N + - 8.0 * (S1 + 1.0 / (N + 1.0)) / (N + 1.0)) + + 1.0 / 27.0 * (448.0 * (1.0 / (N - 1.0) - 1.0 / N) + 344.0 / (N + 1.0)); + + let a_gq_l1 = -(-96.0 + 16.0 * N * (7.0 + N * (21.0 + 10.0 * N + 8.0 * N.powu(2))) + - 48.0 * N * (1.0 + N) * (2.0 + N + N.powu(2)) * S1m) + / (9.0 * (N - 1.0) * (N * (1.0 + N)).powu(2)); + + let a_gq_l2 = 8. * (2. + N + N.powu(2)) / (3. * N * (N.powu(2) - 1.)); + + CF * TR * (a_gq_l2 * L.pow(2) + a_gq_l1 * L + a_gq_l0) +} + +/// |NNLO| gluon-gluon |OME|. +/// +/// Implements Eq. (B.7) of [\[Buza:1996wv\]](crate::bib::Buza1996wv). +pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let S1m = S1 - 1. / N; + + let D1 = -1.0 / N.powu(2); + let D11 = -1.0 / (N + 1.0).powu(2); + let D2 = 2.0 / N.powu(3); + let D21 = 2.0 / (N + 1.0).powu(3); + + let a_gg_f = -15.0 - 8.0 / (N - 1.0) + 80.0 / N - 48.0 / (N + 1.0) - 24.0 / (N + 2.0) + + 4.0 / 3.0 * (-6.0 / N.powu(4) - 6.0 / (N + 1.0).powu(4)) + + 6.0 * D2 + + 10.0 * D21 + + 32.0 * D1 + + 48.0 * D11; + + let a_gg_a = -224.0 / 27.0 * (S1 - 1.0 / N) + + 10.0 / 9.0 + + 4.0 / 3.0 * (S1 + 1.0 / (N + 1.0)) / (N + 1.0) + + 1.0 / 27.0 * (556.0 / (N - 1.0) - 628.0 / N + 548.0 / (N + 1.0) - 700.0 / (N + 2.0)) + + 4.0 / 3.0 * (D2 + D21) + + 1.0 / 9.0 * (52.0 * D1 + 88.0 * D11); + + let a_gg_l0 = TR * (CF * a_gg_f + CA * a_gg_a); + + let a_gg_l1 = 8. / 3. + * ((8. + 2. * N + - 34. * N.powu(2) + - 72. * N.powu(3) + - 77. * N.powu(4) + - 37. * N.powu(5) + - 19. * N.powu(6) + - 11. * N.powu(7) + - 4. * N.powu(8)) + / ((N * (N + 1.)).powu(3) * (-2. + N + N.powu(2))) + + 5. * S1m); + + let a_gg_l2 = 4. / 9. + * (1. + + 6. * (2. + N + N.powu(2)).powu(2) / ((N * (N + 1.)).powu(2) * (-2. + N + N.powu(2))) + - 9. * (-4. - 3. * N + N.powu(3)) / (N * (N + 1.) * (-2. + N + N.powu(2)))) + - 4. * S1m; + + a_gg_l2 * L.pow(2) + a_gg_l1 * L + a_gg_l0 +} + +/// |NNLO| singlet |OME|. +pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar_mass: bool) -> [[Complex; 3]; 3] { + let A_hq_2 = A_hq_ps(c, nf, L); + let A_qq_2 = A_qq_ns(c, nf, L); + let mut A_hg_2 = A_hg(c, nf, L); + let A_gq_2 = A_gq(c, nf, L); + let mut A_gg_2 = A_gg(c, nf, L); + + if is_msbar_mass { + A_hg_2 -= 2.0 * 4.0 * CF * as1::A_hg(c, nf, 1.0); + A_gg_2 -= 2.0 * 4.0 * CF * as1::A_gg(c, nf, 1.0); + } + + [ + [A_gg_2, A_gq_2, Complex::::zero()], + [Complex::::zero(), A_qq_2, Complex::::zero()], + [A_hg_2, A_hq_2, Complex::::zero()], + ] +} + +/// |NNLO| non-singlet |OME|. +pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex; 2]; 2] { + [ + [A_qq_ns(c, nf, L), Complex::::zero()], + [Complex::::zero(), Complex::::zero()], + ] +} + +#[cfg(test)] +mod test { + use super::*; + use crate::{assert_approx_eq_cmplx, cmplx}; + + use num::complex::Complex; + const NF: u8 = 5; + + #[test] + fn test_quark_number_conservation() { + let logs = [0., 100.]; + for L in logs { + let N = cmplx!(1., 0.); + let mut c = Cache::new(N); + let aNSqq2 = A_qq_ns(&mut c, NF, L); + assert_approx_eq_cmplx!(f64, aNSqq2, Complex::zero(), epsilon = 2e-11); + } + } + + #[test] + fn test_momentum_conservation() { + let logs = [0., 100.]; + for L in logs { + let N = cmplx!(2., 0.); + let mut c = Cache::new(N); + let aS2 = A_singlet(&mut c, NF, L, false); + + // gluon momenum conservation + assert_approx_eq_cmplx!( + f64, + aS2[0][0] + aS2[1][0] + aS2[2][0], + Complex::zero(), + epsilon = 2e-6 + ); + + // quark momentum conservation + assert_approx_eq_cmplx!( + f64, + aS2[0][1] + aS2[1][1] + aS2[2][1], + Complex::zero(), + epsilon = 1e-11 + ); + + // additional test, maybe to be moved in a different function + + assert_approx_eq_cmplx!(f64, aS2[1][0], Complex::zero()); + + let aNS2 = A_ns(&mut c, NF, L); + + assert_approx_eq_cmplx!(f64, aNS2[1][1], Complex::zero()); + + for el in aS2.iter() { + assert_approx_eq_cmplx!(f64, el[2], Complex::zero()); + } + } + } + + #[test] + fn pegasus_sign() { + let ref_val = cmplx!(-21133.9, 0.); + let N = cmplx!(2., 0.); + let mut c = Cache::new(N); + let L = 100.; + + let aS2 = A_singlet(&mut c, NF, L, false); + // add macro for this? + assert_approx_eq_cmplx!(f64, aS2[0][0], ref_val, rel = 4e-5); + } + + #[test] + fn Blumlein_2() { + let ref_val_gg = [ + [-9.96091, -30.0093, -36.5914, -40.6765, -43.6823], + [-289.097, -617.811, -739.687, -820.771, -882.573], + ]; + + let ref_val_gq = [ + [5.82716, 1.34435, 0.739088, 0.522278, 0.413771], + [191.506, 60.3468, 37.0993, 27.3308, 21.8749], + ]; + + let ref_val_Hg = [ + [9.96091, 10.6616, 9.27572, 8.25694, 7.49076], + [289.097, 278.051, 223.261, 186.256, 160.027], + ]; + + let ref_val_Hq = [ + [-2.66667, -0.365962, -0.15071, -0.084247, -0.0543195], + [-101.432, -17.2316, -7.25937, -4.04249, -2.58706], + ]; + + let ref_val_qq = [ + [ + -3.16049, -5.0571, -6.43014, -7.51258, -8.409, -9.17546, -9.84569, -10.4416, + -10.9783, + ], + [ + -90.0741, -139.008, -173.487, -200.273, -222.222, -240.833, -256.994, -271.28, + -284.083, + ], + ]; + + for n in 2..11 { + for (i, L) in [0., 10.].iter().enumerate() { + let N = cmplx!(n as f64, 0.); + let mut c = Cache::new(N); + + let aS2 = A_singlet(&mut c, NF, *L, false); + if n % 2 == 0 { + let idx = n / 2 - 1; + assert_approx_eq_cmplx!( + f64, + aS2[0][0], + cmplx!(ref_val_gg[i][idx], 0.), + rel = 2e-6 + ); + assert_approx_eq_cmplx!( + f64, + aS2[0][1], + cmplx!(ref_val_gq[i][idx], 0.), + rel = 4e-6 + ); + assert_approx_eq_cmplx!( + f64, + aS2[2][0], + cmplx!(ref_val_Hg[i][idx], 0.), + rel = 3e-6 + ); + assert_approx_eq_cmplx!( + f64, + aS2[2][1], + cmplx!(ref_val_Hq[i][idx], 0.), + rel = 3e-6 + ); + } + assert_approx_eq_cmplx!( + f64, + aS2[1][1], + cmplx!(ref_val_qq[i][n - 2], 0.), + rel = 4e-6 + ); + } + } + } + + #[test] + fn Hg2_pegasus() { + let L = 0.; + for n in 3..20 { + let N = cmplx!(n as f64, 0.); + + let mut c = Cache::new(N); + + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let S3 = c.get(K::S3); + + let aS2 = A_singlet(&mut c, NF, L, false); + let E2 = 2.0 / N * (ZETA3 - S3 + 1.0 / N * (ZETA2 - S2 - S1 / N)); + + let a_hg_param = -0.006 + 1.111 * (S1.powu(3) + 3.0 * S1 * S2 + 2.0 * S3) / N + - 0.400 * (S1.powu(2) + S2) / N + + 2.770 * S1 / N + - 24.89 / (N - 1.0) + - 187.8 / N + + 249.6 / (N + 1.0) + + 1.556 * 6.0 / N.powu(4) + - 3.292 * 2.0 / N.powu(3) + + 93.68 * 1.0 / N.powu(2) + - 146.8 * E2; + + assert_approx_eq_cmplx!(f64, aS2[2][0], a_hg_param, rel = 7e-4); + } + } + + #[test] + fn msbar_matching() { + let logs = [0., 100.]; + for L in logs { + let N = cmplx!(2., 0.); + let mut c = Cache::new(N); + let aS2 = A_singlet(&mut c, NF, L, false); + assert_approx_eq_cmplx!( + f64, + aS2[0][0] + aS2[1][0] + aS2[2][0], + Complex::zero(), + epsilon = 2e-6 + ); + } + } +} diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs index 4b80263d6..b5c812165 100644 --- a/crates/ekore/src/util.rs +++ b/crates/ekore/src/util.rs @@ -8,23 +8,17 @@ macro_rules! cmplx { }; } -/// Shorthand complex number contructor. +/// Shorthand complex number comparators. #[cfg(test)] #[macro_export] macro_rules! assert_approx_eq_cmplx { - ($size:ty, $ref:expr, $target:expr) => { - use float_cmp::assert_approx_eq; - assert_approx_eq!($size, $ref.re, $target.re); - assert_approx_eq!($size, $ref.im, $target.im); - }; - ($size:ty, $ref:expr, $target:expr, ulps=$ulps:expr) => { - use float_cmp::assert_approx_eq; - assert_approx_eq!($size, $ref.re, $target.re, ulps = $ulps); - assert_approx_eq!($size, $ref.im, $target.im, ulps = $ulps); + ($size:ty, $ref:expr, $target:expr, rel=$rel:expr) => { + assert!($target.norm() > 0.0, "target has norm=0!"); + float_cmp::assert_approx_eq!($size, $ref.re, $target.re, epsilon = $rel * $target.norm()); + float_cmp::assert_approx_eq!($size, $ref.im, $target.im, epsilon = $rel * $target.norm()); }; - ($size:ty, $ref:expr, $target:expr, epsilon=$epsilon:expr) => { - use float_cmp::assert_approx_eq; - assert_approx_eq!($size, $ref.re, $target.re, epsilon = $epsilon); - assert_approx_eq!($size, $ref.im, $target.im, epsilon = $epsilon); + ($size:ty, $ref:expr, $target:expr $(, $set:ident = $val:expr)*) => { + float_cmp::assert_approx_eq!($size, $ref.re, $target.re $(, $set = $val)*); + float_cmp::assert_approx_eq!($size, $ref.im, $target.im $(, $set = $val)*); }; } diff --git a/rustify.sh b/rustify.sh index 297afdb3d..bb4ae2565 100755 --- a/rustify.sh +++ b/rustify.sh @@ -6,4 +6,9 @@ patch -p1 src/eko/evolution_operator/__init__.py.patch patch -p1 src/eko/evolution_operator/operator_matrix_element.py.patch +patch -p1 = 1: +@@ -102,105 +97,6 @@ def build_ome(A, matching_order, a_s, backward_method): + return ome + + +-@nb.njit(cache=True) +-def quad_ker( +- u, +- order, +- mode0, +- mode1, +- is_log, +- logx, +- areas, +- a_s, +- nf, +- L, +- sv_mode, +- Lsv, +- backward_method, +- is_msbar, +- is_polarized, +- is_time_like, +-): +- r"""Raw kernel inside quad. +- +- Parameters +- ---------- +- u : float +- quad argument +- order : tuple(int,int) +- perturbation matching order +- mode0 : int +- pid for first element in the singlet sector +- mode1 : int +- pid for second element in the singlet sector +- is_log : boolean +- logarithmic interpolation +- logx : float +- Mellin inversion point +- areas : tuple +- basis function configuration +- a_s : float +- strong coupling, needed only for the exact inverse +- nf: int +- number of active flavor below threshold +- L : float +- :math:``\ln(\mu_F^2 / m_h^2)`` +- backward_method : InversionMethod or None +- empty or method for inverting the matching condition (exact or expanded) +- is_msbar: bool +- add the |MSbar| contribution +- is_polarized : boolean +- is polarized evolution ? +- is_time_like : boolean +- is time-like evolution ? +- +- Returns +- ------- +- ker : float +- evaluated integration kernel +- """ +- ker_base = QuadKerBase(u, is_log, logx, mode0) +- integrand = ker_base.integrand(areas) +- if integrand == 0.0: +- return 0.0 +- # compute the ome +- if ker_base.is_singlet or ker_base.is_QEDsinglet: +- indices = {21: 0, 100: 1, 90: 2} +- if is_polarized: +- if is_time_like: +- raise NotImplementedError("Polarized, time-like is not implemented") +- A = ome_ps.A_singlet(order, ker_base.n, nf, L) +- else: +- if is_time_like: +- A = ome_ut.A_singlet(order, ker_base.n, L) +- else: +- A = ome_us.A_singlet(order, ker_base.n, nf, L, is_msbar) +- else: +- indices = {200: 0, 91: 1} +- if is_polarized: +- if is_time_like: +- raise NotImplementedError("Polarized, time-like is not implemented") +- A = ome_ps.A_non_singlet(order, ker_base.n, L) +- else: +- if is_time_like: +- A = ome_ut.A_non_singlet(order, ker_base.n, L) +- else: +- A = ome_us.A_non_singlet(order, ker_base.n, nf, L) +- +- # correct for scale variations +- if sv_mode == sv.Modes.exponentiated: +- A = gamma_variation(A, order, nf, Lsv) +- +- # build the expansion in alpha_s depending on the strategy +- ker = build_ome(A, order, a_s, backward_method) +- +- # select the needed matrix element +- ker = ker[indices[mode0], indices[mode1]] +- +- # recombine everything +- return np.real(ker * integrand) +- +- + class OperatorMatrixElement(Operator): + r"""Internal representation of a single |OME|. + +@@ -290,41 +186,14 @@ class OperatorMatrixElement(Operator): + ) + return labels + +- def quad_ker(self, label, logx, areas): +- """Return partially initialized integrand function. +- +- Parameters +- ---------- +- label: tuple +- operator element pids +- logx: float +- Mellin inversion point +- areas : tuple +- basis function configuration +- +- Returns +- ------- +- functools.partial +- partially initialized integration kernel +- """ +- return functools.partial( +- quad_ker, +- order=self.order, +- mode0=label[0], +- mode1=label[1], +- is_log=self.int_disp.log, +- logx=logx, +- areas=areas, +- a_s=self.a_s, +- nf=self.nf, +- L=self.L, +- sv_mode=self.sv_mode, +- Lsv=np.log(self.xif2), +- backward_method=self.backward_method, +- is_msbar=self.is_msbar, +- is_polarized=self.config["polarized"], +- is_time_like=self.config["time_like"], +- ) ++ def update_cfg(self, cfg): ++ """Adjust integration config.""" ++ cfg.is_ome = True ++ cfg.py = ekors.ffi.cast("void *", cb_quad_ker_ome.address) ++ cfg.L = self.L ++ cfg.as1 = self.a_s ++ cfg.as0 = 0.0 ++ cfg.Lsv = np.log(self.xif2) + + @property + def a_s(self): diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index a8ba8788c..45f947a8b 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -7,6 +7,7 @@ from .. import interpolation from .. import scale_variations as sv +from ..io.types import InversionMethod from ..kernels import non_singlet as ns from ..kernels import singlet as s from ..scale_variations import expanded as sv_expanded @@ -38,33 +39,37 @@ def select_singlet_element(ker, mode0, mode1): return ker[j, k] +CB_SIGNATURE = nb.types.double( + nb.types.CPointer(nb.types.double), # re_*_raw + nb.types.CPointer(nb.types.double), # im_*_raw + nb.types.double, # re_n + nb.types.double, # im_n + nb.types.double, # re_jac + nb.types.double, # im_jac + nb.types.uintc, # order_qcd + nb.types.bool_, # is_singlet + nb.types.uintc, # mode0 + nb.types.uintc, # mode1 + nb.types.uintc, # nf + nb.types.bool_, # is_log + nb.types.double, # logx + 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 + nb.types.uintc, # ev_op_iterations + nb.types.uintc, # ev_op_max_order_qcd + nb.types.uintc, # sv_mode_num + nb.types.bool_, # is_threshold + nb.types.double, # Lsv +) + + @nb.cfunc( - nb.types.double( - nb.types.CPointer(nb.types.double), # re_gamma_ns_raw - nb.types.CPointer(nb.types.double), # im_gamma_ns_raw - nb.types.double, # re_n - nb.types.double, # im_n - nb.types.double, # re_jac - nb.types.double, # im_jac - nb.types.uint, # order_qcd - nb.types.bool_, # is_singlet - nb.types.uint, # mode0 - nb.types.uint, # mode1 - nb.types.uint, # nf - nb.types.bool_, # is_log - nb.types.double, # logx - nb.types.CPointer(nb.types.double), # areas_raw - nb.types.uint, # areas_x - nb.types.uint, # areas_y - nb.types.double, # L - nb.types.uint, # method_num - nb.types.double, # as1 - nb.types.double, # as0 - nb.types.uint, # ev_op_iterations - nb.types.uint, # ev_op_max_order_qcd - nb.types.uint, # sv_mode_num - nb.types.bool_, # is_threshold - ), + CB_SIGNATURE, cache=True, nopython=True, ) @@ -93,6 +98,7 @@ def cb_quad_ker_qcd( ev_op_max_order_qcd, _sv_mode_num, is_threshold, + _Lsv, ): """C Callback inside integration kernel.""" # recover complex variables @@ -156,6 +162,132 @@ def cb_quad_ker_qcd( return np.real(res) +@nb.njit(cache=True) +def build_ome(A, matching_order, a_s, backward_method): + r"""Construct the matching expansion in :math:`a_s` with the appropriate + method. + + Parameters + ---------- + A : numpy.ndarray + list of |OME| + matching_order : tuple(int,int) + perturbation matching order + a_s : float + strong coupling, needed only for the exact inverse + backward_method : InversionMethod or None + empty or method for inverting the matching condition (exact or expanded) + + Returns + ------- + ome : numpy.ndarray + matching operator matrix + """ + # to get the inverse one can use this FORM snippet + # Symbol a; + # NTensor c,d,e; + # Local x=-(a*c+a**2* d + a**3 * e); + # Local bi = 1+x+x**2+x**3; + # Print; + # .end + ome = np.eye(len(A[0]), dtype=np.complex_) + A = A[:, :, :] + A = np.ascontiguousarray(A) + if backward_method is InversionMethod.EXPANDED: + # expended inverse + if matching_order[0] >= 1: + ome -= a_s * A[0] + if matching_order[0] >= 2: + ome += a_s**2 * (-A[1] + A[0] @ A[0]) + if matching_order[0] >= 3: + ome += a_s**3 * (-A[2] + A[0] @ A[1] + A[1] @ A[0] - A[0] @ A[0] @ A[0]) + else: + # forward or exact inverse + if matching_order[0] >= 1: + ome += a_s * A[0] + if matching_order[0] >= 2: + ome += a_s**2 * A[1] + if matching_order[0] >= 3: + ome += a_s**3 * A[2] + # need inverse exact ? so add the missing pieces + if backward_method is InversionMethod.EXACT: + ome = np.linalg.inv(ome) + return ome + + +@nb.cfunc( + CB_SIGNATURE, + cache=True, + nopython=True, +) +def cb_quad_ker_ome( + re_ome_raw, + im_ome_raw, + re_n, + im_n, + re_jac, + im_jac, + order_qcd, + is_singlet, + mode0, + mode1, + nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, + L, + _method_num, + as1, + _as0, + _ev_op_iterations, + _ev_op_max_order_qcd, + _sv_mode_num, + _is_threshold, + Lsv, +): + """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) + # TODO recover parameters + sv_mode = sv.Modes.exponentiated + order = (order_qcd, 0) + if is_singlet: + indices = {21: 0, 100: 1, 90: 2} + # reconstruct singlet matrices + re_ome_singlet = nb.carray(re_ome_raw, (order_qcd, 3, 3)) + im_ome_singlet = nb.carray(im_ome_raw, (order_qcd, 3, 3)) + A = re_ome_singlet + im_ome_singlet * 1j + else: + indices = {200: 0, 91: 1} + # construct non-singlet matrices + re_ome_ns = nb.carray(re_ome_raw, (order_qcd, 2, 2)) + im_ome_ns = nb.carray(im_ome_raw, (order_qcd, 2, 2)) + A = re_ome_ns + im_ome_ns * 1j + + # correct for scale variations + if sv_mode == sv.Modes.exponentiated: + A = sv.exponentiated.gamma_variation(A, order, nf, Lsv) + + # TODO recover InversionMethod + backward_method = "" + + # build the expansion in alpha_s depending on the strategy + ker = build_ome(A, order, as1, backward_method) + + # select the needed matrix element + ker = ker[indices[mode0], indices[mode1]] + + # recombine everything + res = ker * pj * jac + 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