diff --git a/docs/lecture02.ipynb b/docs/lecture02.ipynb
index 1c62bf8..bf85121 100644
--- a/docs/lecture02.ipynb
+++ b/docs/lecture02.ipynb
@@ -1,65 +1,1983 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "# Lecture 2 – Kinematics and Lorentz Transformations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "This is a lecture note of some sections and exercises with some answering attempts which based on [Lecture 2](https://indico.ific.uv.es/event/6803/contributions/21220) by Vincent Mathieu "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## Experiments and Frames"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Particle kinematics is the application of special relativity to elementary particle reactions.\n",
+    "In experiments, the \"beam\" collides onto the \"target\" with high energy to create other particles.\n",
+    "\n",
+    "Each particle has a 4-vector $p^{\\mu}$:\n",
+    "\n",
+    "$$\n",
+    "p^{\\mu} = (E, p_x, p_y, p_z) = (E, \\vec{p}) \n",
+    "$$\n",
+    "\n",
+    "Two frames that are commonly used are the Center of Mass (CoM) frame and the laboratory (Lab) frame.\n",
+    "\n",
+    "For a two-particle reaction, $a+b \\to 1 + ... + n$, the frames satisfy the following relations. \n",
+    "- **Lab frame**: $\\vec{p}_b = \\vec{0}$<br>\n",
+    "  Target particle $b$ is at rest, while beam particle $a$ is moving.\n",
+    "- **CoM frame**: $\\vec{p}_a+\\vec{p}_b = \\vec{0}$<br>\n",
+    "  The total momentum is zero, i.e., the momenta of the two initial state particles have opposite direction."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Masses can be defined as a Lorentz invariant quantity (unaffected by reference frame):\n",
+    "\n",
+    "$$\n",
+    " p^2_i = E^2_i = |\\vec{p}_i|^2 = m^2_i\n",
+    "$$\n",
+    "\n",
+    "where $p$ is the 4-vector, $E$ is the energy, $\\vec{p}$ is the three vector, and $m$ is the mass."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## Two particles Final States"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Let us now consider the case of a two-body decay ($n=2$), i.e., $a+b \\rightarrow 1+2$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/29308176/d1cac039-0c13-4d28-b537-4ba2f66064b6\n",
+    ":width: 300px\n",
+    "```\n",
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/29308176/7665969c-31a9-4b65-9ce9-318e2ca6fdb3\n",
+    ":width: 300px\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "When we have a final state with two particles ($n=2$), we can define one $xz$ plane that is spanned by the momenta of the two final state particles in the Lab frame. The CoM frame and the Lab frame can now be related by a boost along the $z$-axis (the direction of particle $a$).\n",
+    "\n",
+    "With this definition of the $xz$ plane, we have:\n",
+    "- $\\vec{z}$ is parallel to the beam $\\vec{p}_a$\n",
+    "- $\\vec{y} $ is parallel to $ \\vec{p}_a \\times \\vec{p}_1$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## Lorentz transformations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Lorentz transformations, can be factorized into rotations and boosts along one axis (or arbitrary direction) and are required for transforming particle properties, such as position and time, between different inertial frames in relative motion."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Boost along one axis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "When we boost in the direction of the $z$ axis, only the energy and the $z$ component of the four-momenta changes, so we can express the Lorentz boost as a 2D matrix:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "$$\\begin{pmatrix}\n",
+    " E^*\n",
+    " \\\\\n",
+    " p_z^*\n",
+    "\\end{pmatrix} =\n",
+    "\\begin{pmatrix}\n",
+    "  \\gamma & \\gamma \\beta \\\\\n",
+    "  \\gamma \\beta & \\gamma\n",
+    "\\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " E^L\\\\\n",
+    " p^L_z\n",
+    "\\end{pmatrix},\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and its inverse relation,"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\\begin{pmatrix}\n",
+    " E^L\n",
+    " \\\\\n",
+    " p_z^L\n",
+    "\\end{pmatrix} =\n",
+    "\\begin{pmatrix}\n",
+    "  \\gamma & -\\gamma \\beta \\\\\n",
+    "  -\\gamma \\beta & \\gamma\n",
+    "\\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " E^*\\\\\n",
+    " p^*_z\n",
+    "\\end{pmatrix},\n",
+    "$$ (inverse_relation)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "with "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\beta = \\frac{p^*_b}{m_b}\n",
+    "$$ (beta_label)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "$$\n",
+    "\\gamma \\beta = \\frac{E^*_b}{m_b}.\n",
+    "$$ (gammabeta_label)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{exercise}\n",
+    ":label: boost-exercise\n",
+    "Check that the above boost definition brings the target at rest in the lab frame.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{solution} boost-exercise\n",
+    ":class: dropdown\n",
+    "\n",
+    "From CoM frame to Lab frame, we apply the boost to the 4-vector of particle b (target)\n",
+    "\n",
+    "\n",
+    "$$\\begin{pmatrix}\n",
+    " E_b^*\n",
+    " \\\\\n",
+    " p_{b,z}^*\n",
+    "\\end{pmatrix} =\n",
+    "\\begin{pmatrix}\n",
+    "  \\gamma & \\gamma \\beta \\\\\n",
+    "  \\gamma \\beta & \\gamma\n",
+    "\\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " E_b^L\\\\\n",
+    " p^L_{b,z}\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "using the relations {eq}`beta_label` and {eq}`gammabeta_label` and subsititute into {eq}`inverse_relation`:\n",
+    "\n",
+    "$$\n",
+    "\\gamma = \\frac{\\gamma \\beta}{\\beta} = \\frac{E_b^* m_b}{m_b p_b^*} = \\frac{E_b^*}{p_b^*}\n",
+    "$$\n",
+    "\n",
+    "$$\\begin{pmatrix}\n",
+    " E_b^L\n",
+    " \\\\\n",
+    " p_{b,z}^L\n",
+    "\\end{pmatrix} =\n",
+    "\\begin{pmatrix}\n",
+    "  \\frac{E_b^*}{p_b^*} & -\\frac{E_b^*}{m_b} \\\\\n",
+    "  -\\frac{E_b^*}{m_b} & \\frac{E_b^*}{p_b^*}\n",
+    "\\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " E_b^*\\\\\n",
+    " p^*_{b,z}\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "Since it's only changing in z-componenet and energy, one yields:\n",
+    "\n",
+    "\n",
+    "$$\\begin{pmatrix}\n",
+    " E_b^L\n",
+    " \\\\\n",
+    " p_{b,z}^L\n",
+    "\\end{pmatrix} =\n",
+    "% \\begin{pmatrix}\n",
+    "%   \\frac{E_b^*}{p_b^*} & -\\frac{E_b^*}{m_b} \\\\\n",
+    "%   -\\frac{E_b^*}{m_b} & \\frac{E_b^*}{p_b^*}\n",
+    "% \\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " \\frac{E_b^*}{p_{b,z}^*} \\cdot E_b^* - \\frac{E_b^*}{m_b}p_{b,z}^* \\\\\n",
+    " -\\frac{E_b^*}{m_b} \\cdot E_b^* + \\frac{E_b^*}{p_{b,z}^*} \\cdot p_{b,z}^*\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "If particle $b$ is at rest after inverse-boost, that means the xyz component of momentum of particle $b$ is zero: $p_{b,x}^L = p_{b,y}^L = p_{b,z}^L=0$.\n",
+    "We have already assume $p_{b,x}^* = p_{b,y}^* = 0$, and thus $p_{b,x}^L = p_{b,y}^L = 0$.\n",
+    "\n",
+    "In order to have $p_{b,z}^L=0$,\n",
+    " $ p_{b,z}^L = -\\frac{E_b^*}{m_b} \\cdot E_b^* + \\frac{E_b^*}{p_{b,z}^*} \\cdot p_{b,z}^* = 0$ must hold.\n",
+    "\n",
+    "Therefore, only when $\\frac{E_b^*}{m_b} = 1$ or $E_b^* = m_b$ \n",
+    "\n",
+    "$ \\rightarrow p_{b,z}^L =0$.\n",
+    "Under this condition, the particle is at rest after inverse-boosted to rest frame.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "<!-- By the boost from CoM frame to Lab frame, and apply if to the 4-vector of particle b (target) -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Rotations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We use active rotations instead of passive rotations. Under an active rotation, the momentum is changed and the axes are fixed"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Example: if we rotate a momentum of unit length over the $z$ axis with an angle $\\theta$, we get:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    " p_x\n",
+    " \\\\\n",
+    " p_y\n",
+    " \\\\\n",
+    " p_z\n",
+    "\\end{pmatrix} =\n",
+    "\\begin{pmatrix}\n",
+    " \\sin(\\theta)\n",
+    " \\\\\n",
+    " 0\n",
+    " \\\\\n",
+    " \\cos(\\theta)\n",
+    "\\end{pmatrix}.\n",
+    "$$ (rotation_label)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "After a rotation of $\\omega$ around the $y$ axis, the rotated momentum forms an angle $\\theta + \\omega$ with the $z$ axis, we get:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    " p_x'\n",
+    " \\\\\n",
+    " p_y'\n",
+    " \\\\\n",
+    " p_z'\n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "R_y(\\omega)\n",
+    "\\begin{pmatrix}\n",
+    " p_x\n",
+    " \\\\\n",
+    " p_y\n",
+    " \\\\\n",
+    " p_z\n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    " \\sin(\\theta+\\omega)\n",
+    " \\\\\n",
+    " 0\n",
+    " \\\\\n",
+    " \\cos(\\theta+\\omega)\n",
+    "\\end{pmatrix}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{exercise}\n",
+    ":label: check-rotation\n",
+    "Check the result of the rotation above\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{solution} check-rotation\n",
+    ":class: dropdown\n",
+    "\n",
+    "Recall the rotation matrices, any rotation can be decomposed into rotations around $z$ and $y$ axes:\n",
+    "\n",
+    "$$\n",
+    "R_z(\\omega)\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "\\cos(\\omega) & -\\sin(\\omega) & 0 \\\\\n",
+    "\\sin(\\omega) & \\cos(\\omega) & 0 \\\\\n",
+    "0 & 0 & 1\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "R_y(\\omega)\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "\\cos(\\omega) & 0 & \\sin(\\omega) \\\\\n",
+    "0 & 1 & 0 \\\\\n",
+    "-\\sin(\\omega) &0 & \\cos(\\omega) \n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "The rotation matrix around $x$ axes is also given for completeness:\n",
+    "\n",
+    "$$\n",
+    "R_x(\\omega)\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "1&0&0 \\\\\n",
+    "0 & \\cos(\\omega) & -\\sin(\\omega) \\\\\n",
+    "0 & \\sin(\\omega) & \\cos(\\omega) \n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "With this, we can rewrite equation {eq}`rotation_label` as:\n",
+    "\n",
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    " p_x' \\\\\n",
+    " p_y' \\\\\n",
+    " p_z'\n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "R_y(\\omega)\n",
+    "\\begin{pmatrix}\n",
+    " p_x \\\\\n",
+    " p_y \\\\\n",
+    " p_z\n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "\\cos(\\omega) & 0 & \\sin(\\omega) \\\\\n",
+    "0 & 1 & 0 \\\\\n",
+    "-\\sin(\\omega) &0 & \\cos(\\omega) \n",
+    "\\end{pmatrix}\n",
+    "\\begin{pmatrix}\n",
+    " p_x \\\\\n",
+    " p_y \\\\\n",
+    " p_z\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    " p_x' \\\\\n",
+    " p_y' \\\\\n",
+    " p_z'\n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "p_x \\cdot \\cos(\\omega)+ p_z \\cdot \\sin(\\omega) \\\\\n",
+    "p_y \\\\\n",
+    "-p_x \\cdot \\sin(\\omega) + p_z \\cdot \\cos(\\omega) \n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "since \n",
+    "$\\begin{pmatrix}\n",
+    " p_x \\\\\n",
+    " p_y \\\\\n",
+    " p_z\n",
+    "\\end{pmatrix}$ =\n",
+    "$\\begin{pmatrix}\n",
+    " \\sin(\\theta) \\\\\n",
+    " 0 \\\\\n",
+    " \\cos(\\theta)\n",
+    "\\end{pmatrix}\n",
+    "$,\n",
+    "and using the trigonometric identities,\n",
+    "thus we have:\n",
+    "\n",
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    " p_x' \\\\\n",
+    " p_y' \\\\\n",
+    " p_z'\n",
+    "\\end{pmatrix} \n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "\\sin(\\theta) \\cdot \\cos(\\omega)+ \\cos(\\theta) \\cdot \\sin(\\omega) \\\\\n",
+    "0 \\\\\n",
+    "-\\sin(\\theta) \\cdot \\sin(\\omega) + \\cos(\\theta) \\cdot \\cos(\\omega) \n",
+    "\\end{pmatrix}\n",
+    "=\n",
+    "\\begin{pmatrix}\n",
+    "\\sin(\\theta+\\omega) \\\\\n",
+    "0 \\\\\n",
+    "\\cos(\\theta+\\omega) \n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Example with data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/29308176/d1cac039-0c13-4d28-b537-4ba2f66064b6\n",
+    ":width: 300px\n",
+    "```\n",
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/29308176/7665969c-31a9-4b65-9ce9-318e2ca6fdb3\n",
+    ":width: 300px\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "#### `Two-particles-1.dat`\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The data are in the lab frame, Compute and plot the $cosθ_L$ distribution.\n",
+    "\n",
+    "Boost to the CoM and Compute and plot the $cosθ^*$ distribution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "remove-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%pip install -q gdown matplotlib numpy particle pylorentz"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "mystnb": {
+     "code_prompt_show": "Import Python libraries"
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "\n",
+    "import gdown\n",
+    "import numpy as np\n",
+    "from IPython.display import display\n",
+    "\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# https://indico.ific.uv.es/event/6803/contributions/21220/attachments/11209/15506/Two-particles-1.dat\n",
+    "filename = gdown.cached_download(\n",
+    "    url=\"https://indico.ific.uv.es/event/6803/contributions/21220/attachments/11209/15506/Two-particles-1.dat\",\n",
+    "    path=\"data/Two-particles-1.dat\",\n",
+    "    md5=\"324f71ee88928915f13de5840ac5ede2\",\n",
+    "    quiet=True,\n",
+    "    verify=False,\n",
+    ")\n",
+    "data = np.loadtxt(filename)\n",
+    "data.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 4-vector for a+b->1+2\n",
+    "pa = data[::3].T\n",
+    "p1 = data[1::3].T\n",
+    "p2 = data[2::3].T\n",
+    "pa.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 4-vector conservation\n",
+    "pb = p1 + p2 - pa\n",
+    "assert pb.std(axis=1).sum().round(5) == 0\n",
+    "pb.mean(axis=1).round(3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m1 = np.sqrt(p1[0] ** 2 - p1[1] ** 2 - p1[2] ** 2 - p1[3] ** 2)\n",
+    "m1.mean(), m1.std()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# p1.shape\n",
+    "np.sqrt((p1[1:] ** 2).sum(axis=0))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "because this decay happend in x-z plane, thus we can find the $\\cos \\theta$ by:\n",
+    "\n",
+    "$$\n",
+    "\\cos \\theta = \\frac{p_{1z}}{|\\vec{p_1}|}\n",
+    "$$(cos_theta_label)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "def p3_length(four_momentum: np.ndarray) -> np.ndarray:\n",
+    "    return np.sqrt(\n",
+    "        four_momentum[1] ** 2 + four_momentum[2] ** 2 + four_momentum[3] ** 2\n",
+    "    )\n",
+    "\n",
+    "\n",
+    "p1_length = p3_length(p1)\n",
+    "cos_theta = p1[3] / p1_length"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "p1_length_vectorized = np.sqrt(np.sum(p1[1:] ** 2, axis=0))\n",
+    "np.testing.assert_array_equal(p1_length, p1_length_vectorized)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "remove-cell"
+    ]
+   },
+   "source": [
+    "##### cos $\\theta$ in lab frame"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# cos(theta) in lab frame\n",
+    "plt.title(r\"$\\cos \\theta$ distribution in lab frame\")\n",
+    "plt.xlabel(r\"$\\cos \\theta$\")\n",
+    "plt.hist(cos_theta, bins=40)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "##### cos $\\theta$ in CoM frame"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Performing boost from calculation via numpy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "The total initial 4-momentum in the lab frame is given by:\n",
+    "\n",
+    "$$\n",
+    "\\vec{P} = \\vec{p}_a + \\vec{p}_b\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# The total initial 4-momentum in lab frame\n",
+    "p_CM_lab = pa + pb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Total Energy in CM frame is the magnitude of the total initial 4-momentum:\n",
+    "\n",
+    "$$\n",
+    "E_{CM} = \\sqrt{\\vec{P} \\cdot \\vec{P}} = \\sqrt{(E_a+m_b)^2 - |\\vec{p_a}|^2} = \\sqrt{s}\n",
+    "$$\n",
+    "\n",
+    "This is valid, as we know that $\\sqrt{s}$ is the total energy in CoM frame of the colliding particles, in this case: for two particle final state of reaction $a+b \\rightarrow 1+2$\n",
+    "\n",
+    "$$\n",
+    "E_{CM} = \\sqrt{s} = \\sqrt{(E_a+E_b)^2 - |\\vec{p_a} + \\vec{p_b}|^2} = \\sqrt{(E_a+m_b)^2 - |\\vec{p_a}|^2}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "pa_length = p3_length(pa)\n",
+    "E_CM = np.sqrt(p_CM_lab[0] ** 2 - pa_length**2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "E_CM_lab = p_CM_lab[0]\n",
+    "βx_CM, βy_CM, βz_CM = β_CM = pa[1:] / E_CM_lab\n",
+    "β2 = βx_CM**2 + βy_CM**2 + βz_CM**2\n",
+    "gamma = 1 / np.sqrt(1 - β2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Find the transformed energy $E_1$ in CoM frame by Lorentz boost for $p^{0}_1$ :\n",
+    "\n",
+    "$$\n",
+    "E_{1,CM} = \\gamma(E_1 - \\vec{\\beta} \\cdot \\vec{p}_1)\n",
+    "$$\n",
+    "\n",
+    "where \n",
+    "\n",
+    "$$\n",
+    "\\gamma = \\frac{1}{\\sqrt{1-\\beta^2}}\n",
+    "$$\n",
+    "\n",
+    "and $\\beta$ is the velocity here"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "E1_CM = gamma * (p1[0] - βx_CM * p1[1] - βy_CM * p1[2] - βz_CM * p1[3])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "And now for the transformed momentum vector $\\vec{p}_1$ in CoM frame by Lorentz boost for $p_1^{i}$, where i = 1,2,3, with:\n",
+    "\n",
+    "$$\n",
+    "\\vec{p}_{1,CM} = \\vec{p}_1 + (\\frac{\\gamma -1}{\\beta^2}(\\vec{\\beta} \\cdot \\vec{p}_1)- \\gamma E_1) \\vec{\\beta}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "p1_vec = p1[1:]\n",
+    "f = β_CM * p1_vec\n",
+    "f1 = (gamma - 1) / β2\n",
+    "f2 = gamma * p1[0]\n",
+    "p1_CM_vec = p1_vec + (f1 * f - f2) * β_CM"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "p1_CM = np.concatenate((E1_CM[None], p1_CM_vec))\n",
+    "p1_CM"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similarily, the equation {eq}`cos_theta_label` is still hold, and thus we find the $\\cos \\theta_{CM}$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "p1_CM_length = p3_length(p1_CM)\n",
+    "cos_theta_CM = -p1_CM[3] / p1_CM_length"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "plt.hist(cos_theta_CM, bins=40)\n",
+    "plt.title(r\"$\\cos \\theta$ distribution in CoM frame\")\n",
+    "plt.xlabel(r\"$\\cos \\theta$\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Perfroming boost using pylorentz package below to test the boost again:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from pylorentz import Momentum4\n",
+    "\n",
+    "p1_Lab = Momentum4(*p1)\n",
+    "pa_Lab = Momentum4(*pa)\n",
+    "pb_Lab = Momentum4(*pb)\n",
+    "p2_Lab = Momentum4(*p2)\n",
+    "\n",
+    "cm = pa_Lab + pb_Lab\n",
+    "p1_CM_pyL = p1_Lab.boost_particle(-cm)\n",
+    "\n",
+    "cos_theta_CM_pyL = p1_CM_pyL.p_z / p1_CM_pyL.p\n",
+    "\n",
+    "# check the cos theta CM in 2 different ways\n",
+    "np.testing.assert_allclose(cos_theta_CM, cos_theta_CM_pyL, rtol=1e-5)\n",
+    "print(\"test past with relative tolerance 1e-5\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input",
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "plt.hist(cos_theta_CM_pyL, bins=40)\n",
+    "plt.title(r\"$\\cos \\theta$ distribution in CoM frame\")\n",
+    "plt.xlabel(r\"$\\cos \\theta$\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "jp-MarkdownHeadingCollapsed": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## Mandelstam variables\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "We have the mass-shell relation:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "$$ p_i^2 = m^2_i $$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "The Mandelstam variables s, t, and u are Lorentz invariants. Their definition are:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "$$ s =(p_a + p_b)^2 = (p_1+p_2)^2 $$\n",
+    "$$ t =(p_a - p_1)^2 = (p_b-p_2)^2 $$\n",
+    "$$ u =(p_a - p_2)^2 = (p_b-p_1)^2 $$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "\n",
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/17490173/4daa6f59-615d-4950-b7ba-e8a93f04253c\n",
+    ":width: 300px\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "s+t+u = m_a^2+m_b^2+m_1^2+m^2_2\n",
+    "$$ (label:mv)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- Check the relation in equation {eq}`label:` by using  -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!--  $$p_a+p_b=p_1+p_2$$ -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{exercise}\n",
+    ":label: mv\n",
+    "Check the relation in equation {eq}`label:mv` by using $p_a+p_b=p_1+p_2$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{solution} mv\n",
+    ":class: dropdown\n",
+    "\n",
+    "$$ s+t+u = (p_a + p_b)^2 + (p_a - p_1)^2 + (p_a - p_2)^2$$\n",
+    "\n",
+    "$$ s+t+u = p_a^2 +p_b^2 + 2 p_a p_b + p_a^2 -2 p_a p_1 + p_1^2 + p_a^2 - 2p_a p_2 +p_2^2$$\n",
+    "\n",
+    "$$\n",
+    "s+t+u = 3p_a^2 + p_b^2 +p_1^2 + p_2^2 + 2p_a p_b - 2p_a p_1 -2 p_a p_2\n",
+    "$$ (label:proving)\n",
+    "\n",
+    "from equation {eq}`label:mv` we have:\n",
+    "\n",
+    "$$ p_b = p_1 + p_2 - p_a $$\n",
+    "\n",
+    "$$ \n",
+    "p_a \\cdot p_b = p_a \\cdot p_1 + p_a \\cdot p_2 - p_a^2 \n",
+    "$$ (label:use)\n",
+    "\n",
+    "By equation {eq}`label:proving` and substitute it into equation {eq}`label:use`:\n",
+    "\n",
+    "$$\n",
+    "s+t+u = 3p_a^2 + p_b^2 +p_1^2 + p_2^2 + 2p_a \\cdot p_1 + 2p_a \\cdot p_2 - 2p_a^2 -2p_a p_1 -2 p_a p_2\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "s+t+u = p_a^2 + p_b^2 +p_1^2 + p_2^2 \n",
+    "$$\n",
+    "\n",
+    "by the mass shell relation, and thus:\n",
+    "\n",
+    "$$s+t+u = m_a^2 + m_b^2 + m_1^2 + m_2^2$$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- As following: -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "% $$ s+t+u = (p_a + p_b)^2 + (p_a - p_1)^2 + (p_a - p_2)^2$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "% $$ s+t+u = p_a^2 +p_b^2 + 2 p_a p_b + p_a^2 -2 p_a p_1 + p_1^2 + p_a^2 - 2p_a p_2 +p_2^2$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "<!-- $$\n",
+    "s+t+u = 3p_a^2 + p_b^2 +p_1^2 + p_2^2 + 2p_a p_b - 2p_a p_1 -2 p_a p_2\n",
+    "$$ (label:proving) -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- from equation {eq}`label:mv` we have: -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "% $$ p_b = p_1 + p_2 - p_a $$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- $$ \n",
+    "p_a \\cdot p_b = p_a \\cdot p_1 + p_a \\cdot p_2 - p_a^2 \n",
+    "$$ (label:use) -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "<!-- By equation {eq}`label:proving` and substitute it into equation {eq}`label:use`: -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "<!-- $$\n",
+    "s+t+u = 3p_a^2 + p_b^2 +p_1^2 + p_2^2 + 2p_a \\cdot p_1 + 2p_a \\cdot p_2 - 2p_a^2 -2p_a p_1 -2 p_a p_2\n",
+    "$$ -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- $$\n",
+    "s+t+u = p_a^2 + p_b^2 +p_1^2 + p_2^2 \n",
+    "$$ -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- by the mass shell relation, and thus: -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- $$s+t+u = m_a^2 + m_b^2 + m_1^2 + m_2^2$$ -->"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "jp-MarkdownHeadingCollapsed": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### CoM Kinematics"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Every frame dependent quantities is expressed with Mandelstam variables\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "```{image} https://user-images.githubusercontent.com/17490173/278340815-e5eb448f-e620-4b67-873a-e18cd146c55a.png\n",
+    ":width: 500px\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$p_a = (E_a^*,0,0,|\\vec{p_a}^*|)$\n",
+    "\n",
+    "$p_b = (E_b^*,0,0,-|\\vec{p_a}^*|)$\n",
+    "\n",
+    "$p_1 = (E_1^*, |\\vec{p_1}^*|sin \\theta^*,0, |\\vec{p_1}^*|cos \\theta^*)$\n",
+    "\n",
+    "$p_2 = (E_2^*, -|\\vec{p_1}^*| sin \\theta^*,0, -|\\vec{p_1}^*|cos \\theta^*)$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Invariant: \n",
+    "$$s = (p_a + p_b)^2 = (E_a^* + E_b^*)^2 = E_{CM}^2$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$ t = (p_a - p_1)^2 = m_a^2 + m_1^2 - 2E_a^* E_1^* + 2 |\\vec{p_a}^*| |\\vec{p_1}^*| cos \\theta^*$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$p_1^2 = [ (p_a+p_b) - p_2 ]^2$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{exercise}\n",
+    ":label: E2*\n",
+    "Check that $E_1^* = \\frac{s+m^2_1-m_2^2}{2 \\sqrt{s}}$ and $E_2^* = \\frac{s+m^2_2-m_1^2}{2 \\sqrt{s}}$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{solution} E2*\n",
+    ":class: dropdown\n",
+    "\n",
+    "$$ \n",
+    "E_2^* = E_a^* + E_b^* - E_1^* = E^*_{CM} - E_1^* = \\sqrt{s} - E_1^* = \\sqrt{s} - \\sqrt{\\vec{p_1}^2 + m_1^2}\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "E_2^{*2} = s + \\vec{p_1}^2 + m_1^2 - 2\\sqrt{s}\\sqrt{\\vec{p}^2+m_1^2}\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "E_2^{*2} = \\vec{p_1}^2 + m_2^2\n",
+    "$$\n",
+    "\n",
+    "from the 2 equations of $E_2^{*2}$ above, we have:\n",
+    "$$\n",
+    "s + \\vec{p_1}^2 + m_1^2 - 2\\sqrt{s}\\sqrt{\\vec{p}^2+m_1^2} = \\vec{p_1}^2 + m_2^2\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "s + m_1^2 - 2\\sqrt{s}\\sqrt{\\vec{p}^2+m_1^2} =  m_2^2\n",
+    "$$\n",
+    "\n",
+    "$$\n",
+    "s + m_1^2 - m_2^2 =  2\\sqrt{s}E_1^* \n",
+    "$$\n",
+    "\n",
+    "We have E_1^* :\n",
+    "$$\n",
+    "E_1^*  = \\frac{s + m_1^2 - m_2^2 }{2\\sqrt{s}}\n",
+    "$$\n",
+    "\n",
+    "\n",
+    "For E_2^* : \n",
+    "$$\n",
+    "E_2^* = \\sqrt{s} - E_1^* = \\sqrt{s} - \\frac{s + m_1^2 - m_2^2 }{2\\sqrt{s}} = \\frac{s + m_2^2 - m_1^2 }{2\\sqrt{s}}\n",
+    "$$\n",
+    "\n",
+    "\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{exercise}\n",
+    ":label: p1*\n",
+    "Check that $|\\vec{p_1}^*| = \\sqrt{(E_1*)^2-m_1^2} = \\frac{\\lambda^{1/2}(s,m_1^2,m_2^2)}{2 \\sqrt{s}}$, where $\\lambda(a,b,c) = a^2+b^2+c^2 - 2(ab+bc+ca)$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{solution} p1*\n",
+    ":class: dropdown\n",
+    "\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "% $$t = m_a^2 +m_1^2 -2 E^*_a E^*_1 + 2|\\vec{p_a}^*| |\\vec{p_1}^*| cos \\theta^*$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$ t_{min,max} = m_a^2 + m_1^2 - 2E_a^* E_1^* + 2 |\\vec{p_a}^*| |\\vec{p_1}^*| $$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Equivalent expression for the scattering angle in the CoM"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "cos \\theta^* = 1 + \\frac{t-t_{min}}{2|\\vec{v}_a^*||\\vec{p}_1^*|}\n",
+    "=\n",
+    "\\frac{s(t-u)+(m_a^2-m-b^2)(m_1^2-m_2^2)}{\\lambda^{1/2}(s,m_a^2,m_b^2)\\lambda^{1/2}(s,m_1^2,m_2^2)}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compute $cos \\theta^L$ as a function of s,t,u"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Lecture 2 – Kinematics"
+    "$$ cos \\theta^L $$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "source": [
+    "## Three particles final states"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Relevant Variables"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "source": [
-    "[Lecture 2](https://indico.ific.uv.es/event/6803/contributions/21220) by Vincent Mathieu contains a few data files containing four-momenta data samples. Our goal in this notebook is to identify which reaction was used to generate these data samples."
+    "## Reference Frames "
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "For scattering study in 3 body final states,\n",
+    "the Gottfried-Jackson frame and the helicity frame are two different  particular center-of-mass reference frames used to analyze the angular distribution and polarization of particles in high-energy physics experiments, they are particularly useful\n",
+    "in the study of scattering processes.\n",
+    "\n",
+    "Gottfried-Jackson and helicity frames\n",
+    "are related by a rotation of angle ω around y"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Gottfried-Jackson Frame"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/17490173/d28e6a93-033c-42ad-abec-0387de3a3eac\n",
+    ":width: 300px\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "It is particularly useful when the particles involved have non-zero spin, and their spin orientations are aligned with their motion, making the analysis of their angular distributions relative to their spin directions more straightforward"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Helicity Frame"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "```{image} https://github.com/ComPWA/strong2020-salamanca/assets/17490173/3fe1ca26-5346-444b-a770-435820cf3da7\n",
+    ":width: 300px\n",
+    "```\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "This frame is often used in the analysis of the decay products of a particle, where it's important to understand the spin orientation relative to the direction of motion."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
     },
     "tags": [
      "remove-cell"
     ]
    },
