{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "hideCode": true, "hideOutput": true, "hidePrompt": true, "jupyter": { "source_hidden": true }, "slideshow": { "slide_type": "skip" }, "tags": [ "remove-cell", "skip-execution" ] }, "outputs": [], "source": [ "# WARNING: advised to install a specific version, e.g. ampform==0.1.2\n", "%pip install -q ampform[doc,viz] IPython" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hideCode": true, "hideOutput": true, "hidePrompt": true, "jupyter": { "source_hidden": true }, "slideshow": { "slide_type": "skip" }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "%config InlineBackend.figure_formats = ['svg']\n", "import os\n", "\n", "from IPython.display import display # noqa: F401\n", "\n", "STATIC_WEB_PAGE = {\"EXECUTE_NB\", \"READTHEDOCS\"}.intersection(os.environ)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{autolink-concat}\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Formulate amplitude model" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "source_hidden": true }, "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "import sympy as sp" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Generate transitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In {doc}`qrules:usage/reaction`, we used {func}`~qrules.generate_transitions` to create a list of allowed {class}`~qrules.transition.StateTransition`s for a specific decay channel:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import qrules\n", "\n", "reaction = qrules.generate_transitions(\n", " initial_state=(\"J/psi(1S)\", [-1, +1]),\n", " final_state=[\"gamma\", \"pi0\", \"pi0\"],\n", " allowed_intermediate_particles=[\"f(0)(980)\", \"f(0)(1500)\"],\n", " allowed_interaction_types=[\"strong\", \"EM\"],\n", " formalism=\"canonical-helicity\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "source_hidden": true }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "import graphviz\n", "\n", "dot = qrules.io.asdot(reaction, collapse_graphs=True)\n", "graphviz.Source(dot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use the {class}`~qrules.transition.ReactionInfo` to formulate an amplitude model. The type of this amplitude model is dependent on the {attr}`~qrules.transition.ReactionInfo.formalism`. The function {func}`.get_builder` helps to get the correct amplitude builder class for this {attr}`~qrules.transition.ReactionInfo.formalism`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ampform import get_builder\n", "\n", "model_builder = get_builder(reaction)\n", "type(model_builder)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we now use the {meth}`.HelicityAmplitudeBuilder.formulate` method of this builder, we get a {class}`.HelicityModel` without any dynamics:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_no_dynamics = model_builder.formulate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{seealso} {doc}`formalism`\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Main expressions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A {class}`.HelicityModel` has a few attributes. The most important is its main {attr}`~.HelicityModel.expression`. This is just a {class}`sympy.Expr `, which we can pull apart by using its {attr}`~sympy.core.basic.Basic.args` (see {doc}`sympy:tutorial/manipulation`):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "full-width" ] }, "outputs": [], "source": [ "intensities = model_no_dynamics.expression.args\n", "intensity_1 = intensities[0]\n", "base, power = intensity_1.args\n", "abs_arg = base.args[0]\n", "amplitudes = abs_arg.args\n", "amplitude = amplitudes[0]\n", "amplitude" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amplitude.doit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parameters and kinematic variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As can be seen, the expression contains several {class}`~sympy.core.symbol.Symbol`s. Some of these represent (kinematic) **variables**, such as the helicity angles $\\phi_{1+2}$ and $\\theta_{1+2}$ (see {func}`.get_helicity_angle_label` for the meaning of their subscripts). Others will later on be interpreted **parameters** when fitting the model to data.\n", "\n", "The {class}`.HelicityModel` comes with expressions for these {attr}`~.HelicityModel.kinematic_variables`, so that it's possible to compute them from 4-momentum data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kinematic_variables = set(model_no_dynamics.kinematic_variables)\n", "kinematic_variables" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "theta = sp.Symbol(\"theta_1+2\", real=True)\n", "expr = model_no_dynamics.kinematic_variables[theta]\n", "expr.doit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The remaining symbols in the {class}`.HelicityModel` are parameters. Each of these parameters comes with suggested parameter values ({attr}`~.HelicityModel.parameter_defaults`), that have been extracted from the {class}`~qrules.transition.ReactionInfo` object where possible:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "full-width" ] }, "outputs": [], "source": [ "dict(model_no_dynamics.parameter_defaults)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Scalar masses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, the {class}`.HelicityAmplitudeBuilder` creates {class}`sympy.Expr `s for each kinematic variable, including all 'invariant' final state masses ($m_0, m_1, \\dots$). However, it often happens that certain particles in a final state are stable. In that case, you may want to substitute these symbols with _scalar_ values. This can be achieved by specifying which final state IDs are to be considered _stable_. Their corresponding mass symbols will then be considered parameters and a scalar suggested parameter value will be provided." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "full-width" ] }, "outputs": [], "source": [ "model_builder.stable_final_state_ids = [0, 1, 2]\n", "model_stable_masses = model_builder.formulate()\n", "dict(model_stable_masses.parameter_defaults)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "set(model_stable_masses.kinematic_variables)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can reset this option as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_builder.stable_final_state_ids = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, it can happen that the initial state mass is fixed (for instance when analyzing generated data or when applying a kinematic fit while reconstructing measured data). In that case, set {attr}`.HelicityAmplitudeBuilder.scalar_initial_state_mass` to {obj}`True`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "full-width" ] }, "outputs": [], "source": [ "model_builder.scalar_initial_state_mass = True\n", "model_stable_masses = model_builder.formulate()\n", "dict(model_stable_masses.parameter_defaults)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Extend kinematic variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The {class}`.HelicityAmplitudeBuilder` by default only generates {attr}`.kinematic_variables` (helicity angles and invariant masses) for the topologies that are available in the {class}`~qrules.transition.ReactionInfo` object that it was created with. If you want to calculate more kinematic variables, you can use the method {meth}`.register_topology` of its helicity {attr}`.HelicityAmplitudeBuilder.adapter` to register more topologies and generate more kinematic variables. This is especially useful when generating data later on with [TensorWaves](https://tensorwaves.rtfd.io).\n", "\n", "To make this a bit easier, there is {meth}`.permutate_registered_topologies`, which is a small convenience function makes it possible to generate permutations of all {attr}`.registered_topologies` and register them as well. Note that initially, only one {class}`~qrules.topology.Topology` is registered in the {attr}`.HelicityAmplitudeBuilder.adapter`, namely the one for the decay $J/\\psi \\to \\gamma f_0, f_0 \\to \\pi^0\\pi^0$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dot = qrules.io.asdot(model_builder.adapter.registered_topologies)\n", "graphviz.Source(dot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now {meth}`.permutate_registered_topologies` to register permutations of this {class}`~qrules.topology.Topology`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_builder.adapter.permutate_registered_topologies()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are now **three** {attr}`.registered_topologies`―one for each permutation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "len(model_builder.adapter.registered_topologies)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "assert len(model_builder.adapter.registered_topologies) == 3\n", "dot = qrules.io.asdot(model_builder.adapter.registered_topologies)\n", "graphviz.Source(dot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And if we {meth}`~.HelicityAmplitudeBuilder.formulate` a new {class}`.HelicityModel`, we see that it has many more {attr}`.kinematic_variables`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "set(model_builder.formulate().kinematic_variables)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "\n", "To register even more topologies, use e.g. {func}`~qrules.topology.create_isobar_topologies` to generate other, non-isomorphic topologies that cannot be created with permutations. This is relevant for more than three final states.\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set dynamics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To set dynamics for specific resonances, use {meth}`.set_dynamics` on the same {class}`.HelicityAmplitudeBuilder` instance. You can set the dynamics to be any kind of {class}`~sympy.core.expr.Expr`, as long as you keep track of which {class}`~sympy.core.symbol.Symbol` names you use (see {doc}`/usage/dynamics/custom`).\n", "\n", "AmpForm does provide a few common {mod}`.dynamics` functions, which can be constructed as {class}`~sympy.core.expr.Expr` with the correct {class}`~sympy.core.symbol.Symbol` names using {meth}`.set_dynamics`. This function takes specific {mod}`.dynamics.builder` functions and classes, such as {class}`.RelativisticBreitWignerBuilder`, which can create {func}`.relativistic_breit_wigner` functions for specific resonances. Here's an example for a relativistic Breit-Wigner _with form factor_ for the intermediate resonances and use a Blatt-Weisskopf barrier factor for the production decay:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ampform.dynamics import PhaseSpaceFactor # optional\n", "from ampform.dynamics.builder import RelativisticBreitWignerBuilder\n", "\n", "bw_builder = RelativisticBreitWignerBuilder(\n", " energy_dependent_width=True,\n", " form_factor=True,\n", " phsp_factor=PhaseSpaceFactor, # optional\n", ")\n", "for name in reaction.get_intermediate_particles().names:\n", " model_builder.set_dynamics(name, bw_builder)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this {class}`.RelativisticBreitWignerBuilder` can also be initialized with a different {class}`.PhaseSpaceFactorProtocol`. This allows us to insert different phase space factors (see {doc}`/usage/dynamics/analytic-continuation` and {func}`.create_analytic_breit_wigner`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{seealso}\n", "{doc}`/usage/dynamics/custom`\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we use the reconfigured {class}`.HelicityAmplitudeBuilder` to generate another {class}`.HelicityModel`, this time with relativistic Breit-Wigner functions and form factors:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = model_builder.formulate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Export" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is no special export function to export an {class}`.HelicityModel`. However, we can just use the built-in {mod}`pickle` module to write the model to disk and load it back:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pickle\n", "\n", "with open(\"helicity_model.pickle\", \"wb\") as stream:\n", " pickle.dump(model, stream)\n", "with open(\"helicity_model.pickle\", \"rb\") as stream:\n", " model = pickle.load(stream)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mathematical formula" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's possible to view the complete amplitude model as an expression. This would, however, clog the screen here, so we instead just view the formula of one of its {attr}`~.HelicityModel.components`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "some_amplitude = list(model.components.values())[0]\n", "some_amplitude.doit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{note}\n", "We use {meth}`~sympy.core.basic.Basic.doit` to evaluate the Wigner-$D$ ({meth}`Rotation.D `) and Clebsch-Gordan ({class}`~sympy.physics.quantum.cg.CG`) functions in the full expression.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The {attr}`.HelicityModel.parameter_defaults` attribute can be used to substitute all parameters with suggested values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "some_amplitude.doit().subs(model.parameter_defaults)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "\n", "To view the full expression for the amplitude model without crashing Jupyter Lab, install [`jupyterlab-katex`](https://pypi.org/project/jupyterlab-katex).\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case ($J/\\psi \\to \\gamma f_0, f_0 \\to \\pi^0\\pi^0$) _without dynamics_, the total intensity is only dependent on the $\\theta$ angle of $\\gamma$ in the center of mass frame (see {func}`.get_helicity_angle_label`):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "no_dynamics = model_no_dynamics.expression.doit()\n", "no_dynamics_substituted = no_dynamics.subs(model.parameter_defaults)\n", "no_dynamics_substituted" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "assert len(no_dynamics_substituted.free_symbols) == 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "theta = next(iter(no_dynamics_substituted.free_symbols))\n", "sp.plot(\n", " no_dynamics_substituted,\n", " (theta, 0, sp.pi),\n", " axis_center=(0, 0),\n", " ylabel=\"$I$\",\n", " ylim=(0, 16),\n", ");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this decay channel, the amplitude model is built up of four components:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "no_dynamics.subs(zip(no_dynamics.args, sp.symbols(\"I_{:4}\")))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This can be nicely visualized as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "plots = []\n", "colors = [\"red\", \"blue\", \"green\", \"purple\"]\n", "\n", "total = 0\n", "for i, intensity in enumerate(no_dynamics.args):\n", " total += intensity.subs(model.parameter_defaults).doit()\n", " plots.append(\n", " sp.plot(\n", " total,\n", " (theta, 0, sp.pi),\n", " axis_center=(0, 0),\n", " ylabel=\"$I$\",\n", " ylim=(0, 16),\n", " line_color=colors[i],\n", " show=False,\n", " label=f\"$I_{i}$\",\n", " legend=True,\n", " )\n", " )\n", "for i in range(1, 4):\n", " plots[0].extend(plots[i])\n", "plots[0].show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "\n", "Use {meth}`.HelicityModel.sum_components` for adding up separate components of the model.\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{tip}\n", "See {doc}`/usage/interactive` for a much more didactic way to visualize the model!\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the model _with dynamics_, we have several free symbols, such as the mass and width of the resonances. For the fitting package these will be considered **parameters** that are to be optimized and (kinematic) **variables** that represent the data set. Examples of parameters are mass ($m_\\text{particle}$) and width ($\\Gamma_\\text{particle}$) of the resonances and certain amplitude coefficients ($C$). Examples of kinematic variables are the helicity angles $\\theta$ and $\\phi$ and the invariant masses ($m_{ij...}$)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "full-width" ] }, "outputs": [], "source": [ "sorted(model.expression.free_symbols, key=lambda s: s.name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's say we want to plot the amplitude model with respect to $m_{3+4}$. We would have to substitute all other free symbols with some value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we can use {attr}`.HelicityModel.parameter_defaults` to substitute the parameters with suggested values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "suggested_expression = model.expression.subs(model.parameter_defaults)\n", "suggested_expression.free_symbols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ideally, we would now 'integrate out' the helicity angles. Here, we however just set these angles to $0$, as computing the integral would take quite some time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "angle = 0\n", "angle_substitutions = {\n", " s: angle\n", " for s in suggested_expression.free_symbols\n", " if s.name.startswith(\"phi\") or s.name.startswith(\"theta\")\n", "}\n", "evaluated_angle_intensity = suggested_expression.subs(angle_substitutions)\n", "evaluated_angle_intensity.free_symbols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By now we are only left with the masses of the final state particles ($m_1$ and $m_2$), since they appear as symbols in the {func}`.relativistic_breit_wigner_with_ff`. Final state particles 3 and 4 are the $\\pi^0$'s, so we can just substitute them with their masses. (Alternatively, see {ref}`usage/amplitude:Scalar masses`.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from qrules import load_pdg\n", "\n", "pi0 = load_pdg()[\"pi0\"]\n", "plotted_intensity = evaluated_angle_intensity.doit().subs(\n", " {\n", " sp.Symbol(\"m_1\", real=True): pi0.mass,\n", " sp.Symbol(\"m_2\", real=True): pi0.mass,\n", " },\n", " simultaneous=True,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{tip}\n", "Use {meth}`~sympy.core.basic.Basic.subs` with `simultaneous=True`, as that avoids a bug later on when plotting with {mod}`matplotlib.pyplot`.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's it! Now we are only left with the invariant mass $m_{3+4}$ of the two pions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "assert len(plotted_intensity.free_symbols) == 1\n", "m = next(iter(plotted_intensity.free_symbols))\n", "m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...and we can plot it with with {func}`sympy.plot `:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sp.plot(\n", " plotted_intensity,\n", " (m, 2 * pi0.mass, 2.5),\n", " axis_center=(2 * pi0.mass, 0),\n", " xlabel=Rf\"$m(\\pi^{0}\\pi^{0})$\",\n", " ylabel=\"$I$\",\n", " backend=\"matplotlib\",\n", ");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The expression itself looks like this (after some rounding of the {class}`float` values in this expression using {doc}`tree traversal `):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "source_hidden": true }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "def round_nested(expression: sp.Expr, n_decimals: int = 2) -> sp.Expr:\n", " for node in sp.preorder_traversal(expression):\n", " if node.free_symbols:\n", " continue\n", " if isinstance(node, (float, sp.Float)):\n", " expression = expression.subs(node, round(node, n_decimals))\n", " if isinstance(node, sp.Pow) and node.args[1] == 1 / 2:\n", " expression = expression.subs(node, round(node.n(), n_decimals))\n", " return expression" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "rounded = round_nested(plotted_intensity)\n", "round_nested(rounded)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }