{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%% md\n"
}
},
"source": [
"# PhaseSpace + DecayLanguage - lightning talk\n",
"Primary author: Simon Thor\n",
"\n",
"Co-authors: Jonas Eschle, Eduardo Rodrigues, Albert Puig\n",
"\n",
"This tutorial/lightning talk is based on the one in the [PhaseSpace documentation](https://github.com/zfit/phasespace/blob/master/docs/GenMultiDecay_Tutorial.ipynb).\n",
"\n",
"PhaseSpace and more specifically the new feature I have developed will:\n",
"- Make it easy to simulate phase space generation\n",
"- Load decays from DecayLanguage to PhaseSpace\n",
"- easily specify resonances at various levels of customizability \n",
"\n",
"## Installation\n",
"In order to use this functionality, you need to install the extra dependencies to PhaseSpace by running\n",
"`pip install \"phasespace[fromdecay]\"`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"# Import libraries\n",
"from pprint import pprint\n",
"from copy import deepcopy\n",
"import zfit\n",
"from particle import Particle\n",
"from decaylanguage import DecFileParser, DecayChainViewer, DecayChain, DecayMode\n",
"import tensorflow as tf\n",
"\n",
"from phasespace.fromdecay import GenMultiDecay\n",
"from phasespace import GenParticle"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Quick introduction to GenParticle\n",
"The PhaseSpace GenParticle class can be used to simulate decays of particles:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"pion = GenParticle('pi-', 139.57018)\n",
"kaon = GenParticle('K+', 493.677)\n",
"kstar = GenParticle('K*', 895.81).set_children(pion, kaon)\n",
"gamma = GenParticle('gamma', 0)\n",
"bz = GenParticle('B0', 5279.58).set_children(kstar, gamma)\n",
"weights, particles = bz.generate(n_events=1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, one cannot simulate a particle which can decay in multiple ways using GenParticle. This is instead done using the GenMultiDecay class."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parse a decay\n",
"The last presentation introduced DecayLanguage. We will now use it to parse a .dec file and simulate the decay."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n\n\n\n",
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"parser = DecFileParser('example_decays.dec')\n",
"parser.parse()\n",
"pi0_chain = parser.build_decay_chains(\"pi0\")\n",
"DecayChainViewer(pi0_chain)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating a GenMultiDecay object\n",
"A regular `phasespace.GenParticle` instance would not be able to simulate this decay, since the $\\pi^0$ particle can decay in four different ways. However, a `GenMultiDecay` object can be created directly from a DecayLanguage dict:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"pi0_decay = GenMultiDecay.from_dict(pi0_chain)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When creating a `GenMultiDecay` object, the DecayLanguage dict is \"unpacked\" into separate GenParticle instances, where each GenParticle instance corresponds to one way that the particle can decay.\n",
"\n",
"These GenParticle instances and the probabilities of that decay mode can be accessed via `GenMultiDecay.gen_particles`. This is a list of tuples, where the first element in the tuple is the probability and the second element is the GenParticle."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"There is a probability of 0.988228297 that pi0 decays into gamma, gamma [0]\n",
"There is a probability of 0.011738247 that pi0 decays into e+, e-, gamma [1]\n",
"There is a probability of 3.3392e-05 that pi0 decays into e+ [0], e+ [1], e- [0], e- [1]\n",
"There is a probability of 6.5e-08 that pi0 decays into e+ [2], e- [2]\n"
]
}
],
"source": [
"for probability, particle in pi0_decay.gen_particles:\n",
" print(f\"There is a probability of {probability} \"\n",
" f\"that pi0 decays into {', '.join(child.name for child in particle.children)}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One can simulate this decay using the `.generate` method, which works the same as the `GenParticle.generate` method.\n",
"\n",
"When calling the `GenMultiDecay.generate` method, it internally calls the generate method on the of the GenParticle instances in `GenMultiDecay.gen_particles`. The outputs are placed in a list, which is returned."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of events for each decay mode: 9881, 119\n"
]
}
],
"source": [
"weights, events = pi0_decay.generate(n_events=10_000)\n",
"print(\"Number of events for each decay mode:\", \", \".join(str(len(w)) for w in weights))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can confirm that the counts above are close to the expected counts based on the probabilities."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Changing mass settings\n",
"Since DecayLanguage dicts do not contain any information about the mass of a particle, the `fromdecay` submodule uses the [particle](https://github.com/scikit-hep/particle) package to find the mass of a particle based on its name.\n",
"The mass can either be a constant value or a function (besides the top particle, which is always a constant).\n",
"These settings can be modified by passing in additional parameters to `GenMultiDecay.from_dict`.\n",
"There are two optional parameters that can be passed to `GenMultiDecay.from_dict`: `tolerance` and `mass_converter`.\n",
"\n",
"### Constant vs variable mass\n",
"If a particle has a width less than `tolerance`, its mass is set to a constant value.\n",
"This will be demonsttrated with the decay below:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n\n\n\n",
"text/plain": [
""
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dsplus_chain = parser.build_decay_chains(\"D*+\", stable_particles=[\"D+\"])\n",
"DecayChainViewer(dsplus_chain)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pi0 width = 7.81e-06\n",
"D0 width = 1.605e-09\n"
]
}
],
"source": [
"print(f\"pi0 width = {Particle.from_evtgen_name('pi0').width}\\n\"\n",
" f\"D0 width = {Particle.from_evtgen_name('D0').width}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\pi^0$ has a greater width than $D^0$.\n",
"If the tolerance is set to a value between their widths, the $D^0$ particle will have a constant mass while $\\pi^0$ will not."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Configuring mass functions\n",
"Consider for example the previous $D^{*+}$ example:\n",
"\n",
"\n",
"\n",
"By default, the mass function used for variable mass is the relativistic Breit-Wigner distribution. This can however be changed. If you want the mother particle to have a specific mass function for a specific decay, you can add a `zfit` parameter to the DecayLanguage dict, e.g., when it decays to two photons:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{'pi0': [{'bf': 0.988228297,\n",
" 'fs': ['gamma', 'gamma'],\n",
" 'model': 'PHSP',\n",
" 'model_params': '',\n",
" 'zfit': 'gauss'},\n",
" {'bf': 0.011738247,\n",
" 'fs': ['e+', 'e-', 'gamma'],\n",
" 'model': 'PI0_DALITZ',\n",
" 'model_params': ''},\n",
" {'bf': 3.3392e-05,\n",
" 'fs': ['e+', 'e+', 'e-', 'e-'],\n",
" 'model': 'PHSP',\n",
" 'model_params': ''},\n",
" {'bf': 6.5e-08,\n",
" 'fs': ['e+', 'e-'],\n",
" 'model': 'PHSP',\n",
" 'model_params': ''}]}\n"
]
}
],
"source": [
"dsplus_custom_mass_func = deepcopy(dsplus_chain)\n",
"dsplus_chain_subset = dsplus_custom_mass_func[\"D*+\"][1][\"fs\"][1]\n",
"# Set the mass function of pi0 to a gaussian distribution when it decays into two photons (gamma)\n",
"dsplus_chain_subset[\"pi0\"][0][\"zfit\"] = \"gauss\"\n",
"pprint(dsplus_chain_subset)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the added `zfit` field to the first decay mode of the $\\pi^0$ particle. This dict can then be passed to `GenMultiDecay.from_dict`, like before."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"GenMultiDecay.from_dict(dsplus_custom_mass_func)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you want all $\\pi^0$ particles to decay with the same mass function, you do not need to specify the `zfit` parameter for each decay in the `dict`. Instead, one can pass the `particle_model_map` parameter to the constructor:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"GenMultiDecay.from_dict(dsplus_chain, particle_model_map={'pi0': 'gauss'}) # pi0 always decays with a gaussian mass distribution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When using `DecayChain`s, the syntax for specifying the mass function becomes cleaner:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dplus_decay = DecayMode(1, \"K- pi+ pi+ pi0\", model=\"PHSP\") # The model parameter will be ignored by GenMultiDecay\n",
"pi0_decay = DecayMode(1, \"gamma gamma\", zfit=\"gauss\") # Make pi0 have a gaussian mass distribution\n",
"dplus_single = DecayChain(\"D+\", {\"D+\": dplus_decay, \"pi0\": pi0_decay})\n",
"GenMultiDecay.from_dict(dplus_single.to_dict())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Extra: Custom mass functions\n",
"The built-in supported mass function names are `gauss`, `bw`, and `relbw`, with `gauss` being the gaussian distribution, `bw` being the Breit-Wigner distribution, and `relbw` being the relativistic Breit-Wigner distribution.\n",
"\n",
"If a non-supported value for the `zfit` parameter is specified, it will automatically use the relativistic Breit-Wigner distribution. This behavior can be changed by changing the value of `GenMultiDecay.DEFAULT_MASS_FUNC` to a different string, e.g., `\"gauss\"`. If an invalid value for the `zfit` parameter is used, a `KeyError` is raised.\n",
"\n",
"It is also possible to add your own mass functions besides the built-in ones. You should then create a function that takes the mass and width of a particle and returns a mass function with the [format](https://phasespace.readthedocs.io/en/stable/usage.html#resonances-with-variable-mass) that is used for all phasespace mass functions. Below is an example of a custom gaussian distribution (implemented in the same way as the built-in gaussian distribution), which uses `zfit` PDFs:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def custom_gauss(mass, width):\n",
" particle_mass = tf.cast(mass, tf.float64)\n",
" particle_width = tf.cast(width, tf.float64)\n",
"\n",
" # This is the actual mass function that will be returned\n",
" def mass_func(min_mass, max_mass, n_events):\n",
" min_mass = tf.cast(min_mass, tf.float64)\n",
" max_mass = tf.cast(max_mass, tf.float64)\n",
" # Use a zfit PDF\n",
" pdf = zfit.pdf.Gauss(mu=particle_mass, sigma=particle_width, obs=\"\")\n",
" iterator = tf.stack([min_mass, max_mass], axis=-1)\n",
" return tf.vectorized_map(\n",
" lambda lim: pdf.sample(1, limits=(lim[0], lim[1])), iterator\n",
" )\n",
"\n",
" return mass_func"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function can then be passed to `GenMultiDecay.from_dict` as a dict, where the key specifies the `zfit` parameter name. In the example below, it is set to `\"custom_gauss\"`. However, this name can be chosen arbitrarily and does not need to be the same as the function name."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dsplus_chain_subset = dsplus_custom_mass_func[\"D*+\"][1][\"fs\"][1]\n",
"# Set the mass function of pi0 to the custom gaussian distribution\n",
"# when it decays into an electron-positron pair and a photon (gamma)\n",
"dsplus_chain_subset[\"pi0\"][1][\"zfit\"] = \"custom_gauss\"\n",
"pprint(dsplus_chain_subset)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"GenMultiDecay.from_dict(dsplus_custom_mass_func, {\"custom_gauss\": custom_gauss})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary\n",
"- PhaseSpace is a package for phase space simulations in Python\n",
"- The newly developed `GenMultiDecay` class makes it easy to load and simulate decays from DecayLanguage in PhaseSpace\n",
"- `GenMultiDecay` also makes it easy to customize resonances by setting the mass to be constant or variable, and which mass function to simulate"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.4 ('pyhep2022-phasespace')",
"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.10.4"
},
"vscode": {
"interpreter": {
"hash": "0d2d0fb24786d24986ed6bd7236c16622606db7bf454c512f84077e919c46993"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}