use particle package to get PDG data

flavio/parameters.py

@@ -9,6 +9,7 @@
 import re
 from flavio.measurements import _fix_correlation_matrix
 from math import sqrt
+from particle import Particle, data as p_data
 
 
 def _read_yaml_object_metadata(obj, constraints):
@@ -91,168 +92,127 @@ def write_file(filename, constraints):
         yaml.dump(constraints.get_yaml_dict(), f)
 
 
-# particles from the PDG data file whose mass we're interested in)
-pdg_include = ['B(s)', 'B(c)', 'B(s)*', 'B*+', 'B*0', 'B+', 'B0', 'D(s)', 'D(s)*', 'D+', 'D0',
-               'H', 'J/psi(1S)', 'K(L)', 'K(S)', 'K*(892)+', 'K*(892)0', 'K+', 'K0',
-               'Lambda', 'Lambda(b)', 'Lambda(c)', 'omega(782)', 'D*(2007)', 'D*(2010)',
-               'W', 'Z',  'e', 'eta', 'f(0)(980)',
-               'mu',  'phi(1020)', 'pi+', 'pi0', 'psi(2S)', 'rho(770)+', 'rho(770)0',
-               't', 'tau', 'u', 'p', 'n']
-# dictionary translating PDG particle names to the ones in the code.
-pdg_translate = {
-    'B(s)': 'Bs',
-    'B(c)': 'Bc',
-    'D(s)': 'Ds',
-    'B(s)*': 'Bs*',
-    'D(s)*': 'Ds*',
-    'D*(2007)' : 'D*0',
-    'D*(2010)' : 'D*+',
-    'J/psi(1S)': 'J/psi',
-    'K(L)': 'KL',
-    'K(S)': 'KS',
-    'K*(892)+': 'K*+',
-    'K*(892)0': 'K*0',
-    'phi(1020)': 'phi',
-    'rho(770)0': 'rho0',
-    'rho(770)+': 'rho+',
-    'f(0)(980)': 'f0',
-    "eta'(958)": "eta'",
-    'omega(782)': 'omega',
-    'Lambda(b)' : 'Lambdab',
-    'Lambda(c)' : 'Lambdac',
-    'Higgs' : 'h', # this is necessary for the 2013 data file
-    'H' : 'h',
-    }
-
-def _read_pdg_masswidth(filename):
-    """Read the PDG mass and width table and return a dictionary.
-
-    Parameters
-    ----------
-    filname : string
-        Path to the PDG data file, e.g. 'data/pdg/mass_width_2015.mcd'
-
-    Returns
-    -------
-    particles : dict
-        A dictionary where the keys are the particle names with the charge
-        appended in case of a multiplet with different masses, e.g. 't'
-        for the top quark, 'K+' and 'K0' for kaons.
-        The value of the dictionary is again a dictionary with the following
-        keys:
-        - 'id': PDG particle ID
-        - 'mass': list with the mass, postitive and negative error in GeV
-        - 'width': list with the width, postitive and negative error in GeV
-        - 'name': same as the key
-    """
-    data = pkgutil.get_data('flavio.physics', filename)
-    lines = data.decode('utf-8').splitlines()
-    particles_by_name = {}
-    for line in lines:
-        if  line.strip()[0] == '*':
-            continue
-        mass = ((line[33:51]), (line[52:60]), (line[61:69]))
-        if mass[0].replace(' ', '') == '':
-            # if mass is empty, go to next line
-            # (necessasry for 2019 neutrino entries)
-            continue
-        mass = [float(m) for m in mass]
-        width = ((line[70:88]), (line[89:97]), (line[98:106]))
-        if  width[0].strip() == '':
-            width = (0,0,0)
+# particles from the PDG data file whose mass and width we're interested in
+pdg_particles = {
+    'Bs': 531,
+    'Bc': 541,
+    'Bs*': 533,
+    'B*+': 523,
+    'B*0': 513,
+    'B+': 521,
+    'B0': 511,
+    'Ds': 431,
+    'Ds*': 433,
+    'D+': 411,
+    'D0': 421,
+    'h': 25,
+    'J/psi': 443,
+    'KL': 130,
+    'KS': 310,
+    'K*+': 323,
+    'K*0': 313,
+    'K+': 321,
+    'K0': 311,
+    'Lambda': 3122,
+    'Lambdab': 5122,
+    'Lambdac': 4122,
+    'omega': 223,
+    'D*0': 423,
+    'D*+': 413,
+    'W': 24,
+    'Z': 23,
+    'e': 11,
+    'eta': 221,
+    'f0': 9010221,
+    'mu': 13,
+    'phi': 333,
+    'pi+': 211,
+    'pi0': 111,
+    'psi(2S)': 100443,
+    'rho+': 213,
+    'rho0': 113,
+    't': 6,
+    'tau': 15,
+    'u': 2,
+    'p': 2212,
+    'n': 2112,
+}
+
+_pdg_tex_regex = re.compile(
+    r"^([A-Za-z\\/]+)" # latin or greek letters or slash
+    r"(?:_\{(.*?)\})*" # _{...}
+    r"(?:\^\{(.*?)\})*" # ^{...}
+    r"(?:\((.*?)\))*" # (...)
+    r"(?:\^\{(.*?)\})*" # ^{...}
+)
+
+def _pdg_tex_simplify(string):
+    m = _pdg_tex_regex.match(string)
+    name = m.group(1)
+    sub = m.group(2)
+    sup = (m.group(3) or '') + (m.group(5) or '')
+    par = m.group(4)
+    if sub or name in {'W', 'Z', 'H', 'e', '\\mu', '\\tau'}:
+        # remove superscripts +-0 and keep only *
+        sup = '*' if '*' in sup else ''
+    if not sub and par and not par.isdigit() and name != 'J/\\psi':
+        # subscript absent and parantheses contain letter but not for 'J/\\psi'
+        sub = par
+    sub_tex = r'_{' + sub + r'}' if sub else ''
+    sup_tex = r'^{' + sup + r'}' if sup else ''
+    return name + sub_tex + sup_tex
+
+def _read_pdg_parameters(particle, flavio_name, tex_name, read_parameter):
+    if read_parameter == 'm':
+        name = 'm_' + flavio_name
+        tex = r'$m_{' + tex_name + '}$'
+        if tex_name =='t':
+            description = r'$' + tex_name + r'$ quark pole mass'
         else:
