Theory citation function (#123)

* Copy over PyBaMM citations code

Also add license notice

* Minimal working version of citations function

By default, just cite the flavio paper.
Removed the read_citation method, as rather than reading the bibtex
from a file, it seems better to just ouput the INSPIRE texkeys and
worry about getting bibtex from INSPIRE later.
(E.g. if you print to a .tex file, you could use INSPIRE's bibliography
generator at https://inspirehep.net/bibliography-generator, or if fixed
David's inspiretools python package.)

* A initial citation reference for testing

The way flavio does Gamma_12 is based on the a,b,c notation originally
introduced in the 2003 BBLN paper

* Add SM_citations function to Observable class

Gives the theory papers to be cited for an observable, using a throwaway
instance.
As part of this, move the registering of the flavio paper out of the
constructor so you don't get that as well.

* Merge _reset method into constructor

Since this only gets used here, just do it in the constructor.
(In the original PyBaMM version the _reset just gets used for testing,
which I think we want to do slightly differently anyway.)

* Add docstring to SM_citations

* Copy over the PyBaMM citation test class

* Revert "Merge _reset method into constructor"

This reverts commit a6a545e.
Yeah, we need the reset method for testing...

* Rejig the print functions to return a list

Return something instead of printing directly, unless outputting to a
file.

* Add some tests of the citation functionality

* Fix typo in docstring

* Citations for B->Xgamma

Note the CP asymmetry needs citations adding, wherever the numbers in
ka1_r, ka2_r, ka1_i came from.

* Citations for meson mixing

Just the Inami-Lim function - everything else worth citing comes in
through parameters (fB, B, RG factors)

* Add citations for W and Z observables

* Improve theory citation function

Allow for arguments to be passed to prediction function.
Also change the name as you are really getting the citations for a
general theory prediction, not just the SM.

* Add citations for beta decay

* Fix failing test

Forgot to update the test when I renamed the function

* Citations for mu and tau decays

* Citations for quark mass conversions

* Citations for K decays

* Add D decay form factor citations

* Add a bunch of b decay citations

* Add bvll citations

* Add B -> P formfactor citations

* Add B -> gamma formfactor citations

* Add citations for B->V formfactors

* Small change to docstrings

* Add citations for Higgs stuff

I'm assuming all the numbers come (indirectly) from 1911.07866 as
mentioned in the release notes for flavio v2

* Adds citations for B->llgamma

* Fix indentation

Got borked when I rebased

* Move copyright notice to appropriate file

* Give citations instance a unique name

* Correct docstring

* Fix code typo

Somehow this find+replace got mucked up

* Improve internals of the citation class

Update tests and the theory_citation function too

* Remove tex citation string method, add properties

* Add convenience `register_citation` method

* use multiprocessing.Array

* raise meaningful error meassge if inspire key not in YAML

* Add docstrings for the main methods

* Update and add tests

Specifically test that a multithreaded computation works and gives the
same citations as single threaded.

* Update tests

Remove the test for a "unknown" inspire key, since it becomes a "known"
inspirekey as soon as we re-run update_citations.py.

Add a test to check that the list of "known" citations matches what's in
the source code.

* improve regexp in `extract_citations` function

Co-authored-by: Peter Stangl <peter.stangl@ph.tum.de>

flavio/__init__.py

@@ -7,6 +7,7 @@
 from . import classes
 from .classes import Measurement, Parameter, ParameterConstraints, Observable, NamedInstanceClass
 from .config import config
+from . import citations
 from flavio.physics.eft import WilsonCoefficients
 from flavio.parameters import default_parameters
 from flavio.functions import sm_prediction, sm_uncertainty, np_uncertainty, sm_error_budget, np_prediction, sm_covariance, combine_measurements

flavio/citations.py

@@ -0,0 +1,86 @@
+# Inspired by pybamm.citations
+# https://github.com/pybamm-team/PyBaMM
+
+from multiprocessing import Array
+import flavio
+import sys
+from itertools import compress
+import ctypes
+import yaml
+import pkgutil
+
+
+class CitationScope:
+
+    def __enter__(self):
+        self._citations_global = flavio.citations
+        flavio.citations = Citations()
+        return flavio.citations
+
+    def __exit__(self, type, value, traceback):
+        flavio.citations = self._citations_global
+
+
+class Citations:
+
+    collect = CitationScope
+    __all_citations = set(yaml.safe_load(
+        pkgutil.get_data('flavio', 'data/citations.yml')
+    ))
+
+    def __init__(self, initial_citations=[]):
+        self._initial_citations = set(initial_citations)
+        self._all_citations = {
+            k:i for i,k in
+            enumerate(sorted(self.__all_citations | self._initial_citations))
+        }
+        self._array = Array(ctypes.c_bool, len(self._all_citations))
+        for inspire_key in self._initial_citations:
+            self.register(inspire_key)
+
+    def __iter__(self):
+        for citation in self.set:
+            yield citation
+
+    def __str__(self):
+        return ",".join(self.set)
+
+    @property
+    def set(self):
+        return set(compress(sorted(self._all_citations.keys()), self._array))
+
+    @property
+    def string(self):
+        return str(self)
+
+    def register(self, inspire_key):
+        """Register a paper to be cited. The intended use is that this method
+        should be called only when the referenced functionality is actually being used.
+        Parameters
+        ----------
+        key : str
+            The INSPIRE texkey for the paper to be cited
+        """
+        try:
+            self._array[self._all_citations[inspire_key]] = True
+        except KeyError:
+            from flavio.util import get_datapath
+            yaml_path = get_datapath('flavio', 'data/citations.yml')
+            raise KeyError(
+                f'The inspire key must be contained in {yaml_path}. '
+                f'The key `{inspire_key}` was not found there.'
+            )
+
+
+    def clear(self):
+        """Clear the list of cited papers (including any default citations)."""
+        self._array[:] = [False]*len(self._array)
+
+    def reset(self):
+        """Reset the list of cited papers back to only the default ones."""
+        self.clear()
+        for inspire_key in self._initial_citations:
+            self.register(inspire_key)
+
+
+sys.modules[__name__] = Citations(["Straub:2018kue"])

flavio/classes.py

@@ -697,6 +697,17 @@ def get_measurements(self):
                             break
         return ms
 
+    def theory_citations(self, *args, **kwargs):
+        """Return a set of theory papers (in the form of INSPIRE texkeys) to
+        cite for the theory prediction for an observable.
+
+        Arguments are passed to the observable and are necessary,
+        depending on the observable (e.g. $q^2$-dependent observables).
+        """
+        with flavio.citations.collect() as citations:
+            flavio.sm_prediction(self.name, *args, **kwargs)
+        return citations.set
+
 
 
 class AuxiliaryQuantity(NamedInstanceClass):

flavio/data/citations.yml

@@ -0,0 +1,42 @@
+- Antonelli:2010yf
+- Asatrian:2003vq
+- Awramik:2003rn
+- Beneke:2001at
+- Beneke:2003az
+- Beneke:2004dp
+- Benson:2003kp
+- Bernard:2009zm
+- Bjorn:2016zlr
+- Bobeth:2013uxa
+- Boer:2014kda
+- Bourrely:2008za
+- Brignole:2004ah
+- Brivio:2017vri
+- Brod:2010hi
+- Buras:2006gb
+- Buras:2015yba
+- Caprini:1997mu
+- Chobanova:2017rkj
+- Cirigliano:2007ga
+- Descotes-Genon:2015hea
+- Falkowski:2019hvp
+- Freitas:2014hra
+- Gambino:2011cq
+- Gonzalez-Alonso:2018omy
+- Gorbahn:2006bm
+- Greub:2008cy
+- Guadagnoli:2016erb
+- Gubernari:2018wyi
+- Huber:2015sra
+- Inami:1980fz
+- Kitano:2002mt
+- Kozachuk:2017mdk
+- Kruger:2002gf
+- Kuno:1999jp
+- Melikhov:2004mk
+- Misiak:2006ab
+- Misiak:2015xwa
+- Pich:2013lsa
+- Seidel:2004jh
+- Straub:2015ica
+- fakename:2020abc

flavio/data/update_citations.py

@@ -0,0 +1,25 @@
+#! /usr/bin/env python3
+
+"""Standalone helper script to update the list of citations referenced in the
+flavio source code."""
+
+import argparse
+import logging
+import sys
+import yaml
+import pkgutil
+logging.basicConfig(level=logging.INFO)
+
+def main(argv):
+    parser = argparse.ArgumentParser(description='Update the list of citations referenced in the flavio source code.')
+    args = parser.parse_args()
+
+    from flavio.util import get_datapath, extract_citations
+
+    filename = get_datapath('flavio', 'data/citations.yml')
+    with open(filename, 'w') as f:
+        f.write(yaml.dump(sorted(extract_citations())))
+        logging.info(f"Saved updated citations to {filename}")
+
+if __name__ == '__main__':
+    main(sys.argv)

flavio/physics/bdecays/bllgamma.py

@@ -52,7 +52,8 @@ def getfft(s, par, B, ff, ff0, lep, wc):
 
     #Add light meson resonances
     resonances = {'Bs': ['phi'], 'B0': ['rho0', 'omega']}
-    fVtofVemfactors = {'phi': -1/3, 'rho0': 1/sqrt(2), 'omega': 1/(3*sqrt(2))}
+    flavio.citations.register("Kozachuk:2017mdk")
+    fVtofVemfactors = {'phi': -1/3, 'rho0': 1/sqrt(2), 'omega': 1/(3*sqrt(2))} # Given in Sec. 8.A.3 of 1712.07926
     #We use the general B->rho and B->omega parameters, hence the isotopic factors
     resgV = {('Bs','phi'): -par['Bs->phi BSZ a0_T1'],
              ('B0','rho0'): 1/sqrt(2)*par['B->rho BSZ a0_T1'],
@@ -66,6 +67,7 @@ def getfft(s, par, B, ff, ff0, lep, wc):
     return (fftv, ffta)
 
 def getF1(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
     scale = config['renormalization scale']['bllgamma']
     mb = running.get_mb(par, scale, nf_out=5)
     mbh = mb/par['m_'+B]
@@ -75,6 +77,7 @@ def getF1(s, par, B, ff, ff0, lep, wc):
     return (abs(wc['C9_'+label])**2 + abs(wc['C10_'+label])**2)*ff['v']**2 + 4*mbh**2/s**2*abs(wc['C7_'+bq]*fftv)**2 + 4*mbh/s*ff['v']*_Re(wc['C7_'+bq]*ffta*_Co(wc['C9_'+label]))
 
 def getF2(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
     scale = config['renormalization scale']['bllgamma']
     mb = running.get_mb(par, scale, nf_out=5)
     mbh = mb/par['m_'+B]
@@ -85,6 +88,8 @@ def getF2(s, par, B, ff, ff0, lep, wc):
 
 
 def B10(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
+    flavio.citations.register("Guadagnoli:2016erb")
     F1 = getF1(s, par, B, ff, ff0, lep, wc)
     F2 = getF2(s, par, B, ff, ff0, lep, wc)
     mlh = par['m_'+lep]/par['m_'+B]
@@ -94,11 +99,15 @@ def B10(s, par, B, ff, ff0, lep, wc):
 # B11 vanishes for the BR estimation
 
 def B12(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
+    flavio.citations.register("Guadagnoli:2016erb")
     F1 = getF1(s, par, B, ff, ff0, lep, wc)
     F2 = getF2(s, par, B, ff, ff0, lep, wc)
     return s*(F1+F2)
 
 def B120(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
+    flavio.citations.register("Guadagnoli:2016erb")
     scale = config['renormalization scale']['bllgamma']
     mb = running.get_mb(par, scale, nf_out=5)
     mbh = mb/par['m_'+B]
@@ -110,18 +119,21 @@ def B120(s, par, B, ff, ff0, lep, wc):
 
 
 def dG1dsMN(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
     mlh = par['m_'+lep]/par['m_'+B]
     pref = prefactor(s, par, B, ff, lep, wc)*(1-s)**3*sqrt(1-4*mlh**2/s)
     return pref*(B10(s, par, B, ff, ff0, lep, wc) + (s - 4*mlh**2)/(3*s)*B12(s, par, B, ff, ff0, lep, wc))
 
 def dG2dsMN(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
     mlh = par['m_'+lep]/par['m_'+B]
     bq = meson_quark[B]
     label = bq+lep+lep
     pref = 16*prefactor(s, par, B, ff, lep, wc)*(par['f_'+B]/par['m_'+B]*abs(wc['C10_'+label])*mlh)**2/(1-s)
     return pref*(-8*sqrt(s-4*mlh**2)*sqrt(s)+4*atanh(sqrt(1-4*mlh**2/s))*(1+s-mlh**2*(1-s)**2))
 
 def dG12dsMN(s, par, B, ff, ff0, lep, wc):
+    flavio.citations.register("Melikhov:2004mk")
     mlh = par['m_'+lep]/par['m_'+B]
     pref = 4*prefactor(s, par, B, ff, lep, wc)*par['f_'+B]/par['m_'+B]*(1-s)
     return pref*atanh(sqrt(1-4*mlh**2/s))*B120(s, par, B, ff, ff0, lep, wc)

flavio/physics/bdecays/bvll/observables_bs.py

@@ -34,6 +34,7 @@ def bsvll_obs(function, q2, wc_obj, par, B, V, lep):
 
 def S_theory_num_Bs(y, J, J_bar, J_h, i):
     # (42) of 1502.05509
+    flavio.citations.register("Descotes-Genon:2015hea")
     return 1/(1-y**2) * (J[i] + J_bar[i]) - y/(1-y**2) * J_h[i]
 
 def S_experiment_num_Bs(y, J, J_bar, J_h, i):
@@ -51,6 +52,7 @@ def S_experiment_Bs(y, J, J_bar, J_h, i):
 
 def dGdq2_ave_Bs(y, J, J_bar, J_h):
     # (48) of 1502.05509
+    flavio.citations.register("Descotes-Genon:2015hea")
     return (1/(1-y**2) * (observables.dGdq2(J) + observables.dGdq2(J_bar))
             - y/(1-y**2) * observables.dGdq2(J_h))/2.
 

flavio/physics/bdecays/bvll/qcdf.py

@@ -37,6 +37,8 @@ def get_input(par, B, V, scale):
 # see eqs. (18), (50), (79) of hep-ph/0106067v2
 # and eqs. (47), (48) of hep-ph/0412400
 def Cq34(q2, par, wc, B, V, scale):
+    flavio.citations.register("Beneke:2001at")
+    flavio.citations.register("Beneke:2004dp")
     # this is -C_q^12 (for q=u) or C_q^34 - eps_u * C_q^12 (for q=d,s) of hep-ph/0412400
     mB, mb, mc, alpha_s, q, eq, ed, eu, eps_u, qiqj = get_input(par, B, V, scale)
     T_t = wc['C3_'+qiqj] + 4/3.*(wc['C4_'+qiqj] + 12*wc['C5_'+qiqj] + 16*wc['C6_'+qiqj])
@@ -63,16 +65,18 @@ def T_para_minus_WA(q2, par, wc, B, V, scale):
     return -eq * 4*mB/mb * Cq34(q2, par, wc, B, V, scale)
 
 
-# B->V, weak annihaltion at O(1/mb)
+# B->V, weak annihilation at O(1/mb)
 
 # (51) of hep-ph/0412400
 
 def T_perp_WA_PowC_1(q2, par, wc, B, V, scale):
+    flavio.citations.register("Beneke:2004dp")
     mB, mb, mc, alpha_s, q, eq, ed, eu, eps_u, qiqj = get_input(par, B, V, scale)
     # NB, the remaining prefactors are added below in the function T_perp
     return -eq * 4/mb * (wc['C3_'+qiqj] + 4/3.*(wc['C4_'+qiqj] + 3*wc['C5_'+qiqj] + 4*wc['C6_'+qiqj]))
 
 def T_perp_WA_PowC_2(q2, par, wc, B, V, scale):
+    flavio.citations.register("Beneke:2004dp")
     mB, mb, mc, alpha_s, q, eq, ed, eu, eps_u, qiqj = get_input(par, B, V, scale)
     # NB, the remaining prefactors are added below in the function T_perp
     return eq * 2/mb * Cq34(q2, par, wc, B, V, scale)
@@ -84,11 +88,13 @@ def T_perp_WA_PowC_2(q2, par, wc, B, V, scale):
 
 # chromomagnetic dipole contribution
 def T_perp_plus_O8(q2, par, wc, B, V, u, scale):
+    flavio.citations.register("Beneke:2001at")
     mB, mb, mc, alpha_s, q, eq, ed, eu, eps_u, qiqj = get_input(par, B, V, scale)
     ubar = 1 - u
     return - (alpha_s/(3*pi)) * 4*ed*wc['C8eff_'+qiqj]/(u + ubar*q2/mB**2)
 
 def T_para_minus_O8(q2, par, wc, B, V, u, scale):
+    flavio.citations.register("Beneke:2001at")
     mB, mb, mc, alpha_s, q, eq, ed, eu, eps_u, qiqj = get_input(par, B, V, scale)
     ubar = 1 - u
     return (alpha_s/(3*pi)) * eq * 8 * wc['C8eff_'+qiqj]/(ubar + u*q2/mB**2)
@@ -148,6 +154,7 @@ def En_V(mB, mV, q2):
 
 # (27) of hep-ph/0106067v2
 def t_perp(q2, u, mq, par, B, V):
+    flavio.citations.register("Beneke:2001at")
     mB = par['m_'+B]
     mV = par['m_'+V]
     EV = En_V(mB, mV, q2)
@@ -162,6 +169,7 @@ def t_perp(q2, u, mq, par, B, V):
 
 # (28) of hep-ph/0106067v2
 def t_para(q2, u, mq, par, B, V):
+    flavio.citations.register("Beneke:2001at")
     mB = par['m_'+B]
     mV = par['m_'+V]
     EV = En_V(mB, mV, q2)
@@ -177,6 +185,7 @@ def B0diffBFS(q2, u, mq, mB):
 
 # (29) of hep-ph/0106067v2
 def B0(s, mq):
+    flavio.citations.register("Beneke:2001at")
     if s==0.:
         return -2.
     if 4*mq**2/s == 1.:
@@ -188,6 +197,7 @@ def B0(s, mq):
 
 # (30), (31) of hep-ph/0106067v2
 def i1_bfs(q2, u, mq, mB):
+    flavio.citations.register("Beneke:2001at")
     ubar = 1 - u
     iepsilon = 1e-8j
     mq2 = mq**2 - iepsilon
@@ -201,6 +211,7 @@ def i1_bfs(q2, u, mq, mB):
 
 # (32) of hep-ph/0106067v2
 def L1(x):
+    flavio.citations.register("Beneke:2001at")
     if x == 0.:
         return -(pi**2/6.)
     elif x == 1.:
@@ -219,12 +230,14 @@ def phiV(u, a1, a2):
 # moments of the B meson light-cone distribution amplitude as in
 # eq. (54) and surroundings of hep-ph/0106067v2
 def lB_minus(q2, par, B):
+    flavio.citations.register("Beneke:2001at")
     mB = par['m_'+B]
     mb = running.get_mb_pole(par)
     LambdaBar = mB - mb
     w0 = 2*LambdaBar/3
     return 1/(exp(-q2/mB/w0)/w0 * (-ei(q2/mB/w0) + 1j*pi))
 def lB_plus(par, B):
+    flavio.citations.register("Beneke:2001at")
     mB = par['m_'+B]
     mb = running.get_mb_pole(par)
     LambdaBar = mB - mb

flavio/physics/bdecays/bxgamma.py

@@ -40,6 +40,7 @@ def BRBXgamma(wc_obj, par, q, E0):
     cutoff $E_0$ in GeV
     (currently works only for `E0=1.6`).
     See arXiv:1503.01789 and references therein."""
+    flavio.citations.register("Misiak:2015xwa")
     scale = flavio.config['renormalization scale']['bxgamma']
     alphaem = 1/137.035999139 # this is alpha_e(0), a constant for our purposes
     bq = 'b' + q
@@ -55,6 +56,7 @@ def BRBXgamma(wc_obj, par, q, E0):
     # central value of non-perturbative correction (M. Misiak, private communication)
     P_nonpert = 0.00381745
     # eq. (2.1) of hep-ph/0609241
+    flavio.citations.register("Misiak:2006ab")
     return BRSL * abs(xi_t)**2/abs(Vcb)**2 * 6*alphaem/pi/C * (PE0 + P_nonpert + P_uncertainty)
 
 def PE0_BR_BXgamma(wc, par, q, E0):