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):
@@ -90,169 +91,166 @@ def write_file(filename, constraints):
     with open(filename, 'w') as f:
         yaml.dump(constraints.get_yaml_dict(), f)
 
+class FlavioParticle(Particle):
+    """This class extends the `particle.Particle` class.
 
-# 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',
+    Additional class methods
+    ------------------------
+    - from_flavio_name(flavio_name)
+      returns a class instance for a given `flavio_name`
+    - flavio_all()
+      returns a set of all class instances used in flavio
+
+    Additional properties
+    ---------------------
+    - flavio_name
+      the particle name as used in flavio if defined, otherwise `None`
+    - latex_name_simplified
+      a simplified version of the latex name returned by `latex_name`
+    - flavio_m
+      a tuple with data on the particle mass as used in flavio, containing
+      entries `name`, `tex`, `description`, `central`, `right`, `left`
+    - flavio_tau
+      a tuple with data on the particle lifetime as used in flavio, containing
+      entries `name`, `tex`, `description`, `central`, `right`, `left`
+    """
+
+    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_particles_inv = {v:k for k,v in PDG_PARTICLES.items()}
+    _pdg_tex_regex = re.compile(
+        r"^([A-Za-z\\/]+)" # latin or greek letters or slash
+        r"(?:_\{(.*?)\})*" # _{...}
+        r"(?:\^\{(.*?)\})*" # ^{...}
+        r"(?:\((.*?)\))*" # (...)
+        r"(?:\^\{(.*?)\})*" # ^{...}
+    )
 
-def _read_pdg_masswidth(filename):
-    """Read the PDG mass and width table and return a dictionary.
+    @classmethod
+    def from_flavio_name(cls, flavio_name):
+        return cls.from_pdgid(cls.PDG_PARTICLES[flavio_name])
 
-    Parameters
-    ----------
-    filname : string
-        Path to the PDG data file, e.g. 'data/pdg/mass_width_2015.mcd'
+    @classmethod
+    def flavio_all(cls):
+        return {particle for particle in cls.all() if particle.flavio_name}
 
-    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)
-        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
+    @property
+    def flavio_name(self):
+        return self._pdg_particles_inv.get(self.pdgid, None)
+
+    @property
+    def latex_name_simplified(self):
+        m = self._pdg_tex_regex.match(self.latex_name)
+        if m is None:
+            return self.latex_name
+        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 _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
+    @property
+    def flavio_m(self):
+        name = 'm_' + self.flavio_name
+        tex = r'$m_{' + self.latex_name_simplified + '}$'
+        pole_mass = ' quark pole' if self.name == 't' else ''
+        description = r'${}${} mass'.format(
+            self.latex_name_simplified, pole_mass
+        )
+        central = self.mass*1e-3
+        right = self.mass_upper*1e-3
+        left = self.mass_lower*1e-3
+        return name, tex, description, central, right, left
+
+    @property
+    def flavio_tau(self):
+        if {self.width, self.width_upper, self.width_lower} & {None, 0}:
+            return None
+        name = 'tau_' + self.flavio_name
+        tex = r'$\tau_{' + self.latex_name_simplified + '}$'
+        description = r'${}$ lifetime'.format(self.latex_name_simplified)
+        G_central = self.width*1e-3
+        G_right = self.width_upper*1e-3
+        G_left = self.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))
+    FlavioParticle.load_table(p_data.open_text(p_data, "particle{}.csv".format(year)))
+    for particle in FlavioParticle.flavio_all():
+        for data in (particle.flavio_m, particle.flavio_tau):
+            if data is None:
+                continue
+            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))