diff --git a/data/crosscheck.json b/data/crosscheck.json
index 1343d36c..b9f4c359 100644
--- a/data/crosscheck.json
+++ b/data/crosscheck.json
@@ -1,4 +1,162 @@
{
+ "chains_LS": {
+ "ArD(1232)1": {
+ "A++": "(0.5027960066411277-0.5328624371775472j)",
+ "A+-": "(-0.5468818114853441+0.579584505578918j)",
+ "A-+": "(0.5468818114853441-0.579584505578918j)",
+ "A--": "(0.5027960066411277-0.5328624371775472j)"
+ },
+ "ArD(1232)2": {
+ "A++": "(-0.18048909207761316+0.19128206313914814j)",
+ "A+-": "(0.6898181946421187-0.7310682653626869j)",
+ "A-+": "(0.6898181946421187-0.7310682653626869j)",
+ "A--": "(0.18048909207761316-0.19128206313914814j)"
+ },
+ "ArD(1600)1": {
+ "A++": "(-0.29762380281702644-0.08430714758864405j)",
+ "A+-": "(0.3237198431886241+0.09169930744368116j)",
+ "A-+": "(-0.3237198431886241-0.09169930744368116j)",
+ "A--": "(-0.29762380281702644-0.08430714758864405j)"
+ },
+ "ArD(1600)2": {
+ "A++": "(0.14354077202512136+0.040660434204416125j)",
+ "A+-": "(-0.548603990834672-0.1554016754884764j)",
+ "A-+": "(-0.548603990834672-0.1554016754884764j)",
+ "A--": "(-0.14354077202512136-0.040660434204416125j)"
+ },
+ "ArD(1700)1": {
+ "A++": "(-0.04216368675863489-0.003922071424037337j)",
+ "A+-": "(-0.22655107034528332-0.021073808943055453j)",
+ "A-+": "(-0.22655107034528332-0.021073808943055453j)",
+ "A--": "(0.04216368675863489+0.003922071424037337j)"
+ },
+ "ArD(1700)2": {
+ "A++": "(-0.10534889124997639-0.009799567060888748j)",
+ "A+-": "(0.33638118886167023+0.031290220325618845j)",
+ "A-+": "(-0.33638118886167023-0.031290220325618845j)",
+ "A--": "(-0.10534889124997639-0.009799567060888748j)"
+ },
+ "ArK(892)1": {
+ "A++": "(0.21951295940623478+2.3869430597726082j)",
+ "A+-": "(-0.8577326150431347-9.326824795017867j)",
+ "A-+": "(-0.8577326150431347-9.326824795017867j)",
+ "A--": "(-0.21951295940623478-2.3869430597726082j)"
+ },
+ "ArK(892)2": {
+ "A++": "(0.2195491881362379+2.387337004238936j)",
+ "A+-": "(-0.8578741764498475-9.32836410740379j)",
+ "A-+": "(0.8578741764498475+9.32836410740379j)",
+ "A--": "(0.2195491881362379+2.387337004238936j)"
+ },
+ "ArK(892)3": {
+ "A++": "(0.3104894394702698+3.3762043693498582j)",
+ "A+-": "(0.6066086475725121+6.596149517722416j)",
+ "A-+": "(-0.6066086475725121-6.596149517722416j)",
+ "A--": "(0.3104894394702698+3.3762043693498582j)"
+ },
+ "ArK(892)4": {
+ "A++": "(0.3106291951709562+3.37772404682457j)",
+ "A+-": "(0.6068816907289188+6.5991185381800985j)",
+ "A-+": "(0.6068816907289188+6.5991185381800985j)",
+ "A--": "(-0.3106291951709562-3.37772404682457j)"
+ },
+ "ArK(1430)1": {
+ "A++": "(0.643091013936314+0.0514357428458752j)",
+ "A+-": "0j",
+ "A-+": "0j",
+ "A--": "(0.643091013936314+0.0514357428458752j)"
+ },
+ "ArK(1430)2": {
+ "A++": "(-0.6430910139363138-0.051435742845875196j)",
+ "A+-": "0j",
+ "A-+": "0j",
+ "A--": "(0.6430910139363138+0.051435742845875196j)"
+ },
+ "ArK(700)1": {
+ "A++": "(-1.070937377519722+2.282901856136404j)",
+ "A+-": "0j",
+ "A-+": "0j",
+ "A--": "(-1.070937377519722+2.282901856136404j)"
+ },
+ "ArK(700)2": {
+ "A++": "(1.0709373775197217-2.2829018561364034j)",
+ "A+-": "0j",
+ "A-+": "0j",
+ "A--": "(-1.0709373775197217+2.2829018561364034j)"
+ },
+ "ArL(1405)1": {
+ "A++": "(-0.39844416319798226+0.09729547088787999j)",
+ "A+-": "(-0.009584479616780904+0.0023404194204913516j)",
+ "A-+": "(0.009584479616780904-0.0023404194204913516j)",
+ "A--": "(-0.39844416319798226+0.09729547088787999j)"
+ },
+ "ArL(1405)2": {
+ "A++": "(0.16326954794324627-0.039868541231199255j)",
+ "A+-": "(0.31138652527942895-0.07603699941803088j)",
+ "A-+": "(0.31138652527942895-0.07603699941803088j)",
+ "A--": "(-0.16326954794324627+0.039868541231199255j)"
+ },
+ "ArL(1520)1": {
+ "A++": "(0.11738701288405862-0.13599929167700006j)",
+ "A+-": "(-0.5996268565848272+0.6947005955981446j)",
+ "A-+": "(-0.5996268565848272+0.6947005955981446j)",
+ "A--": "(-0.11738701288405862+0.13599929167700006j)"
+ },
+ "ArL(1520)2": {
+ "A++": "(0.3300063714856392-0.3823304782043633j)",
+ "A+-": "(0.27812655900829075-0.3222248704722396j)",
+ "A-+": "(-0.27812655900829075+0.3222248704722396j)",
+ "A--": "(0.3300063714856392-0.3823304782043633j)"
+ },
+ "ArL(1600)1": {
+ "A++": "(-0.43178233077154826+0.4757752694838285j)",
+ "A+-": "(0.11019851926355992-0.12142629853716529j)",
+ "A-+": "(0.11019851926355992-0.12142629853716529j)",
+ "A--": "(0.43178233077154826-0.4757752694838285j)"
+ },
+ "ArL(1600)2": {
+ "A++": "(0.10230973280389569-0.11273374852697748j)",
+ "A+-": "(-0.38914817177003963+0.42879725060122775j)",
+ "A-+": "(0.38914817177003963-0.42879725060122775j)",
+ "A--": "(0.10230973280389569-0.11273374852697748j)"
+ },
+ "ArL(1670)1": {
+ "A++": "(-0.8175657439628587+0.06182679192908746j)",
+ "A+-": "(-0.019666349597137256+0.0014872287807133656j)",
+ "A-+": "(0.019666349597137256-0.0014872287807133656j)",
+ "A--": "(-0.8175657439628587+0.06182679192908746j)"
+ },
+ "ArL(1670)2": {
+ "A++": "(0.34527078011554824-0.026110419668424232j)",
+ "A+-": "(0.6584979860302582-0.04979760743236993j)",
+ "A-+": "(0.6584979860302582-0.04979760743236993j)",
+ "A--": "(-0.34527078011554824+0.026110419668424232j)"
+ },
+ "ArL(1690)1": {
+ "A++": "(0.11030832871921545-0.07151111105183935j)",
+ "A+-": "(-0.5634680939565115+0.36528728074250294j)",
+ "A-+": "(-0.5634680939565115+0.36528728074250294j)",
+ "A--": "(-0.11030832871921545+0.07151111105183935j)"
+ },
+ "ArL(1690)2": {
+ "A++": "(0.333286554170557-0.2160642996236766j)",
+ "A+-": "(0.2808910690356812-0.1820971513620747j)",
+ "A-+": "(-0.2808910690356812+0.1820971513620747j)",
+ "A--": "(0.333286554170557-0.2160642996236766j)"
+ },
+ "ArL(2000)1": {
+ "A++": "(1.036313677049281+1.1059495419849636j)",
+ "A+-": "(0.024928279120843153+0.0266033532961457j)",
+ "A-+": "(-0.024928279120843153-0.0266033532961457j)",
+ "A--": "(1.036313677049281+1.1059495419849636j)"
+ },
+ "ArL(2000)2": {
+ "A++": "(-0.5292968632065758-0.5648633578822808j)",
+ "A+-": "(-1.0094712281097769-1.0773033948117194j)",
+ "A-+": "(-1.0094712281097769-1.0773033948117194j)",
+ "A--": "(0.5292968632065758+0.5648633578822808j)"
+ }
+ },
"chains": {
"AiD(1232)1": {
"A++": "(-0.5177095697674613-0.4884981303199729j)",
diff --git a/docs/cross-check.ipynb b/docs/cross-check.ipynb
index a58ce961..de6fec7d 100644
--- a/docs/cross-check.ipynb
+++ b/docs/cross-check.ipynb
@@ -6,13 +6,8 @@
"tags": []
},
"source": [
- "# Cross-check with LHCb data"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
+ "# Cross-check with LHCb data\n",
+ "\n",
"```{autolink-concat}\n",
"```"
]
@@ -36,13 +31,15 @@
"import json\n",
"import logging\n",
"import os\n",
+ "from functools import lru_cache\n",
+ "from textwrap import dedent\n",
"\n",
"import numpy as np\n",
"import sympy as sp\n",
- "from IPython.display import Markdown, Math, display\n",
+ "from IPython.display import Markdown, Math\n",
"from tqdm.auto import tqdm\n",
"\n",
- "from polarimetry.amplitude import simplify_latex_rendering\n",
+ "from polarimetry.amplitude import AmplitudeModel, simplify_latex_rendering\n",
"from polarimetry.data import create_data_transformer\n",
"from polarimetry.io import (\n",
" as_latex,\n",
@@ -53,22 +50,28 @@
")\n",
"from polarimetry.lhcb import (\n",
" get_conversion_factor,\n",
+ " get_conversion_factor_ls,\n",
" load_model,\n",
" load_model_builder,\n",
" parameter_key_to_symbol,\n",
")\n",
"from polarimetry.lhcb.particle import load_particles\n",
"\n",
- "mute_jax_warnings()\n",
"\n",
- "model_file = \"../data/model-definitions.yaml\"\n",
- "particles = load_particles(\"../data/particle-definitions.yaml\")\n",
- "model = load_model(model_file, particles, model_id=0)\n",
- "simplify_latex_rendering()\n",
+ "@lru_cache(maxsize=None)\n",
+ "def load_model_cached(model_id: int | str) -> AmplitudeModel:\n",
+ " return load_model(MODEL_FILE, PARTICLES, model_id)\n",
+ "\n",
"\n",
+ "mute_jax_warnings()\n",
+ "simplify_latex_rendering()\n",
"NO_TQDM = \"EXECUTE_NB\" in os.environ\n",
"if NO_TQDM:\n",
- " logging.getLogger().setLevel(logging.ERROR)"
+ " logging.getLogger().setLevel(logging.ERROR)\n",
+ "\n",
+ "MODEL_FILE = \"../data/model-definitions.yaml\"\n",
+ "PARTICLES = load_particles(\"../data/particle-definitions.yaml\")\n",
+ "DEFAULT_MODEL = load_model_cached(model_id=0)"
]
},
{
@@ -87,34 +90,50 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Lineshape comparison"
+ "## Lineshape comparison\n",
+ "\n",
+ "We compute a few lineshapes for the following point in phase space and compare it with the values from {cite}`LHCb-PAPER-2022-002`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
+ "jupyter": {
+ "source_hidden": true
+ },
+ "mystnb": {
+ "code_prompt_show": "Load phase space point"
+ },
"tags": [
- "hide-cell"
+ "hide-input"
]
},
"outputs": [],
"source": [
"σ1, σ2, σ3 = sp.symbols(\"sigma1:4\", nonnegative=True)\n",
- "lineshape_vars = {k: v for k, v in crosscheck_data[\"mainvars\"].items()}\n",
+ "lineshape_vars = crosscheck_data[\"mainvars\"]\n",
"lineshape_subs = {\n",
" σ1: lineshape_vars[\"m2kpi\"],\n",
" σ2: lineshape_vars[\"m2pk\"],\n",
- " **model.parameter_defaults,\n",
- "}"
+ " **DEFAULT_MODEL.parameter_defaults,\n",
+ "}\n",
+ "lineshape_vars"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The lineshapes are computed for the following decay chains:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "jupyter": {
- "source_hidden": true
+ "mystnb": {
+ "code_prompt_show": "Load selected decay chains"
},
"tags": [
"hide-input"
@@ -122,16 +141,23 @@
},
"outputs": [],
"source": [
- "K892_chain = model.decay.find_chain(\"K(892)\")\n",
- "L1405_chain = model.decay.find_chain(\"L(1405)\")\n",
- "L1690_chain = model.decay.find_chain(\"L(1690)\")\n",
+ "K892_chain = DEFAULT_MODEL.decay.find_chain(\"K(892)\")\n",
+ "L1405_chain = DEFAULT_MODEL.decay.find_chain(\"L(1405)\")\n",
+ "L1690_chain = DEFAULT_MODEL.decay.find_chain(\"L(1690)\")\n",
"Math(as_latex([K892_chain, L1405_chain, L1690_chain]))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "mystnb": {
+ "code_prompt_show": "Values for LHCb-PAPER-2022-002"
+ },
+ "tags": [
+ "hide-input"
+ ]
+ },
"outputs": [],
"source": [
"crosscheck_data[\"lineshapes\"]"
@@ -141,8 +167,8 @@
"cell_type": "code",
"execution_count": null,
"metadata": {
- "jupyter": {
- "source_hidden": true
+ "mystnb": {
+ "code_prompt_show": "Values as computed by this framework"
},
"tags": [
"hide-input"
@@ -150,10 +176,8 @@
},
"outputs": [],
"source": [
- "amplitude_builder = load_model_builder(model_file, particles, model_id=0)\n",
- "build_dynamics = lambda c: amplitude_builder.dynamics_choices.get_builder(c)(c)[\n",
- " 0\n",
- "].doit()\n",
+ "builder = load_model_builder(MODEL_FILE, PARTICLES, model_id=0)\n",
+ "build_dynamics = lambda c: builder.dynamics_choices.get_builder(c)(c)[0].doit()\n",
"K892_bw_val = build_dynamics(K892_chain).xreplace(lineshape_subs).n()\n",
"L1405_bw_val = build_dynamics(L1405_chain).xreplace(lineshape_subs).n()\n",
"L1690_bw_val = build_dynamics(L1690_chain).xreplace(lineshape_subs).n()\n",
@@ -167,8 +191,11 @@
"jupyter": {
"source_hidden": true
},
+ "mystnb": {
+ "code_prompt_show": "Assert that these values are equal"
+ },
"tags": [
- "remove-input"
+ "hide-input"
]
},
"outputs": [],
@@ -191,47 +218,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Amplitude comparison"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### SymPy expressions"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "half = sp.Rational(1, 2)\n",
- "amplitude_exprs = {\n",
- " (λ_Λc, λ_p): amplitude_builder.formulate_aligned_amplitude(λ_Λc, λ_p, 0, 0)[0]\n",
- " for λ_Λc in [-half, +half]\n",
- " for λ_p in [-half, +half]\n",
- "}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "unfolded_amplitude_exprs = {\n",
- " k: perform_cached_doit(expr.doit().xreplace(model.amplitudes))\n",
- " for k, expr in tqdm(amplitude_exprs.items(), disable=NO_TQDM)\n",
- "}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Numerical functions"
+ "## Amplitude comparison\n",
+ "\n",
+ "The amplitude for each decay chain and each outer state helicity combination are evaluated on the following point in phase space:"
]
},
{
@@ -241,161 +230,255 @@
"jupyter": {
"source_hidden": true
},
- "tags": []
+ "mystnb": {
+ "code_prompt_show": "Load phase space point as in DPD coordinates"
+ },
+ "tags": [
+ "hide-input"
+ ]
},
"outputs": [],
"source": [
- "%%time\n",
- "production_couplings = {\n",
- " symbol: value\n",
- " for symbol, value in model.parameter_defaults.items()\n",
- " if isinstance(symbol, sp.Indexed)\n",
- " if \"production\" in str(symbol)\n",
- "}\n",
- "fixed_parameters = {\n",
- " s: v\n",
- " for s, v in model.parameter_defaults.items()\n",
- " if s not in production_couplings\n",
+ "amplitude_vars = {k: v for k, v in crosscheck_data[\"chainvars\"].items()}\n",
+ "transformer = create_data_transformer(DEFAULT_MODEL)\n",
+ "input_data = {\n",
+ " str(σ1): amplitude_vars[\"m2kpi\"],\n",
+ " str(σ2): amplitude_vars[\"m2pk\"],\n",
+ " str(σ3): amplitude_vars[\"m2ppi\"],\n",
"}\n",
- "amplitude_funcs = {\n",
- " k: perform_cached_lambdify(\n",
- " expr.xreplace(fixed_parameters),\n",
- " parameters=production_couplings,\n",
- " backend=\"numpy\",\n",
- " )\n",
- " for k, expr in unfolded_amplitude_exprs.items()\n",
- "}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Input data"
+ "input_data = {k: float(v) for k, v in transformer(input_data).items()}\n",
+ "display_latex({sp.Symbol(k): v for k, v in input_data.items()})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
+ "jupyter": {
+ "source_hidden": true
+ },
+ "mystnb": {
+ "code_prompt_show": "Code for creating functions for each sub-amplitude"
+ },
"tags": [
+ "scroll-input",
"hide-input"
]
},
"outputs": [],
"source": [
- "amplitude_vars = {k: v for k, v in crosscheck_data[\"chainvars\"].items()}\n",
- "display(amplitude_vars)"
+ "@lru_cache(maxsize=None)\n",
+ "def create_amplitude_functions(\n",
+ " model_id: int | str,\n",
+ ") -> dict[tuple[sp.Rational, sp.Rational], sp.Expr]:\n",
+ " model = load_model(MODEL_FILE, PARTICLES, model_id)\n",
+ " production_couplings = get_production_couplings(model_id)\n",
+ " fixed_parameters = {\n",
+ " s: v\n",
+ " for s, v in model.parameter_defaults.items()\n",
+ " if s not in production_couplings\n",
+ " }\n",
+ " exprs = formulate_amplitude_expressions(model_id)\n",
+ " return {\n",
+ " k: perform_cached_lambdify(\n",
+ " expr.xreplace(fixed_parameters),\n",
+ " parameters=production_couplings,\n",
+ " backend=\"numpy\",\n",
+ " )\n",
+ " for k, expr in tqdm(exprs.items(), desc=\"Performing doit\", disable=NO_TQDM)\n",
+ " }\n",
+ "\n",
+ "\n",
+ "@lru_cache(maxsize=None)\n",
+ "def formulate_amplitude_expressions(\n",
+ " model_id: int | str,\n",
+ ") -> dict[tuple[sp.Rational, sp.Rational], sp.Expr]:\n",
+ " builder = load_model_builder(MODEL_FILE, PARTICLES, model_id)\n",
+ " half = sp.Rational(1, 2)\n",
+ " exprs = {\n",
+ " (λ_Λc, λ_p): builder.formulate_aligned_amplitude(λ_Λc, λ_p, 0, 0)[0]\n",
+ " for λ_Λc in [-half, +half]\n",
+ " for λ_p in [-half, +half]\n",
+ " }\n",
+ " model = load_model(MODEL_FILE, PARTICLES, model_id)\n",
+ " return {\n",
+ " k: perform_cached_doit(expr.doit().xreplace(model.amplitudes))\n",
+ " for k, expr in tqdm(exprs.items(), desc=\"Lambdifying\", disable=NO_TQDM)\n",
+ " }\n",
+ "\n",
+ "\n",
+ "@lru_cache(maxsize=None)\n",
+ "def get_production_couplings(model_id: int | str) -> dict[sp.Indexed, complex]:\n",
+ " model = load_model(MODEL_FILE, PARTICLES, model_id)\n",
+ " return {\n",
+ " symbol: value\n",
+ " for symbol, value in model.parameter_defaults.items()\n",
+ " if isinstance(symbol, sp.Indexed)\n",
+ " if \"production\" in str(symbol)\n",
+ " }"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
+ "jupyter": {
+ "source_hidden": true
+ },
+ "mystnb": {
+ "code_prompt_show": "Code for creating a comparison table"
+ },
"tags": [
- "hide-cell"
+ "hide-input",
+ "scroll-input"
]
},
"outputs": [],
"source": [
- "transformer = create_data_transformer(model)\n",
- "input_data = {\n",
- " str(σ1): amplitude_vars[\"m2kpi\"],\n",
- " str(σ2): amplitude_vars[\"m2pk\"],\n",
- " str(σ3): amplitude_vars[\"m2ppi\"],\n",
- "}\n",
- "input_data = {k: float(v) for k, v in transformer(input_data).items()}"
+ "def plusminus_to_helicity(plusminus: str) -> sp.Rational:\n",
+ " half = sp.Rational(1, 2)\n",
+ " if plusminus == \"+\":\n",
+ " return +half\n",
+ " if plusminus == \"-\":\n",
+ " return -half\n",
+ " raise NotImplementedError(plusminus)\n",
+ "\n",
+ "\n",
+ "def create_comparison_table(\n",
+ " model_id: int | str, decimals: int | None = None\n",
+ ") -> Markdown:\n",
+ " min_ls = not is_ls_model(model_id)\n",
+ " amplitude_funcs = create_amplitude_functions(model_id)\n",
+ " real_amp_crosscheck = {\n",
+ " k: v\n",
+ " for k, v in get_amplitude_crosscheck_data(model_id).items()\n",
+ " if k.startswith(\"Ar\")\n",
+ " }\n",
+ " production_couplings = get_production_couplings(model_id)\n",
+ " couplings_to_zero = {str(symbol): 0 for symbol in production_couplings}\n",
+ "\n",
+ " src = \"\"\n",
+ " if decimals is not None:\n",
+ " src += dedent(\n",
+ " f\"\"\"\n",
+ " :::{{tip}}\n",
+ " Computed amplitudes are equal to LHCb amplitudes up to **{decimals} decimals**.\n",
+ " :::\n",
+ " \"\"\"\n",
+ " )\n",
+ " src += dedent(\n",
+ " \"\"\"\n",
+ " | | Computed | Expected | Difference |\n",
+ " | ---:| --------:| --------:| ----------:|\n",
+ " \"\"\"\n",
+ " )\n",
+ " for i, (amp_identifier, entry) in enumerate(real_amp_crosscheck.items()):\n",
+ " coupling = parameter_key_to_symbol(\n",
+ " amp_identifier.replace(\"Ar\", \"A\"),\n",
+ " min_ls,\n",
+ " particle_definitions=PARTICLES,\n",
+ " )\n",
+ " src += f\"| **`{amp_identifier}`** | ${sp.latex(coupling)}$ |\\n\"\n",
+ " for matrix_key, expected in entry.items():\n",
+ " matrix_suffix = matrix_key[1:] # ++, +-, -+, --\n",
+ " λ_Λc, λ_p = map(plusminus_to_helicity, matrix_suffix)\n",
+ " func = amplitude_funcs[(λ_Λc, -λ_p)]\n",
+ " func.update_parameters(couplings_to_zero)\n",
+ " func.update_parameters({str(coupling): 1})\n",
+ " computed = complex(func(input_data))\n",
+ " computed *= determine_conversion_factor(coupling, λ_p, min_ls)\n",
+ " expected = complex(expected)\n",
+ " if abs(expected) != 0.0:\n",
+ " diff = abs(computed - expected) / abs(expected)\n",
+ " if diff < 1e-6:\n",
+ " diff = f\"{diff:.2e}\"\n",
+ " else:\n",
+ " diff = f'{diff:.2e}'\n",
+ " else:\n",
+ " diff = \"\"\n",
+ " src += (\n",
+ " f\"| `{matrix_key}` | {computed:>.6f} | {expected:>.6f} | {diff} |\\n\"\n",
+ " )\n",
+ " if decimals is not None:\n",
+ " np.testing.assert_array_almost_equal(\n",
+ " computed,\n",
+ " expected,\n",
+ " decimal=decimals,\n",
+ " err_msg=f\" {amp_identifier} {matrix_key}\",\n",
+ " )\n",
+ " return Markdown(src)\n",
+ "\n",
+ "\n",
+ "def determine_conversion_factor(\n",
+ " coupling: sp.Indexed, λ_p: sp.Rational, min_ls: bool\n",
+ ") -> int:\n",
+ " resonance_name = coupling.indices[0]\n",
+ " resonance = PARTICLES[str(resonance_name)]\n",
+ " if min_ls:\n",
+ " factor = get_conversion_factor(resonance)\n",
+ " else:\n",
+ " _, L, S = coupling.indices\n",
+ " factor = get_conversion_factor_ls(resonance, L, S)\n",
+ " half = sp.Rational(1, 2)\n",
+ " factor *= int((-1) ** (half + λ_p)) # # additional sign flip for amplitude\n",
+ " return factor\n",
+ "\n",
+ "\n",
+ "def is_ls_model(model_id: int | str) -> bool:\n",
+ " if isinstance(model_id, int):\n",
+ " return model_id == 17\n",
+ " return \"LS coupling\" in model_id\n",
+ "\n",
+ "\n",
+ "def get_amplitude_crosscheck_data(model_id: int | str) -> dict[str, complex]:\n",
+ " if is_ls_model(model_id):\n",
+ " return crosscheck_data[\"chains_LS\"]\n",
+ " return crosscheck_data[\"chains\"]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Default model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "jupyter": {
- "source_hidden": true
- },
"tags": [
- "remove-input"
+ "hide-input",
+ "hide-output"
]
},
"outputs": [],
"source": [
- "display_latex({sp.Symbol(k): v for k, v in input_data.items()})"
+ "create_comparison_table(model_id=0, decimals=13)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Comparison table"
+ "### LS-model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "jupyter": {
- "source_hidden": true
- },
"tags": [
- "hide-input"
+ "hide-input",
+ "hide-output"
]
},
"outputs": [],
"source": [
- "def plusminus_to_helicity(plusminus: str) -> sp.Rational:\n",
- " if plusminus == \"+\":\n",
- " return +half\n",
- " if plusminus == \"-\":\n",
- " return -half\n",
- " raise NotImplementedError(plusminus)\n",
- "\n",
- "\n",
- "amplitude_decimals = 13\n",
- "real_amp_crosscheck = {\n",
- " k: v for k, v in crosscheck_data[\"chains\"].items() if k.startswith(\"Ar\")\n",
- "}\n",
- "couplings_to_zero = {str(symbol): 0 for symbol in production_couplings}\n",
- "\n",
- "src = f\"\"\"\n",
- ":::{{tip}}\n",
- "Computed amplitudes are equal to LHCb amplitudes up to **{amplitude_decimals} decimals**.\n",
- ":::\n",
- "\n",
- "| | Computed | Expected | Difference |\n",
- "| ---:| --------:| --------:| ----------:|\n",
- "\"\"\"\n",
- "for i, (amp_identifier, entry) in enumerate(real_amp_crosscheck.items()):\n",
- " resonance_name = amp_identifier[2:-1]\n",
- " subsystem_identifier = resonance_name[0]\n",
- " coupling = parameter_key_to_symbol(amp_identifier.replace(\"Ar\", \"A\"))\n",
- " src += f\"| **`{amp_identifier}`** | ${sp.latex(coupling)}$ |\\n\"\n",
- " for matrix_key, expected in entry.items():\n",
- " matrix_suffix = matrix_key[1:] # ++, +-, -+, --\n",
- " λ_Λc, λ_p = map(plusminus_to_helicity, matrix_suffix)\n",
- " func = amplitude_funcs[(λ_Λc, -λ_p)]\n",
- " func.update_parameters(couplings_to_zero)\n",
- " func.update_parameters({str(coupling): 1})\n",
- " computed = complex(func(input_data))\n",
- " resonance = model.decay.find_chain(amp_identifier[2:-1]).resonance\n",
- " computed *= get_conversion_factor(resonance, -λ_p)\n",
- " expected = complex(expected)\n",
- " if abs(expected) != 0.0:\n",
- " diff = abs(computed - expected) / abs(expected)\n",
- " if diff < 1e-6:\n",
- " diff = f\"{diff:.2e}\"\n",
- " else:\n",
- " diff = f'{diff:.2e}'\n",
- " else:\n",
- " diff = \"\"\n",
- " np.testing.assert_array_almost_equal(\n",
- " computed,\n",
- " expected,\n",
- " decimal=amplitude_decimals,\n",
- " err_msg=f\" {amp_identifier} {matrix_key}\",\n",
- " )\n",
- " src += f\"| `{matrix_key}` | {computed:>.6f} | {expected:>.6f} | {diff} |\\n\"\n",
- "Markdown(src)"
+ "create_comparison_table(\n",
+ " \"Alternative amplitude model obtained using LS couplings\",\n",
+ " decimals=13,\n",
+ ")"
]
}
],
diff --git a/src/polarimetry/lhcb/__init__.py b/src/polarimetry/lhcb/__init__.py
index e72fbd88..8109efb1 100644
--- a/src/polarimetry/lhcb/__init__.py
+++ b/src/polarimetry/lhcb/__init__.py
@@ -451,20 +451,15 @@ def uncertainty(self) -> ParameterType:
)
-def get_conversion_factor(
- resonance: Particle, proton_helicity: sp.Rational | None = None
-) -> Literal[-1, 1]:
+def get_conversion_factor(resonance: Particle) -> Literal[-1, 1]:
# https://github.com/ComPWA/polarimetry/issues/5#issue-1220525993
half = sp.Rational(1, 2)
- factor = 1
- if proton_helicity is not None:
- factor = int((-1) ** (half - proton_helicity)) # two-particle convention
if resonance.name.startswith("D"):
- return int(-resonance.parity * factor * (-1) ** (resonance.spin - half))
+ return int(-resonance.parity * (-1) ** (resonance.spin - half))
if resonance.name.startswith("K"):
- return factor
+ return 1
if resonance.name.startswith("L"):
- return int(-resonance.parity * factor)
+ return int(-resonance.parity)
raise NotImplementedError(f"No conversion factor implemented for {resonance.name}")
@@ -472,13 +467,11 @@ def get_conversion_factor_ls(
resonance: Particle, L: sp.Rational, S: sp.Rational
) -> Literal[-1, 1]:
# https://github.com/ComPWA/polarimetry/issues/192#issuecomment-1321892494
- if resonance.name.startswith("K"):
- return int((-1) ** (L + S - sp.Rational(1, 2)))
- if resonance.name.startswith("L"):
- return int(-resonance.parity * (-1) ** (L + S - resonance.spin))
- if resonance.name.startswith("D"):
- return int(-resonance.parity * (-1) ** (L + S - sp.Rational(1, 2)))
- raise NotImplementedError(f"No conversion factor implemented for {resonance.name}")
+ if resonance.name.startswith("K") and resonance.spin == 0:
+ return 1
+ half = sp.Rational(1, 2)
+ CG_flip_factor = int((-1) ** (L + S - half))
+ return get_conversion_factor(resonance) * CG_flip_factor
def parameter_key_to_symbol(
diff --git a/src/polarimetry/lhcb/dynamics.py b/src/polarimetry/lhcb/dynamics.py
index 21b5847f..9afe79f0 100644
--- a/src/polarimetry/lhcb/dynamics.py
+++ b/src/polarimetry/lhcb/dynamics.py
@@ -3,7 +3,13 @@
import sympy as sp
from polarimetry.decay import Particle, ThreeBodyDecayChain
-from polarimetry.dynamics import BreitWignerMinL, BuggBreitWigner, FlattéSWave, Q
+from polarimetry.dynamics import (
+ BlattWeisskopf,
+ BreitWignerMinL,
+ BuggBreitWigner,
+ FlattéSWave,
+ Q,
+)
from .particle import PARTICLE_TO_ID, K, Σ, p, π
@@ -48,21 +54,33 @@ def formulate_flatte_1405(
) -> tuple[FlattéSWave, dict[sp.Symbol, float]]:
s = _get_mandelstam_s(decay)
m1, m2 = map(_to_mass_symbol, decay.decay_products)
- mass = sp.Symbol(f"m_{{{decay.resonance.name}}}")
+ m_res = sp.Symbol(f"m_{{{decay.resonance.name}}}")
Γ1 = sp.Symbol(Rf"\Gamma_{{{decay.resonance.name} \to {p.latex} {K.latex}}}")
Γ2 = sp.Symbol(Rf"\Gamma_{{{decay.resonance.name} \to {Σ.latex} {π.latex}}}")
+ m_top = _to_mass_symbol(decay.parent)
+ m_spec = _to_mass_symbol(decay.spectator)
mπ = _to_mass_symbol(π)
mΣ = sp.Symbol(f"m_{{{Σ.name}}}")
parameter_defaults = {
- mass: decay.resonance.mass,
+ m_res: decay.resonance.mass,
Γ1: decay.resonance.width,
Γ2: decay.resonance.width,
m1: decay.decay_products[0].mass,
m2: decay.decay_products[1].mass,
+ m_top: decay.parent.mass,
+ m_spec: decay.spectator.mass,
mπ: π.mass,
mΣ: Σ.mass,
}
- dynamics = FlattéSWave(s, mass, (Γ1, Γ2), (m1, m2), (mπ, mΣ))
+ l_prod = decay.incoming_ls.L
+ R_prod = sp.Symbol(R"R_{\Lambda_c}")
+ q = Q(s, m_top, m_spec)
+ q0 = Q(m_res**2, m_top, m_spec)
+ dynamics = sp.Mul(
+ (q / q0) ** l_prod,
+ BlattWeisskopf(q * R_prod, l_prod) / BlattWeisskopf(q0 * R_prod, l_prod),
+ FlattéSWave(s, m_res, (Γ1, Γ2), (m1, m2), (mπ, mΣ)),
+ )
return dynamics, parameter_defaults
diff --git a/tests/lhcb.py/test_lhcb.py b/tests/lhcb.py/test_lhcb.py
index 90691ab4..60f1128e 100644
--- a/tests/lhcb.py/test_lhcb.py
+++ b/tests/lhcb.py/test_lhcb.py
@@ -31,30 +31,30 @@ def test_get_conversion_factor_ls():
assert items == [
"L(1405) L=0 S=1/2 factor=+1",
"L(1405) L=1 S=1/2 factor=-1",
- "L(1520) L=1 S=3/2 factor=-1",
- "L(1520) L=2 S=3/2 factor=+1",
+ "L(1520) L=1 S=3/2 factor=+1",
+ "L(1520) L=2 S=3/2 factor=-1",
"L(1600) L=0 S=1/2 factor=-1",
"L(1600) L=1 S=1/2 factor=+1",
"L(1670) L=0 S=1/2 factor=+1",
"L(1670) L=1 S=1/2 factor=-1",
- "L(1690) L=1 S=3/2 factor=-1",
- "L(1690) L=2 S=3/2 factor=+1",
+ "L(1690) L=1 S=3/2 factor=+1",
+ "L(1690) L=2 S=3/2 factor=-1",
"L(2000) L=0 S=1/2 factor=+1",
"L(2000) L=1 S=1/2 factor=-1",
- "D(1232) L=1 S=3/2 factor=-1",
- "D(1232) L=2 S=3/2 factor=+1",
- "D(1600) L=1 S=3/2 factor=-1",
- "D(1600) L=2 S=3/2 factor=+1",
- "D(1700) L=1 S=3/2 factor=+1",
- "D(1700) L=2 S=3/2 factor=-1",
+ "D(1232) L=1 S=3/2 factor=+1",
+ "D(1232) L=2 S=3/2 factor=-1",
+ "D(1600) L=1 S=3/2 factor=+1",
+ "D(1600) L=2 S=3/2 factor=-1",
+ "D(1700) L=1 S=3/2 factor=-1",
+ "D(1700) L=2 S=3/2 factor=+1",
"K(700) L=0 S=1/2 factor=+1",
- "K(700) L=1 S=1/2 factor=-1",
+ "K(700) L=1 S=1/2 factor=+1",
"K(892) L=0 S=1/2 factor=+1",
"K(892) L=1 S=1/2 factor=-1",
"K(892) L=1 S=3/2 factor=+1",
"K(892) L=2 S=3/2 factor=-1",
"K(1430) L=0 S=1/2 factor=+1",
- "K(1430) L=1 S=1/2 factor=-1",
+ "K(1430) L=1 S=1/2 factor=+1",
]
@@ -87,16 +87,16 @@ def test_load_model_parameters():
assert parameters[gamma] == 0.847475
assert parameters[H_prod[Str("K(892)"), 0, +h]] == (1.0 + 0.0j) * +1
assert parameters[H_prod[Str("K(700)"), 0, +h]] == (-0.000167 - 0.68489j) * +1
- assert parameters[H_prod[Str("K(700)"), 1, +h]] == (-0.631117 + 0.040435j) * -1
+ assert parameters[H_prod[Str("K(700)"), 1, +h]] == (-0.631117 + 0.040435j) * +1
assert parameters[H_prod[Str("K(892)"), 1, +h]] == (0.341792 - 0.064047j) * -1
assert parameters[H_prod[Str("K(892)"), 1, 3 * h]] == (-0.755199 - 0.592176j) * +1
assert parameters[H_prod[Str("K(892)"), 2, 3 * h]] == (0.093754 + 0.379956j) * -1
assert parameters[H_prod[Str("K(1430)"), 0, +h]] == (-1.352114 - 3.150814j) * +1
- assert parameters[H_prod[Str("K(1430)"), 1, +h]] == (0.598156 - 0.955655j) * -1
+ assert parameters[H_prod[Str("K(1430)"), 1, +h]] == (0.598156 - 0.955655j) * +1
assert parameters[H_prod[Str("L(1405)"), 0, +h]] == (-1.224670 - 0.039521j) * +1
assert parameters[H_prod[Str("L(1405)"), 1, +h]] == (-1.811842 + 1.625622j) * -1
- assert parameters[H_prod[Str("L(1520)"), 1, 3 * h]] == (0.191708 + 0.167003j) * -1
- assert parameters[H_prod[Str("L(1520)"), 2, 3 * h]] == (0.115638 + 0.242542j) * +1
+ assert parameters[H_prod[Str("L(1520)"), 1, 3 * h]] == (0.191708 + 0.167003j) * +1
+ assert parameters[H_prod[Str("L(1520)"), 2, 3 * h]] == (0.115638 + 0.242542j) * -1
def _load_builder(model_choice: int | str) -> DalitzPlotDecompositionBuilder:
diff --git a/tests/test_io.py b/tests/test_io.py
index 920d9e84..9ef9533f 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -84,24 +84,24 @@ def test_get_readable_hash(assumptions, expected_hash, caplog: LogCaptureFixture
@pytest.mark.parametrize(
("model_id", "intensity_hash", "polarimetry_hash"),
[
- (0, 22280271, 33322092),
- (1, 22280271, 33322092),
- (2, 22280271, 33322092),
- (3, 22280271, 33322092),
- (4, 22280271, 33322092),
- (5, 22280271, 33322092),
- (6, 22280271, 33322092),
- (7, 21170486, 82478999),
- (8, 84314069, 10809627),
- (9, 57893100, 12766508),
- (10, 97621809, 15885652),
- (11, 22280271, 33322092),
- (12, 14785489, 56034405),
- (13, 87547025, 38143557),
- (14, 22280271, 33322092),
- (15, 90746888, 33929793),
- (16, 22280271, 33322092),
- (17, 63839931, 41166578),
+ (0, 72624733, 50993229),
+ (1, 72624733, 50993229),
+ (2, 72624733, 50993229),
+ (3, 72624733, 50993229),
+ (4, 72624733, 50993229),
+ (5, 72624733, 50993229),
+ (6, 72624733, 50993229),
+ (7, 42667828, 59774787),
+ (8, 42079337, 75917912),
+ (9, 49930455, 66970314),
+ (10, 62675824, 87880477),
+ (11, 72624733, 50993229),
+ (12, 57901794, 70431636),
+ (13, 89118263, 27958629),
+ (14, 72624733, 50993229),
+ (15, 55246840, 89733841),
+ (16, 72624733, 50993229),
+ (17, 96889798, 22186476),
],
)
def test_get_readable_hash_large(model_id, intensity_hash, polarimetry_hash):