-   "outputs": [],
    "source": [
-    "%pip install -q gdown matplotlib numpy particle"
+    "## The Decay"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "source": [
-    "## `Three-particles-1.dat`"
+    "### Decay into two particles"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
     },
-    "mystnb": {
-     "code_prompt_show": "Import Python libraries"
+    "tags": []
+   },
+   "source": [
+    "### Decay into three particles\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "$$\\theta^*_{12}+\\theta^*_{23}+\\theta^*_{31}=2 \\pi$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
     },
     "tags": [
-     "hide-cell"
+     "remove-cell"
     ]
    },
-   "outputs": [],
    "source": [
-    "import warnings\n",
+    ":::{exercise}\n",
+    ":label: 2pi\n",
+    "Check relation above using $s_{12}+s_{23}+s_{31}$\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Dalitz Plot"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "### Exercise: `Three-particles-1.dat`"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[Lecture 2](https://indico.ific.uv.es/event/6803/contributions/21220) by Vincent Mathieu contains a few data files containing four-momenta data samples. \n",
+    "Our goal in this notebook is to identify which reaction was used to generate these data samples."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{exercise}\n",
+    ":label: three-particles-1\n",
+    "File: Three-Particles-1.dat is \n",
+    "The data are in the lab frame, \n",
+    "with Format:\n",
     "\n",
-    "import gdown\n",
-    "import numpy as np\n",
-    "from IPython.display import display\n",
+    "Ea, pa,x, pa,y, pa,z\n",
     "\n",
-    "warnings.filterwarnings(\"ignore\")"
+    "E1, p1,x, p1,y, p1,z\n",
+    "\n",
+    "E2, p2,x, p2,y, p2,z \n",
+    "\n",
+    "E3, p3,x, p3,y, p3,z\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    "\n",
+    "Three questions to be found out.\n",
+    "\n",
+    "Which reaction is it?\n",
+    "\n",
+    "What resonances are included? \n",
+    "\n",
+    "What is the spin of the resonances?\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    ":::{solution-start} three-particles-1\n",
+    ":class: dropdown\n",
+    ":::"
    ]
   },
   {
@@ -86,6 +2004,21 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
+   "source": [
+    "print(data)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
    "source": [
     "n_final_state = 3\n",
     "pa, p1, p2, p3 = (data[i::4].T for i in range(n_final_state + 1))\n",
@@ -194,6 +2127,14 @@
     "So this is a photon&nbsp;$\\gamma$ hitting a proton&nbsp;$p$ and producing a meson&nbsp;$\\eta$, pion&nbsp;$\\pi^0$, and proton&nbsp;$p$."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{solution-end} \n",
+    ":::"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -262,9 +2203,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
     "tags": [
      "hide-input",
      "full-width"
@@ -290,6 +2228,45 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input",
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2, ax3) = plt.subplots(figsize=(12, 4), ncols=3)\n",
+    "fig.suptitle(\"1D histogram of $s_{12}, s_{23}$, and $s_{31}$\")\n",
+    "ax1.hist(s12, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax1.set_xlabel(R\"$s_{12}$\")\n",
+    "ax1.set_ylabel(\"counts\")\n",
+    "ax1.axvline(R12, c=\"C0\", ls=\"dashed\", label=\"$R_{12}$\")\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.hist(s23, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax2.set_xlabel(R\"$s_{23}$\")\n",
+    "ax2.set_ylabel(\"counts\")\n",
+    "ax2.axvline(R23, c=\"C1\", ls=\"dashed\", label=\"$R_{23}$\")\n",
+    "ax2.legend()\n",
+    "\n",
+    "ax3.hist(s31, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax3.set_xlabel(R\"$s_{31}$\")\n",
+    "ax3.set_ylabel(\"counts\")\n",
+    "ax3.axvline(R31, c=\"C2\", ls=\"dashed\", label=\"$R_{31}$\")\n",
+    "ax3.legend()\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -372,7 +2349,42 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## `Three-particles-2.dat`"
+    "### Exercise: `Three-particles-2.dat`"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{exercise}\n",
+    ":label: three-particles-2\n",
+    "File: Three-Particles-2.dat is also\n",
+    "The data are in the lab frame, \n",
+    "with Format:\n",
+    "\n",
+    "Ea, pa,x, pa,y, pa,z\n",
+    "\n",
+    "E1, p1,x, p1,y, p1,z\n",
+    "\n",
+    "E2, p2,x, p2,y, p2,z \n",
+    "\n",
+    "E3, p3,x, p3,y, p3,z\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    ".    .    .    .\n",
+    "\n",
+    "\n",
+    "Three questions to be found out.\n",
+    "\n",
+    "Which reaction is it?\n",
+    "\n",
+    "What resonances are included? \n",
+    "\n",
+    "What is the spin of the resonances?\n",
+    ":::"
    ]
   },
   {
@@ -435,9 +2447,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
     "mystnb": {
      "code_prompt_show": "Identify final state particles"
     },
@@ -489,9 +2498,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
     "tags": [
      "hide-input"
     ]
@@ -513,9 +2519,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
     "tags": [
      "hide-input",
      "full-width"
@@ -540,6 +2543,48 @@
     "fig.tight_layout()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide-input",
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "fig, (ax1, ax2, ax3) = plt.subplots(figsize=(12, 4), ncols=3)\n",
+    "fig.suptitle(\"1D histogram of $s_{12}, s_{23}$, and $s_{31}$\")\n",
+    "ax1.hist(s12, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax1.set_xlabel(R\"$s_{12}$\")\n",
+    "ax1.set_ylabel(\"counts\")\n",
+    "ax1.axvline(R12, c=\"C0\", ls=\"dashed\", label=\"$R_{12}$\")\n",
+    "ax1.legend()\n",
+    "\n",
+    "ax2.hist(s23, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax2.set_xlabel(R\"$s_{23}$\")\n",
+    "ax2.set_ylabel(\"counts\")\n",
+    "ax2.axvline(R23, c=\"C1\", ls=\"dashed\", label=\"$R_{23}$\")\n",
+    "ax2.legend()\n",
+    "\n",
+    "ax3.hist(s31, bins=100, color=\"black\", histtype=\"step\")\n",
+    "ax3.set_xlabel(R\"$s_{31}$\")\n",
+    "ax3.set_ylabel(\"counts\")\n",
+    "ax3.axvline(R31, c=\"C2\", ls=\"dashed\", label=\"$R_{31}$\")\n",
+    "ax3.legend()\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -561,7 +2606,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,