Properly deal with CP conjugate decay definitions, even with aliases

decaylanguage/dec/dec.py

@@ -6,7 +6,7 @@
 from lark import Lark
 from lark import Tree, Transformer, Visitor
 
-from particle import Particle
+from particle import Particle, ParticleNotFound
 
 from ..data import open_text
 from .. import data
@@ -100,6 +100,7 @@ def parse(self, include_cdecays=True):
         include_cdecays: boolean, optional, default=True
             Choose whether or not to consider charge-conjugate decays,
             which are specified via "CDecay <MOTHER>".
+            Make sure you understand the consequences of ignoring CP conj. decays!
         """
         # Has a file been parsed already?
         if self._parsed_decays is not None:
@@ -120,9 +121,13 @@ def parse(self, include_cdecays=True):
         dec_file = open(self._dec_file_name).read()
         self._parsed_dec_file = parser.parse(dec_file)
 
-        # At last, find all particle decays defined in the .dec decay file
+        # At last, find all particle decays defined in the .dec decay file ...
         self._find_parsed_decays()
 
+        # ... and create on the fly the CP conjugate decays, if requested
+        if self._include_cdecays:
+            self._add_charge_conjugate_decays()
+
     def grammar(self):
         """
         Access the internal Lark grammar definition file,
@@ -270,18 +275,66 @@ def _find_parsed_decays(self):
         (Tree(decay, [Tree(particle, [Token(LABEL, <MOTHER1>]), ...),
         Tree(decay, [Tree(particle, [Token(LABEL, <MOTHER2>]), ...)).
 
+        Duplicate definitions (a bug, of course) are removed, issuing a warning.
+
         Note
         ----
         1) Method not meant to be used directly!
         2) CP conjugates need to be dealt with differently,
-        see 'list_charge_conjugate_decays()'.
+        see 'self._add_charge_conjugate_decays()'.
         """
         if self._parsed_dec_file is not None:
             self._parsed_decays = get_decays(self._parsed_dec_file)
 
         # Check for duplicates - should be considered a bug in the .dec file!
         self._check_parsed_decays()
 
+    def _add_charge_conjugate_decays(self):
+        """
+        If requested (see the 'self._include_cdecays' class attribute),
+        create the Lark Tree instances describing the CP conjugate decays
+        specified in the input parsed file via the statements of the form
+        "CDecay <MOTHER>".
+        These are added to the internal list of decays stored in the class
+        in variable 'self._parsed_decays', performing a CP transformation
+        on each CP-related decay, which is cloned.
+
+        Note
+        ----
+        Method not meant to be used directly!
+        """
+        # Cross-check - make sure CP conjugate decays are not defined
+        # with both 'Decay' and 'CDecay' statements!
+        mother_names_decays = [get_decay_mother_name(tree)
+                               for tree in self._parsed_decays]
+        mother_names_cdecays = self.list_charge_conjugate_decays()
+        duplicates = [n for n in mother_names_cdecays if n in mother_names_decays]
+        if len(duplicates) > 0:
+            msg = """The following particles are defined in the input .dec file with both 'Decay' and 'CDecay': {0}!
+The 'CDecay' definition(s) will be ignored ...""".format(', '.join(d for d in duplicates))
+            warnings.warn(msg)
+
+        # If that's the case, proceed using the decay definitions specified
+        # via the 'Decay' statement, hence discard/remove the definition
+        # via the 'CDecay' statement.
+        for d in duplicates:
+            mother_names_cdecays.remove(d)
+
+        # At last, create the CP conjugate decays:
+        # First, make a (deep) copy of the list of Tree instances
+        # describing the parsed decays.
+        # By construction, there are no CP conjugate decays in there.
+        cdecays = [ tree.__deepcopy__(None) for tree in self._parsed_decays]
+
+        # Take care of CP conjugate decays defined via aliases,
+        # passing them as CP conjugates to be processed manually ...
+        dict_cdecay_names = self.dict_charge_conjugates()
+
+        # Finally, perform all particle -> CP(particle) replacements in one go!
+        [CPConjugateReplacement(charge_conj_defs=dict_cdecay_names).visit(t)
+        for t in cdecays]
+
+
     def _check_parsing(self):
         """Has the .parse() method been called already?"""
         if self._parsed_dec_file is None:
@@ -291,35 +344,47 @@ def _check_parsed_decays(self):
         """
         Is the number of decays parsed consistent with the number of
         decay mother names? An inconsistency can arise if decays are redefined.
+
+        Duplicates are removed, starting from the second occurrence.
         """
+        # Issue a helpful warning if duplicates are found
         lmn = self.list_decay_mother_names()
+        duplicates = []
         if self.number_of_decays != len(set(lmn)):
-            warnings.warn("Input .dec file redefines decays for particle(s) {0}!".format(set([n for n in lmn if lmn.count(n)>1])))
+            duplicates = set([n for n in lmn if lmn.count(n)>1])
+            msg = """The following particle(s) is(are) redefined in the input .dec file with 'Decay': {0}!
+All but the first occurence will be discarded/removed ...""".format(', '.join(d for d in duplicates))
+            warnings.warn(msg)
+
+        # Create a list with all occurrences to remove
+        # (duplications means multiple instances to remove)
+        duplicates_to_remove = []
+        for item in duplicates:
+            c = lmn.count(item)
+            if c>1:
+                duplicates_to_remove.extend([item]*(c-1))
+
+        # Actually remove all but the first occurence of duplicate decays
+        for tree in reversed(self._parsed_decays):
+            val = tree.children[0].children[0].value
+            if val in duplicates_to_remove:
+                duplicates_to_remove.remove(val)
+                self._parsed_decays.remove(tree)
 
     @property
     def number_of_decays(self):
         """Return the number of particle decays defined in the parsed .dec file."""
         self._check_parsing()
 
-        n = len(self._parsed_decays)
-
-        if self._include_cdecays:
-            n += len(self.list_charge_conjugate_decays())
-
-        return n
+        return len(self._parsed_decays)
 
     def list_decay_mother_names(self):
         """
         Return a list of all decay mother names found in the parsed decay file.
         """
         self._check_parsing()
 
-        names = [get_decay_mother_name(d) for d in self._parsed_decays]
-
-        if self._include_cdecays:
-            names += self.list_charge_conjugate_decays()
-
-        return names
+        return [get_decay_mother_name(d) for d in self._parsed_decays]
 
     def _find_decay_modes(self, mother):
         """
@@ -476,6 +541,12 @@ class CPConjugateReplacement(Visitor):
     (search done via the Particle class in the particle package),
     its CP conjugate name is denoted as 'CPConj(UNKOWN)'.
 
+    Parameters
+    ----------
+    charge_conj_defs: dict, optional, default={}
+        Dictionary with the charge conjugate particle definitions
+        in the parsed file. Argument to be passed to the class constructor.
+
     Examples
     --------
     >>> from lark import Tree, Token
@@ -489,13 +560,19 @@ class CPConjugateReplacement(Visitor):
     [Tree(value, [Token(SIGNED_NUMBER, '1.0')]), Tree(particle, [Token(LABEL, 'K+')]),
     Tree(particle, [Token(LABEL, 'pi-')]), Tree(model, [Token(MODEL_NAME, 'PHSP')])])])
     """
+    def __init__(self, charge_conj_defs=dict()):
+        self.charge_conj_defs = charge_conj_defs
+
     # Method for the rule (here, a replacement) we wish to implement
     def particle(self, tree):
         assert tree.data == 'particle'
-        try:
-            tree.children[0].value = Particle.from_dec(tree.children[0].value).invert().name
-        except:
-            tree.children[0].value = 'CPConj({0})'.format(val)
+        if tree.children[0].value in self.charge_conj_defs:
+            tree.children[0].value = self.charge_conj_defs[tree.children[0].value]
+        else:
+            try:
+                tree.children[0].value = Particle.from_dec(tree.children[0].value).invert().name
+            except ParticleNotFound:
+                tree.children[0].value = 'CPConj({0})'.format(tree.children[0].value)
 
 
 def get_decay_mother_name(decay_tree):

tests/data/duplicate-decays.dec

@@ -0,0 +1,40 @@
+# File for testing purposes.
+# Example bug taken from an old DECAY.DEC file ;-) :
+# a particle decay is duplicated and the CP conjugate definition with
+# 'CDecay anti-Sigma(1775)0' is effectively a different way of defining what is
+# already defined via 'Decay anti-Sigma(1775)0' !
+
+Decay Sigma(1775)0
+0.215  p+             K-              PHSP;
+0.215  n0             anti-K0         PHSP;
+0.20   Lambda0        pi0             PHSP;
+0.02   Sigma+         pi-             PHSP;
+0.02   Sigma-         pi+             PHSP;
+0.055  Sigma*+        pi-             PHSP;
+0.055  Sigma*-        pi+             PHSP;
+0.22   Lambda(1520)0  pi0             PHSP;
+Enddecay
+
+Decay anti-Sigma(1775)0
+0.215  anti-p-             K+         PHSP;
+0.215  anti-n0             K0         PHSP;
+0.20   anti-Lambda0        pi0        PHSP;
+0.02   anti-Sigma+         pi-        PHSP;
+0.02   anti-Sigma-         pi+        PHSP;
+0.055  anti-Sigma*+        pi-        PHSP;
+0.055  anti-Sigma*-        pi+        PHSP;
+0.22   anti-Lambda(1520)0  pi0        PHSP;
+Enddecay
+
+Decay Sigma(1775)0    # PDG 3216
+  0.23   Lambda0 pi0 PHSP;
+  0.23   Lambda(1405)0 pi0 PHSP;
+  0.2    p+       K-        PHSP;
+  0.1    n0       K_S0      PHSP;
+  0.1    n0       K_L0      PHSP;
+  0.05   Sigma*+  pi-       PHSP;
+  0.05   Sigma*-  pi+       PHSP;
+  0.02   Sigma+   pi-       PHSP;
+  0.02   Sigma-   pi+       PHSP;
+Enddecay
+CDecay anti-Sigma(1775)0

tests/test_dec.py

@@ -5,9 +5,15 @@
 except ImportError:
     from pathlib import Path
 
+from lark import Tree, Token
+
 from decaylanguage import data
+
 from decaylanguage.dec.dec import DecFileParser
 from decaylanguage.dec.dec import DecFileNotParsed, DecayNotFound
+from decaylanguage.dec.dec import CPConjugateReplacement
+from decaylanguage.dec.dec import get_decay_mother_name
+from decaylanguage.dec.dec import get_final_state_particle_names
 
 # New in Python 3
 try:
@@ -150,9 +156,32 @@ def test_decay_mode_details():
     assert p._decay_mode_details(tree_Dp) == output
 
 
+def test_duplicate_decay_definitions():
+    p = DecFileParser(DIR / 'data/duplicate-decays.dec')
+    p.parse()
+
+    assert p.number_of_decays == 2
+
+    assert p.list_decay_mother_names() == ['Sigma(1775)0', 'anti-Sigma(1775)0']
+
+
 def test_build_decay_chain():
     p = DecFileParser.from_file(DIR / 'data/test_example_Dst.dec')
     p.parse()
 
     output = {'D+': [{'bf': 1.0, 'fs': ['K-', 'pi+', 'pi+', 'pi0'], 'm': 'PHSP', 'mp': ''}]}
     assert p.build_decay_chain('D+', stable_particles=['pi0']) == output
+
+
+def test_Lark_CPConjugateReplacement_Visitor():
+    t = Tree('decay', [Tree('particle', [Token('LABEL', 'D0')]),
+            Tree('decayline', [Tree('value', [Token('SIGNED_NUMBER', '1.0')]),
+            Tree('particle', [Token('LABEL', 'K-')]),
+            Tree('particle', [Token('LABEL', 'pi+')]),
+            Tree('model', [Token('MODEL_NAME', 'PHSP')])])])
+
+    CPConjugateReplacement().visit(t)
+
+    assert get_decay_mother_name(t) == 'D~0'
+
+    assert get_final_state_particle_names(t.children[1]) == ['K+', 'pi-']