-            width = [float(w) for w in width]
-        ids = line[0:32].split()
-        charges = line[107:128].split()[1].split(',')
-        if len(ids) != len(charges):
-            raise ValueError()
-        for i, id_ in enumerate(ids):
-            particle = {}
-            particle_charge = charges[i].strip()
-            particle[particle_charge] = {}
-            particle[particle_charge]['id'] = id_.strip()
-            particle[particle_charge]['mass']  = mass
-            particle[particle_charge]['charge']  = particle_charge
-            particle[particle_charge]['width'] = width
-            particle_name = line[107:128].split()[0]
-            particle[particle_charge]['name'] = particle_name
-            if particle_name in particles_by_name.keys():
-                particles_by_name[particle_name].update(particle)
-            else:
-                particles_by_name[particle_name] = particle
-    result = { k + kk: vv for k, v in particles_by_name.items() for kk, vv in v.items() if len(v) > 1}
-    result.update({ k: list(v.values())[0] for k, v in particles_by_name.items() if len(v) == 1})
-    return result
-
-def _pdg_particle_string_to_tex(string):
-    regex = re.compile(r"^([A-Za-z]+)([\*+-0]*)(?:\((\d*[A-Za-z]+)\))*(?:\((\d+)\))*([\*+-0]*)$")
-    m = regex.match(string)
-    if m is None:
-        if string=='J/psi(1S)':
-            return r'J/\psi'
-        return string
-    sup = m.group(2) + m.group(5)
-    sub = m.group(3)
-    out = m.group(1)
-    if len(out) > 1:
-        out = "\\" + out
-    if sup is not None and sup != '':
-        out += r'^{' + sup + r'}'
-    if sub is not None and sub != '':
-        out += r'_{' + sub + r'}'
-    return out
+            description = r'$' + tex_name + r'$ mass'
+        central = particle.mass*1e-3
+        right = particle.mass_upper*1e-3
+        left = particle.mass_lower*1e-3
+    elif read_parameter == 'tau':
+        name = 'tau_' + flavio_name
+        tex = r'$\tau_{' + tex_name + '}$'
+        description = r'$' + tex_name + r'$ lifetime'
+        G_central = particle.width*1e-3
+        G_right = particle.width_upper*1e-3
+        G_left = particle.width_lower*1e-3
+        central = 1/G_central # life time = 1/width
+        right = G_right/G_central**2
+        left = G_left/G_central**2
+    return (name, tex, description, central, right, left)
 
 def read_pdg(year, constraints):
     """Read particle masses and widths from the PDG data file of a given year."""
