03_Scikit-HEP_particles-decays.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1><b><span style=\"color:blue\">Particles &amp; decays</span></b></h1></center>\n",
    "\n",
    "&nbsp;<br>\n",
    "**Particles and decays are central \"concepts\" in Particle Physics.** Could one really do without dedicated packages in the ecosystem? ...\n",
    "Of course not!\n",
    "\n",
    "### **Quick intro to the following packages**\n",
    "- [Particle](https://github.com/scikit-hep/particle) - *Particle Data Group* particle data, Monte Carlo identification codes, and more.\n",
    "- [DecayLanguage](https://github.com/scikit-hep/decaylanguage) - Decay files (notably for EvtGen), universal description of decay chains."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;\n",
    "<center><img src=\"images/logo_particle.png\" alt=\"Particle package logo\" style=\"width: 150px;\"/></center>\n",
    "\n",
    "<center><h2><b><span style=\"color:green\">PDG particle data, MC identification codes</span></b></h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The [Particle Data Group](https://pdg.lbl.gov/) (PDG) is an \"international collaboration that provides a comprehensive summary of Particle Physics and related areas of Cosmology: the Review of Particle Physics.\"**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"images/intro_PDG.jpg\" style=\"width:55%\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Review of Particle Physics** is a document whose importance is impossible to oversell.\n",
    "Some interesting facts about it and its more recent edition [https://pdg.lbl.gov/#about]:\n",
    "\n",
    "<center><img src=\"images/intro_PDG_about.png\" style=\"width:55%\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Package motivation - particle data\n",
    "\n",
    "- The PDG provides a series of downloadable <span style=\"color:green\">*Computer Readable Files* and in particular a table of particle masses, widths, etc. and PDG Monte Carlo particle ID numbers</span> (PDG IDs).\n",
    "The most recent file is [here](http://pdg.lbl.gov/2021/html/computer_read.html).\n",
    "- It <span style=\"color:green\">also provided an experimental file with extended information</span>\n",
    "(spin, quark content, P and C parities, etc.) until 2008 only, see [here](http://pdg.lbl.gov/2008/html/computer_read.html) (not widely known!).\n",
    "\n",
    "- But <span style=\"color:green\"><i>anyone</i> wanting to use these data</span>, the only readily available,\n",
    "<span style=\"color:green\">has to parse the file programmatically</span>.\n",
    "- Why not make a Python package to deal with all these data, for everyone?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Package motivation - MC identification codes\n",
    "\n",
    "- The <span style=\"color:green\">C++ HepPID and HepPDT libraries provide functions for processing particle ID codes</apan>\n",
    "in the standard particle (aka PDG) numbering scheme.\n",
    "- Different event generators may have their separate set of particle IDs: Geant3, etc.\n",
    "- Again, why not make a package providing all functionality/conversions, Python-ically, for everyone?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Pythonic interface to**\n",
    "- Particle Data Group (PDG) particle data table.\n",
    "- Particle MC identification codes, with inter-MC converters.\n",
    "- With various extra goodies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Package, in short\n",
    "\n",
    "- <span style=\"color:green\">Particle</span> - loads extended <b>PDG data tables</b> and implements search and manipulations / display.\n",
    "- <span style=\"color:green\">PDGID</span> - find out as much as possible from the PDG ID number. <b>No table lookup</b>.\n",
    "- <span style=\"color:green\">Converters for MC IDs</span> used in [Pythia](https://pythia.org/) and Geant3.\n",
    "- Flexible / advanced usage programmatically.\n",
    "- Basic usage via the command line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **1. `PDGID` class and MC ID classes**\n",
    "\n",
    "\n",
    "- Classes `PDGID`, `PythiaID`, `Geant3ID`.\n",
    "- Converters in module `particle.converters`: `Geant2PDGIDBiMap`, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### PDG IDs module overview\n",
    "\n",
    "- <span style=\"color:green\">Process and query PDG IDs</span>, and more – no look-up table needed.\n",
    "  - Current version of package reflects the latest version of the\n",
    "    <span style=\"color:green\">HepPID & HepPDT utility functions</span> defined in the C++ HepPID and HepPDT versions 3.04.01\n",
    "  - It contains more functionality than that available in the C++ code … and minor fixes too.\n",
    "- Definition of a <span style=\"color:green\">PDGID class, PDG ID literals</span>,\n",
    "and set of standalone HepPID <span style=\"color:green\">functions to query PDG IDs</span>\n",
    "(is_meson, has_bottom, j_spin, charge, etc.).\n",
    "   - All PDGID class functions are available standalone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### PDGID class\n",
    "- Wrapper class `PDGID` for PDG IDs.\n",
    "- Behaves like an int, with extra goodies.\n",
    "- Large spectrum of properties and methods, with a Pythonic interface, and yet more!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from particle import PDGID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pid = PDGID(211)\n",
    "pid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDGID(99999999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from particle.pdgid import is_meson\n",
    "\n",
    "pid.is_meson, is_meson(pid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To print all `PDGID` properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pid.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### MC ID classes and converters\n",
    "\n",
    "- <span style=\"color:green\">Classes for MC IDs</span> used in Pythia and Geant3: `PythiaID` and `Geant3ID`.\n",
    "- <span style=\"color:green\">ID converters</span> in module `particle.converters`: `Geant2PDGIDBiMap`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle import PythiaID, Geant3ID\n",
    "\n",
    "pyid = PythiaID(10221)\n",
    "\n",
    "pyid.to_pdgid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Conversions are directly available via mapping classes.\n",
    "\n",
    "E.g., bi-directional map Pythia ID - PDG ID:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle.converters import Pythia2PDGIDBiMap\n",
    "\n",
    "Pythia2PDGIDBiMap[PDGID(9010221)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pythia2PDGIDBiMap[PythiaID(10221)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **2. `Particle` class**\n",
    "\n",
    "There are various ways to create a particle. The often used method is via its PDG ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle import Particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Particle.from_pdgid(211)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Searching\n",
    "\n",
    "<span style=\"color:green\">Simple and natural API</span> to deal with the PDG particle data table,<br>with <span style=\"color:green\">powerful 1-line search and look-up utilities!</span>\n",
    "\n",
    "- `Particle.find(…)` – search a single match (exception raised if multiple particles match the search specifications).\n",
    "- `Particle.findall(…)` – search a list of candidates.\n",
    "\n",
    "- Search methods that can query any particle property!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Particle.find('J/psi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can specify search terms as keywords - _any particle property_:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can directly check the numeric charge:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Particle.findall('pi', charge=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Or use a **lambda function** for the ultimate in generality! For example, to find all the neutral particles with a bottom quark between 5.2 and 5.3 GeV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepunits import GeV, s  # Units are good. Use them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Particle.findall(lambda p:\n",
    "                     p.pdgid.has_bottom\n",
    "                     and p.charge==0\n",
    "                     and 5.2*GeV < p.mass < 5.3*GeV\n",
    "                )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Another lambda function example: You can use the width or the lifetime:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Particle.findall(lambda p: p.lifetime > 1000*s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If you want infinite lifetime, you could just use the keyword search instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Particle.findall(lifetime=float('inf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Trivially find all pseudoscalar charm mesons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle import SpinType\n",
    "\n",
    "Particle.findall(lambda p: p.pdgid.is_meson and p.pdgid.has_charm and p.spin_type==SpinType.PseudoScalar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Display\n",
    "\n",
    "Nice display in Jupyter notebooks, as well as `str` and `repr` support:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = Particle.from_pdgid(-415)\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(repr(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Full descriptions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You may find LaTeX or HTML to be more useful in your program; both are supported:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p.latex_name, '\\n', p.html_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Particle properties\n",
    "\n",
    "You can do things to particles, like **invert** them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "~p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are a plethora of properties you can access:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.spin_type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can quickly access the PDGID of a particle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.pdgid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **3. Literals**\n",
    "\n",
    "They provide a <span style=\"color:green\">handy way to manipulate things with human-readable names!</span>\n",
    "\n",
    "`Particle` defines <span style=\"color:green\">literals for all particles</span>, with easily recognisable (programmatic friendly) names.\n",
    "- Literals are dynamically generated on import for both `PDGID` and `Particle` classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Particle literals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle import literals as lp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lp.phi_1020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### PDGID literals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from particle.pdgid import literals as lid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lid.phi_1020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **4. Data files, stored in `particle/data/`**\n",
    "\n",
    "- <b>PDG particle data files</b>\n",
    "  - Original PDG data files, which are in a fixed-width format - simply for bookkeeping and reference.\n",
    "  - Code rather uses “digested forms” of these, produced within `Particle`, stored as CSV, for optimised querying.\n",
    "  - Latest PDG data (2020) used by default.\n",
    "  - Advanced usage: user can load older PDG tables, load a “user table” with new particles, append to default table.\n",
    "\n",
    "- <b>Other data files</b>\n",
    "  - CSV file for mapping of PDG IDs to particle LaTeX names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Dump table contents**\n",
    "\n",
    "The package provides the 2 methods `Particle.to_list(...)` and `Particle.to_dict(...)`, which make it easy to dump (selected) particle properties in an easy way. No need to dig into the package installation directory to inspect the particle data table ;-).\n",
    "\n",
    "Tabular output can be formatted with the powerful package `tabulate`, for example (other similar libraries exist)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tabulate import tabulate\n",
    "\n",
    "query_as_list = Particle.to_list(filter_fn=lambda p: p.pdgid.is_lepton and p.charge!=0, exclusive_fields=['pdgid', 'name', 'mass', 'charge'])\n",
    "\n",
    "print(tabulate(query_as_list, headers=\"firstrow\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fancy creating tables of particle properties in, say, HTML or reStructuredText format, as below? Then check out the *exercises/particle.ipynb* notebook in the exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table>\n",
    "<tr style=\"background: white;\">\n",
    "    <td align=\"center\"><img src=\"images/Scikit-HEP_gallery_Particle.jpg\" width=\"80%\"></td>\n",
    "    <td align=\"center\"><img src=\"images/Scikit-HEP_gallery_Particle_ex-table-rst.png\" width=\"80%\"></td>\n",
    "</tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **5. Advanced usage**\n",
    "\n",
    "You can:\n",
    "\n",
    "* Extend or replace the default particle data table in `Particle`.\n",
    "* Adjust properties for a particle.\n",
    "* Make custom particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;<br>\n",
    "<center>\n",
    "    <img src=\"images/logo_decaylanguage.png\" style=\"width: 200px;\"/>\n",
    "    <h2><b><span style=\"color:green\">Decay files, universal description of decay chains</span></b></h2>\n",
    "</center>\n",
    "\n",
    "`DecayLanguage` is designed for the manipulation of decay structures in Python. The current package has:\n",
    "\n",
    "- Decay file parsers:\n",
    "  - Read *.dec DecFiles*, such as [EvtGen](https://evtgen.hepforge.org/)<sup>(*)</sup> decay files typically used in Flavour Physics experiments.\n",
    "  - Manipulate and visualise them in Python.\n",
    "- Amplitude Analysis decay language:\n",
    "  - Input based on AmpGen generator, output format for GooFit C++ program.\n",
    "  \n",
    "> <sup>(*)</sup> *\"EvtGen is a Monte Carlo event generator that simulates the decays of heavy flavour particles, primarily B and D mesons. It contains a range of decay models for intermediate and final states containing scalar, vector and tensor mesons or resonances, as well as leptons, photons and baryons. Decay amplitudes are used to generate each branch of a given full decay tree, taking into account angular and time-dependent correlations which allows for the simulation of CP-violating processes...\"*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Package motivation\n",
    "\n",
    "- Ability to describe decay-tree-like structures.\n",
    "- Provide a translation of decay amplitude models from AmpGen to GooFit.\n",
    "  - Idea is to generalise this to other decay descriptions.\n",
    "- Any experiment uses event generators which, among many things, need to describe particle decay chains.\n",
    "- Programs such as EvtGen rely on so-called .dec decay files.\n",
    "- Many experiments need decay data files.\n",
    "- Why not make a Python package to deal with decay files, for everyone?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Package, in short\n",
    "\n",
    "- Tools to parse decay files and programmatically manipulate them, query, display information.\n",
    "  - Descriptions and parsing built atop the [Lark parser](https://github.com/lark-parser/lark/).\n",
    "- Tools to translate decay amplitude models from AmpGen to GooFit, and manipulate them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **1. Decay files**\n",
    "\n",
    "#### *Master file” DECAY.DEC\n",
    "\n",
    "<span style=\"color:green\">Gigantic file defining decay modes for all relevant particles, including decay model specifications.</span>\n",
    "The LHCb experiment uses one. Belle II as well, and others.\n",
    "\n",
    "#### User .dec files\n",
    "- Needed to produce specific MC samples.\n",
    "- Typically contain a single decay chain (except if defining inclusive samples)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Example user decay file:**\n",
    "\n",
    "<small>\n",
    "<pre>\n",
    "# Decay file for [B_c+ -> (B_s0 -> K+ K-) pi+]cc\n",
    "\n",
    "Alias      B_c+sig        B_c+\n",
    "Alias      B_c-sig        B_c-\n",
    "ChargeConj B_c+sig        B_c-sig\n",
    "Alias      MyB_s0         B_s0\n",
    "Alias      Myanti-B_s0    anti-B_s0\n",
    "ChargeConj MyB_s0         Myanti-B_s0\n",
    "\n",
    "Decay B_c+sig\n",
    "  1.000     MyB_s0     pi+     PHOTOS PHSP;\n",
    "Enddecay\n",
    "CDecay B_c-sig\n",
    "\n",
    "Decay MyB_s0\n",
    "    1.000     K+     K-     SSD_CP 20.e12 0.1 1.0 0.04 9.6 -0.8 8.4 -0.6;\n",
    "Enddecay\n",
    "CDecay Myanti-B_s0\n",
    "</pre>\n",
    "</small>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **2. Decay file parsing**\n",
    "\n",
    "- Parsing should be simple\n",
    "- Parsing should be (reasonably) fast!\n",
    "\n",
    "After parsing, many queries are possible!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from decaylanguage import DecFileParser"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The LHCb \"master decay file\"\n",
    "\n",
    "It's a big file! ~ 500 particle decays defined, thousands of decay modes, over 11k lines in total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp = DecFileParser('data/DECAY_LHCB.DEC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "dfp.parse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's parse and play with a small decay file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('data/Dst.dec') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dfp_Dst = DecFileParser('data/Dst.dec')\n",
    "dfp_Dst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp_Dst.parse()\n",
    "dfp_Dst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It can be handy to **parse from a multi-line string** rather than a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = \"\"\"\n",
    "# Decay file for [B_c+ -> (B_s0 -> K+ K-) pi+]cc\n",
    "\n",
    "Alias      B_c+sig        B_c+\n",
    "Alias      B_c-sig        B_c-\n",
    "ChargeConj B_c+sig        B_c-sig\n",
    "Alias      MyB_s0         B_s0\n",
    "Alias      Myanti-B_s0    anti-B_s0\n",
    "ChargeConj MyB_s0         Myanti-B_s0\n",
    "\n",
    "Decay B_c+sig\n",
    "  1.000     MyB_s0     pi+     PHOTOS PHSP;\n",
    "Enddecay\n",
    "CDecay B_c-sig\n",
    "\n",
    "Decay MyB_s0\n",
    "    1.000     K+     K-     SSD_CP 20.e12 0.1 1.0 0.04 9.6 -0.8 8.4 -0.6;\n",
    "Enddecay\n",
    "CDecay Myanti-B_s0\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dfp = DecFileParser.from_string(s)\n",
    "dfp.parse()\n",
    "dfp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Decay file information\n",
    "\n",
    "Several methods are available, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp_Dst.list_decay_mother_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp_Dst.print_decay_modes('D*+')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Info such as particle aliases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp.dict_aliases()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfp.dict_charge_conjugates()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **3.  Display of decay chains**\n",
    "\n",
    "The parser can provide a simple `dict` representation of any decay chain found in the input decay file(s). Being generic and simple, that is what is used as input information for the viewer class (see below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc = dfp_Dst.build_decay_chains('D+')\n",
    "dc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from decaylanguage import DecayChainViewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "DecayChainViewer(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dc = dfp_Dst.build_decay_chains('D*+', stable_particles=['D+', 'D0', 'pi0'])\n",
    "DecayChainViewer(dc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### **Charge conjugation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_cc = dfp_Dst.build_decay_chains('D*-', stable_particles=['D-', 'anti-D0', 'pi0'])\n",
    "DecayChainViewer(dc_cc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### **Parsing several files**\n",
    "\n",
    "Typically useful when the user decay file needs information from the master decay file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = u\"\"\"\n",
    "Alias      MyXic+              Xi_c+\n",
    "Alias      MyantiXic-          anti-Xi_c-\n",
    "ChargeConj MyXic+              MyantiXic-\n",
    "\n",
    "Decay Xi_cc+sig\n",
    "  1.000       MyXic+    pi-    pi+       PHSP;\n",
    "Enddecay\n",
    "CDecay anti-Xi_cc-sig\n",
    "\n",
    "Decay MyXic+\n",
    "  1.000       p+    K-    pi+       PHSP;\n",
    "Enddecay\n",
    "CDecay MyantiXic-\n",
    "\n",
    "End\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dfp = DecFileParser.from_string(s)\n",
    "dfp.parse()\n",
    "dfp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the subtletly: 3, not 4 decays, are found! This is because the file contains no statement\n",
    "`ChargeConj anti-Xi_cc-sigXi_cc+sig`, hence the parser cannot know to which particle (matching `Decay` statement) the charge-conjugate decay of `anti-Xi_cc-sig` relates to (code does not rely on position of statements to guess ;-))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "d = dfp.build_decay_chains('Xi_cc+sig')\n",
    "DecayChainViewer(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As said in the warning, the information provided is not enough for the anti-Xi_cc-sig to make sense:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from decaylanguage.dec.dec import DecayNotFound\n",
    "\n",
    "try:\n",
    "    d = dfp.build_decay_chains('anti-Xi_cc-sig')\n",
    "except DecayNotFound:\n",
    "    print(\"Decays of particle 'anti-Xi_cc-sig' not found in .dec file!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "But the missing information is easily providing **parsing two files simultaneously ...!** (Any number of files is allowed.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from tempfile import NamedTemporaryFile\n",
    "\n",
    "with NamedTemporaryFile(delete=False) as tf:\n",
    "    tf.write(s.encode('utf-8'))\n",
    "    \n",
    "dfp = DecFileParser(tf.name, 'data/DECAY_LHCB.DEC')\n",
    "dfp.parse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dc = dfp.build_decay_chains('Xi_cc+sig')\n",
    "\n",
    "DecayChainViewer(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dc_cc = dfp.build_decay_chains('anti-Xi_cc-sig')\n",
    "\n",
    "DecayChainViewer(dc_cc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Want to save a graph? Try for example \n",
    "```python\n",
    "dcv = DecayChainViewer(...)\n",
    "dcv.graph.save(...)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **4. Representation of decay chains**\n",
    "\n",
    "<span style=\"color:green\">The universal (and digital) representation of decay chains is of interest well outside the context of decay file parsing!</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Building blocks\n",
    "\n",
    "- A <span style=\"color:green\">daughters list</span> - list of final-state particles.\n",
    "- A <span style=\"color:green\">decay mode</span> - typically a branching fraction and a list of final-state particles (may also contain _any_ metadata such as decay model and optional decay-model parameters, as defined for example in .dec decay files).\n",
    "- A <span style=\"color:green\">decay chain</span> - can be seen as a mother particle and a list of decay modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from decaylanguage.decay.decay import DaughtersDict, DecayMode, DecayChain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Daughters list** (actually a ``Counter`` dictionary, internally):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Constructor from a dictionary\n",
    "dd = DaughtersDict({'K+': 1, 'K-': 2, 'pi+': 1, 'pi0': 1})\n",
    "\n",
    "# Constructor from a list of particle names\n",
    "dd = DaughtersDict(['K+', 'K-', 'K-', 'pi+', 'pi0'])\n",
    "\n",
    "# Constructor from a string representing the final state\n",
    "dd = DaughtersDict('K+ K- pi0')\n",
    "dd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decay Modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# A 'default' and hence empty, decay mode\n",
    "dm = DecayMode()\n",
    "\n",
    "# Decay mode with minimal input information\n",
    "dd = DaughtersDict('K+ K-')\n",
    "dm = DecayMode(0.5, dd)\n",
    "\n",
    "# Decay mode with decay model information and user metadata\n",
    "dm = DecayMode(0.2551,                                              # branching fraction\n",
    "               'pi- pi0 nu_tau',                                    # final-state particles\n",
    "               model='TAUHADNU',                                    # decay model\n",
    "               model_params=[-0.108, 0.775, 0.149, 1.364, 0.400],   # decay-model parameters\n",
    "               study='toy', year=2019                               # user metadata\n",
    "              )\n",
    "dm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(dm.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Various manipulations are available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = DecayMode.from_pdgids(0.5, [321, -321])\n",
    "print(dm)\n",
    "\n",
    "dm = DecayMode(1.0, 'K+ K+ pi-')\n",
    "dm.charge_conjugate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Decay chains "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm1 = DecayMode(0.0124, 'K_S0 pi0', model='PHSP')\n",
    "dm2 = DecayMode(0.692, 'pi+ pi-')\n",
    "dm3 = DecayMode(0.98823, 'gamma gamma')\n",
    "dc = DecayChain('D0', {'D0':dm1, 'K_S0':dm2, 'pi0':dm3})\n",
    "\n",
    "dc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc.decays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Flatten the decay chain, i.e. replace all intermediate, decaying particles, with their final states:\n",
    "- The BF is now the *visible BF*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dc.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course you can sill just as easily visualise decays defined via this `DecayChain` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "DecayChainViewer(dc.to_dict())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}