Skip to content

Commit

Permalink
docs: illustrate use of intensity components
Browse files Browse the repository at this point in the history
  • Loading branch information
@redeboer
redeboer committed Mar 22, 2021
1 parent 94f4c02 commit bb74f62
Show file tree
Hide file tree
Showing 2 changed files with 221 additions and 5 deletions.
3 changes: 3 additions & 0 deletions cspell.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@
"codemirror",
"colab",
"commitlint",
"darkred",
"doctest",
"doctests",
"dtype",
Expand All @@ -138,6 +139,7 @@
"imag",
"iminuit",
"indeterministic",
"isinstance",
"isort",
"jupyterlab",
"keras",
Expand Down Expand Up @@ -184,6 +186,7 @@
"seealso",
"sharex",
"sharey",
"startswith",
"subsys",
"suptitle",
"tolist",
Expand Down
223 changes: 218 additions & 5 deletions docs/usage/step3.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,10 @@
"source": [
"import pickle\n",
"\n",
"import expertsystem as es\n",
"from tensorwaves.model import SympyModel\n",
"\n",
"result = es.io.load(\"transitions.json\")\n",
"with open(\"helicity_model.pickle\", \"rb\") as stream:\n",
" model = pickle.load(stream)\n",
"with open(\"data_set.pickle\", \"rb\") as stream:\n",
Expand Down Expand Up @@ -389,8 +391,8 @@
" callback=callbacks,\n",
" use_analytic_gradient=False,\n",
")\n",
"result = minuit2.optimize(estimator, initial_parameters)\n",
"result"
"fit_result = minuit2.optimize(estimator, initial_parameters)\n",
"fit_result"
]
},
{
Expand All @@ -413,7 +415,7 @@
},
"outputs": [],
"source": [
"optimized_parameters = result[\"parameter_values\"]\n",
"optimized_parameters = fit_result[\"parameter_values\"]\n",
"for p in optimized_parameters:\n",
" print(p)\n",
" print(f\" initial: {initial_parameters[p]:.3}\")\n",
Expand Down Expand Up @@ -515,6 +517,217 @@
"compare_model(\"m_12\", data_set, phsp_set, intensity)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Intensity components"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that the {class}`~expertsystem.amplitude.helicity.HelicityModel` contains {attr}`~expertsystem.amplitude.helicity.HelicityModel.components` attribute. This is {obj}`dict` of component names to {class}`sympy.Expr <sympy.core.expr.Expr>`s helps us to identify amplitudes and intensities in the total amplitude, for instance:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"intensity_component = sum(\n",
" abs(expr) ** 2\n",
" for name, expr in model.components.items()\n",
" if \"f_{0}(1370)\" in name\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{margin}\n",
"The (analytical) intensity function you see here looks different than the data sample, because is not plotted over the phase-space as domain.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"jupyter": {
"source_hidden": true
},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"import sympy as sp\n",
"\n",
"symbol_replacements = {\n",
" **model.parameters,\n",
" sp.Symbol(\"m_012\"): 3.097,\n",
" sp.Symbol(\"m_0\"): 0,\n",
" sp.Symbol(\"m_1\"): 0.1349766,\n",
" sp.Symbol(\"m_2\"): 0.1349766,\n",
" sp.Symbol(\"phi_1+2\"): 0,\n",
" sp.Symbol(\"phi_1,1+2\"): 0,\n",
" sp.Symbol(\"theta_1+2\"): 0,\n",
" sp.Symbol(\"theta_1,1+2\"): 0,\n",
"}\n",
"p1 = sp.plot(\n",
" intensity_component.subs(symbol_replacements).doit(),\n",
" (sp.Symbol(\"m_12\"), 0, 2.5),\n",
" xlabel=R\"$m_{\\pi^0\\pi^0}$ [GeV]\",\n",
" show=False,\n",
")\n",
"p2 = sp.plot(\n",
" model.expression.subs(symbol_replacements).doit(),\n",
" (sp.Symbol(\"m_12\"), 0, 2.5),\n",
" line_color=\"darkred\",\n",
" show=False,\n",
")\n",
"p2.append(p1[0])\n",
"p2.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Just like in {ref}`usage.step2:2.2 Generate intensity-based sample`, these _intensity components_ can each be expressed in a computational backend. After that we, update the parameters with the optimized parameter values we found in {ref}`usage/step3:Perform fit`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"from typing import Sequence, Union\n",
"\n",
"import sympy as sp\n",
"from expertsystem.amplitude.helicity import HelicityModel\n",
"\n",
"\n",
"def add_components(\n",
" model: HelicityModel, components: Union[str, Sequence[str]]\n",
") -> sp.Expr:\n",
" if isinstance(components, str):\n",
" components = [components]\n",
" if any(map(lambda c: c.startswith(\"I\"), components)) and any(\n",
" map(lambda c: c.startswith(\"A\"), components)\n",
" ):\n",
" raise ValueError(\"Cannot add both intensities and amplitudes\")\n",
" if all(map(lambda c: c.startswith(\"I\"), components)):\n",
" return sum(model.components[c] for c in components)\n",
" if all(map(lambda c: c.startswith(\"A\"), components)):\n",
" return sum(abs(model.components[c]) ** 2 for c in components)\n",
" raise ValueError('Not all component names started with either \"A\" or \"I\"')\n",
"\n",
"\n",
"def create_intensity_component( # noqa: F841\n",
" model: HelicityModel,\n",
" components: Union[str, Sequence[str]],\n",
" backend: str,\n",
") -> LambdifiedFunction:\n",
" added_components = add_components(model, components)\n",
" sympy_model = SympyModel(\n",
" expression=added_components,\n",
" parameters=model.parameters,\n",
" )\n",
" return LambdifiedFunction(sympy_model, backend=backend)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"intensity_component = create_intensity_component(\n",
" model,\n",
" [\n",
" \"A[J/\\\\psi(1S)_{+1} \\\\to f_{0}(1710)_{0} \\\\gamma_{+1};f_{0}(1710)_{0} \\\\to \\\\pi^{0}_{0} \\\\pi^{0}_{0}]\",\n",
" \"A[J/\\\\psi(1S)_{+1} \\\\to f_{0}(1710)_{0} \\\\gamma_{-1};f_{0}(1710)_{0} \\\\to \\\\pi^{0}_{0} \\\\pi^{0}_{0}]\",\n",
" \"A[J/\\\\psi(1S)_{-1} \\\\to f_{0}(1710)_{0} \\\\gamma_{+1};f_{0}(1710)_{0} \\\\to \\\\pi^{0}_{0} \\\\pi^{0}_{0}]\",\n",
" \"A[J/\\\\psi(1S)_{-1} \\\\to f_{0}(1710)_{0} \\\\gamma_{-1};f_{0}(1710)_{0} \\\\to \\\\pi^{0}_{0} \\\\pi^{0}_{0}]\",\n",
" ],\n",
" backend=\"numpy\",\n",
")\n",
"intensity_component.update_parameters(latest_parameters)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or, generically for all resonances in this decay channel:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"resonances = [p.latex for p in result.get_intermediate_particles()]\n",
"components = {\n",
" f\"${p}$\": create_intensity_component(\n",
" model,\n",
" components=[c for c in model.components if p in c],\n",
" backend=\"numpy\",\n",
" )\n",
" for p in resonances\n",
"}\n",
"for c in components.values():\n",
" c.update_parameters(latest_parameters)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"jupyter": {
"source_hidden": true
},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, figsize=(10, 6))\n",
"bins = 150\n",
"ax.hist(\n",
" phsp_set[\"m_12\"],\n",
" weights=intensity(phsp_set),\n",
" bins=bins,\n",
" alpha=0.2,\n",
" label=\"full intensity\",\n",
")\n",
"ax.set_xlim(0.25, 2.5)\n",
"ax.set_xlabel(R\"$m_{\\pi^0\\pi^0}$ [GeV]\")\n",
"ax.set_yticks([])\n",
"for label, component in components.items():\n",
" ax.hist(\n",
" phsp_set[\"m_12\"],\n",
" weights=component(phsp_set),\n",
" bins=bins,\n",
" histtype=\"step\",\n",
" label=label,\n",
" )\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -641,8 +854,8 @@
" ),\n",
" use_analytic_gradient=False,\n",
")\n",
"result = scipy.optimize(estimator, initial_parameters)\n",
"result"
"fit_result = scipy.optimize(estimator, initial_parameters)\n",
"fit_result"
]
},
{
Expand Down

0 comments on commit bb74f62

Please sign in to comment.