-    particles = _read_pdg_masswidth('data/pdg/mass_width_' + str(year) + '.mcd')
-    for particle in pdg_include:
-        parameter_name = 'm_' + pdg_translate.get(particle, particle) # translate if necessary
-        tex_name = _pdg_particle_string_to_tex(particle)
-        try:
-            # if parameter already exists, remove existing constraints on it
-            m = Parameter[parameter_name]
-            constraints.remove_constraint(parameter_name)
-        except KeyError:
-            # otherwise, create it
-            m = Parameter(parameter_name)
-        m.tex = r'$m_{' + tex_name + '}$'
-        if tex_name =='t':
-            m.description = r'$' + tex_name + r'$ quark pole mass'
-        else:
-            m.description = r'$' + tex_name + r'$ mass'
-        m_central, m_right, m_left = particles[particle]['mass']
-        m_left = abs(m_left) # make left error positive
-        if m_right == m_left:
-            constraints.add_constraint([parameter_name],
-                NormalDistribution(m_central, m_right))
-        else:
-            constraints.add_constraint([parameter_name],
-                AsymmetricNormalDistribution(m_central,
-                right_deviation=m_right, left_deviation=m_left))
-        if particles[particle]['width'][0] == 0: # 0 is for particles where the width is unknown (e.g. B*)
-            continue
-        G_central, G_right, G_left = particles[particle]['width']
-        G_left = abs(G_left) # make left error positive
-        parameter_name = 'tau_' + pdg_translate.get(particle, particle) # translate if necessary
-        try:
-            # if parameter already exists, remove existing constraints on it
-            tau = Parameter[parameter_name]
-            constraints.remove_constraint(parameter_name)
-        except KeyError:
-            # otherwise, create it
-            tau = Parameter(parameter_name)
-        tau.tex = r'$\tau_{' + tex_name + '}$'
-        tau.description = r'$' + tex_name + r'$ lifetime'
-        tau_central = 1/G_central # life time = 1/width
-        tau_left = G_left/G_central**2
-        tau_right = G_right/G_central**2
-        if tau_left == tau_right:
-            constraints.add_constraint([parameter_name],
-                NormalDistribution(tau_central, tau_right))
-        else:
-            constraints.add_constraint([parameter_name],
-                AsymmetricNormalDistribution(tau_central,
-                    right_deviation=tau_right, left_deviation=tau_left))
+    Particle.load_table(p_data.open_text(p_data, "particle{}.csv".format(year)))
+    for flavio_name, pdgid in pdg_particles.items():
+        particle = Particle.from_pdgid(pdgid)
+        tex_name = _pdg_tex_simplify(particle.latex_name)
+        read_parameters = ['m']
+        if {particle.width, particle.width_upper, particle.width_lower}.isdisjoint({None, 0}):
+            read_parameters += ['tau']
+        for read_parameter in read_parameters:
+            data = _read_pdg_parameters(particle, flavio_name, tex_name, read_parameter)
+            (name, tex, description, central, right, left) = data
+            try:
+                # if parameter already exists, remove existing constraints on it
+                p = Parameter[name]
+                constraints.remove_constraint(name)
+            except KeyError:
+                # otherwise, create it
+                p = Parameter(name)
+            p.tex = tex
+            p.description = description
+            if right == left:
+                constraints.add_constraint([name],
+                    NormalDistribution(central, right))
+            else:
+                constraints.add_constraint([name],
+                    AsymmetricNormalDistribution(central,
+                    right_deviation=right, left_deviation=left))