3.0 Resonance Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Pkg\n",
    "Pkg.activate(@__DIR__)\n",
    "Pkg.instantiate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DelimitedFiles # read and write plain txt files\n",
    "using TypedTables # form tables\n",
    "using Plots # plotting \n",
    "using Optim # fitting\n",
    "using Parameters # @unpack, working with the named tupples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theme(:wong)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resonance analysis\n",
    "\n",
    "1. Intro: complex functions\n",
    "2. A pole of $a_2(1320)$\n",
    "3. Advancing the model: good fit\n",
    "4. Extending the range: the pole of $a_2'(1700)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.0 Intro: Scattering amplitude\n",
    "The observed intensity (cross section / differential width) of the QM process:\n",
    "$$\n",
    "\\begin{align}\n",
    "I(s) = \\frac{1}{J} \\int |M|^2 d\\Phi\n",
    "\\end{align}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "M = \\left\\langle \\text{assymp. final state}| \\hat{T} | \\text{assymp. initial state}\\right\\rangle\n",
    "$$\n",
    "(there is also delta function of the what we will often omit):\n",
    "$$\n",
    "(2\\pi)^4\\delta^4(\\text{energy-mom. conservation})\n",
    "$$\n",
    "\n",
    "$M$ is the matrix element, or **scattering amplitude** - a dynamic function containing all information about the scattering.\n",
    "It can also be called\n",
    " * transition amplitude\n",
    " * production amplitude: the particle species in the initial state are different from the final one\n",
    "\n",
    "Two things to distinguish:\n",
    " - **Elastic** process: $\\pi\\pi\\to\\pi\\pi$ in $P$-wave below $1\\,$GeV - no other possible outcomes of the initial states\n",
    " - **Inelastic** process: $\\pi\\pi\\to K\\bar{K}$ in $S$-wave above $1\\,$GeV - there are other possibilities for the final state when you bring two particles together.\n",
    " \n",
    " ***Wave***??? - *the quantum numbers of the system) - the total angural momenum $J$. Concerved!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The phase space\n",
    "$\\Phi$ is the number of the possible configrations of the system.\n",
    "\n",
    "Two particle system of $m_1 = 0$, $m_2 = 0$, with the energy $E = 5\\,$GeV\n",
    " - $E_1 = 2.5\\,$GeV, $E_2=2.5\\,$GeV\n",
    " - the phase space is equal to the possible configutations are their angle\n",
    "\n",
    "$$\n",
    "\\Phi = n_\\Phi d\\Omega = n_\\Phi (4\\pi)\n",
    "$$\n",
    "the full solid angle, ($n_\\Phi = 1/(8\\pi)/(4\\pi)$ - the normalization)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Φ = 1/(8π)\n",
    "M = 3+1im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Intensity = abs2(M)*Φ # not the probability yet - more like cross section (the flux is missing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What the $M$ depends on? \n",
    "  - most importantly, on the **energy of the system**. The energy is conserved between the initial and the final state.\n",
    "$$\n",
    "s = (p_1+p_2)^2,\\qquad\\text{still the energy of the system, but in the Lorentz invariant form}\n",
    "$$\n",
    "Analog of the frequency: the system might start resonating if we ajust the \"frequency\" well to fit the mass of the other particle.\n",
    "\n",
    "Mathematically, $M(s)$ is a complex function. The poles of $M$ are the resonances: they correspond to other particles.\n",
    "$$\n",
    "\\begin{align}\n",
    "    M(s) \\xrightarrow[s\\to s_\\text{pole}^R]{} \\frac{\\text{res}}{s_\\text{pole}^R-s},\n",
    "    \\qquad\n",
    "    s_\\text{pole}^R = (m_R - \\frac{i\\Gamma}{2})^2.\n",
    "\\end{align}\n",
    "$$\n",
    "$s_\\text{pole}^R$ is universal propery of the resonance: does not matter which system resonates - the pole is fixed for $R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mres(s) = 1/(s^2+3^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(Mres, -30, 30, lab=\"my first amplitude\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    xv = range(-10,10, length=100)\n",
    "    yv = range(-10,10, length=100)\n",
    "    calv = [Mres(x+1im*y) for y in yv, x in xv]\n",
    "    plot3d(xv,yv,log.(abs2.(calv)))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise:*** Come up with an amplitude that\n",
    " - have a pole at the at $(m-i\\Gamma/2) = (1.3-0.01i)$\n",
    " - have two poles below the real axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.0 Example: COMPASS $\\eta\\pi$ scattering. $a_2\\to \\eta\\pi$\n",
    "We have been given a clean sample of the data for $\\eta\\pi$ scattering with the orbital momentum $L=2$:\n",
    "$$\n",
    "J^P:\\qquad \\underbrace{0^-}_{\\eta} \\otimes \\underbrace{0^-}_{\\pi} \\otimes \\underbrace{2^+}_{L=2} = 2^+\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const tdata = let\n",
    "    M = readdlm(joinpath(\"data\", \"EtaPi-2pp_arxiv1408.4286.txt\"))\n",
    "    Table(mηπ = M[:,1], I = M[:,2], δI = M[:,3])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A(s; m, Γ) = m*Γ/(m^2-s-1im*m*Γ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black), lab=\"data\")\n",
    "plot!(mηπ->abs2(10*A(mηπ^2; m=1.1, Γ=0.4)), 1.0,2.0, l=:red, lab=\"simple model\", title=\"hands adjustment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Exercise:*** adjust parameters $m$, $\\Gamma$, and normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Advancing the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const mπ = 0.140; # GeV\n",
    "const mη = 0.547; # GeV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ(x,y,z) = x^2+y^2+z^2-2x*y-2y*z-2z*x\n",
    "function model(mηπ; pars)\n",
    "    @unpack N0, m, Γ = pars\n",
    "    #\n",
    "    psq(x) = λ(x^2,mη^2,mπ^2)/(4x^2)\n",
    "    bwf(z) = 1/(9+3z+z^2)\n",
    "    FF(x;R = 3) = bwf(R^2*psq(x))/bwf(R^2*psq(m))\n",
    "    #\n",
    "    return abs2(N0*A(mηπ^2; m=m, Γ=Γ))*(psq(mηπ)/psq(m))^2*FF(mηπ)*sqrt(λ(mηπ^2,mη^2,mπ^2))/mηπ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let pars = (N0 = 350, m=1.32,Γ=0.12)\n",
    "    scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black), lab=\"data\")\n",
    "    plot!(mηπ->model(mηπ; pars=pars), 1.0,2.0, l=:red, lab=\"model\", title=\"hands adjustment\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss(; data, pars) = sum((y - model(e; pars=pars))^2/δy^2 for (e, y, δy) in data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss(; data=filter(x->1.0<x.mηπ<1.5, tdata), pars=(m=1.32,Γ=0.12,N0=350))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fr = Optim.optimize(x->loss(; data=filter(x->1.0<x.mηπ<1.5, tdata), pars=(m=x[1],Γ=x[2],N0=x[3])),\n",
    "    [1.32,0.12,350])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let pars=NamedTuple{(:m,:Γ,:N0)}(Optim.minimizer(fr))\n",
    "    scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black), label=\"data\")\n",
    "    plot!(mηπ->model(mηπ;pars=pars), 1.0,1.5, l=:red, label=\"model\", title=\"fit\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestmodel(mηπ) = model(mηπ; pars=NamedTuple{(:m,:Γ,:N0)}(Optim.minimizer(fr)))\n",
    "bestA(mηπ) = A(mηπ^2; NamedTuple{(:m,:Γ)}(Optim.minimizer(fr)[1:2])...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    plot()\n",
    "    plot!(x->real(bestA(x)), 1.0,1.5, lab=\"Re\")\n",
    "    plot!(x->imag(bestA(x)), 1.0,1.5, lab=\"Im\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    realmv = range(1.0,1.5, length=50)\n",
    "    imagmv = range(-0.3,0, length=50)\n",
    "    calv = [bestA(r+1im*i) for i in imagmv, r in realmv]\n",
    "    calv = log10.(abs2.(calv))\n",
    "    contour(realmv, imagmv, calv)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pole_fr = Optim.optimize(x->abs2(1/bestA(x[1]+1im*x[2])), [1.3, -0.06])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Optim.minimizer(pole_fr) .* [1, -2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $a_2'$ getting the higher-mass resonances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. A general expression of the amplitude: numerator over the denominator (hehe)\n",
    "$$\n",
    "\\begin{align}\n",
    "A = \\frac{N(s,\\dots)}{D(s)}\n",
    "\\end{align}\n",
    "$$\n",
    "The numerator is production dependent, the denominator contains the resonance poles.\n",
    "\n",
    "The general scattering theory controls only imaginary part of the denominator!\n",
    "In the $\\eta\\pi$ example, Im$D=\\text{const}$, does not change with $s$:\n",
    "$$\n",
    "\\begin{align}\n",
    "A^{(\\text{simple})} = \\frac{1}{(m^2-s)/(m\\Gamma) - i}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "2. Incorporating the \"unknown physics\" in the parametrization of the numerator and the **demonimator**\n",
    " - once should do both, leaving as little freedom as possible in $D$,\n",
    " while $N$ must be a smouth function. No poles\n",
    "$$\n",
    "\\begin{align}\n",
    "A^{(\\text{adv.})} = \\frac{N(s)}{R(s) - i}\n",
    "\\end{align}\n",
    "$$\n",
    "Let's put a polynomical in $s$ in the $R(s)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A2(s; ps) = 1/(sum(p*s^(i-1) for (i,p) in enumerate(ps))-1im)\n",
    "function model2(mηπ; pars)\n",
    "    @unpack N0, ps = pars\n",
    "    #\n",
    "    psq(x) = λ(x^2,mη^2,mπ^2)/(4x^2)\n",
    "    bwf(z) = 1/(9+3z+z^2)\n",
    "    FF(x;R = 3) = bwf(R^2*psq(x))\n",
    "    #\n",
    "    return abs2(N0*A2(mηπ^2; ps=ps))*(psq(mηπ))^2*FF(mηπ)*sqrt(λ(mηπ^2,mη^2,mπ^2))/mηπ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let \n",
    "    m, Γ, N0 = Optim.minimizer(fr)\n",
    "    pars=(ps = [m^2/(m*Γ), -1/(m*Γ), 0], N0=17*N0)\n",
    "    scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black), label=\"data\")\n",
    "    plot!(mηπ->model2(mηπ;pars=pars), 1.0,1.5, l=:red, label=\"model\", title=\"fit\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function loss2(; data, pars)\n",
    "    return sum((y - model2(e; pars=pars))^2/δy^2 for (e, y, δy) in data)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fr2 = Optim.optimize(x->loss2(; data=filter(x->1.0<x.mηπ<2.0, tdata), pars=(N0=x[1],ps=x[2:end])),\n",
    "#     [350,1.32/0.12,-1/(1.32*0.12),rand(),rand()])\n",
    "    [350,rand(4)...])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function makeanalysis(settings)\n",
    "    data = settings[\"data\"]\n",
    "    Nps = settings[\"Nps\"]\n",
    "    Natt = settings[\"Natt\"]\n",
    "    frs = [Optim.optimize(x->loss2(; data=data, pars=(N0=x[1],ps=x[2:end])),\n",
    "        [350,rand(Nps)...]) for e in 1:Natt]\n",
    "    settings[\"fit_resultls\"] = frs\n",
    "    return settings\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = Dict(\n",
    "    \"data\"=>filter(x->1.0<x.mηπ<2.0, tdata),\n",
    "    \"Nps\"=>6,\n",
    "    \"Natt\"=>500)\n",
    "makeanalysis(settings);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestmin = findmin(Optim.minimum.(settings[\"fit_resultls\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestpars = Optim.minimizer(settings[\"fit_resultls\"][bestmin[2]])\n",
    "bestmodel2(mηπ) = model2(mηπ; pars=(N0=bestpars[1], ps=bestpars[2:end]))\n",
    "bestA2(mηπ) = A2(mηπ^2; ps=bestpars[2:end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let \n",
    "    pars=(N0=bestpars[1], ps=bestpars[2:end])\n",
    "    scatter(tdata.mηπ, tdata.I, yerr=tdata.δI, lc=:black, m=(2,:black), label=\"data\")\n",
    "    xmin, xmax = settings[\"data\"].mηπ[[1,end]]\n",
    "    plot!(bestmodel2, xmin, xmax, l=:red, label=\"model\", title=\"fit\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    realmv = range(1.0,2.2, length=50)\n",
    "    imagmv = range(-0.5,0, length=50)\n",
    "    calv = [bestA2(r+1im*i) for i in imagmv, r in realmv]\n",
    "    calv = log10.(abs2.(calv))\n",
    "    contour(realmv, imagmv, calv)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2′pole_fr = Optim.optimize(x->abs2(1/bestA2(x[1]+1im*x[2])), [1.9, -0.2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Optim.minimizer(a2′pole_fr) .* [1, -2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See mode details in the research paper [Tensow wave by JPAC and COMPASS](https://inspirehep.net/literature/1665092)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.0",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}