2.0 The partial-wave analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partial Wave Analysis of $\\eta\\pi$"
   ]
  },
  {
   "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 TypedTables\n",
    "using DelimitedFiles\n",
    "using Plots\n",
    "using PartialWaveFunctions\n",
    "using Cuba\n",
    "using QuadGK\n",
    "using Optim\n",
    "using Test\n",
    "using ForwardDiff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = let\n",
    "    M = readdlm(joinpath(\"data\",\"data_metapi_costheta_phi.txt\"))\n",
    "    Table(mηπ=M[:,1], cosθ=M[:,2], ϕ=M[:,3])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basis function are the orthonormal rotation functions\n",
    "$$\n",
    "\\Psi_{LM}(\\Omega) = \\sqrt{\\frac{2L+1}{2\\pi}}d_{M0}^{L}(\\theta)\\,\\sin(M\\phi)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ΨLM(cosθ,ϕ,L,M) = sqrt((2L+1)/(4π))*wignerd(L,M,0,cosθ)*sin(M*ϕ)*sqrt(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can check orthogonality of the functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let M=1, L=5\n",
    "    4π*cuhre((x,f)->f[1] = abs2(ΨLM(2x[1]-1,π*(2*x[2]-1),L,M)),2,1).integral[1]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const LMs = NamedTuple{(:L,:M)}.([(1,1),(2,1),(4,1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model of the amplitude is a truncated series of the partial waves.\n",
    "\n",
    "We truncate at $L<5$ and $M<2$. $L=0$ and $M=0$ are forbidden by the selection rule (*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model(cosθ,ϕ; pars) = sum(c*ΨLM(cosθ,ϕ,L,M) for (c,(L,M)) in zip(pars, LMs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@test model(0.1,0.1; pars=rand(length(LMs)) + 1im*rand(length(LMs))) != 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the model in $\\cos\\theta \\times \\phi$ coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let pars = rand(length(LMs)) + 1im*rand(length(LMs))\n",
    "    cosθv = range(-1,1,length=100)\n",
    "    ϕv = range(-π,π,length=100)\n",
    "    calv = [model(cosθ,ϕ; pars=pars) for ϕ in ϕv, cosθ in cosθv]\n",
    "    heatmap(cosθv, ϕv, abs2.(calv), lab=\"\",\n",
    "        color=cgrad(:viridis, scale=:exp))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellh(;data, pars) = -sum(log, abs2(model(cosθ,ϕ; pars=pars)) for (_,cosθ,ϕ) in data) + sum(abs2,pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@time ellh(;data=data, pars=rand(length(LMs)) + 1im*rand(length(LMs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let mηπ_lims = (1.5,1.54)\n",
    "    ldata = filter(x->1.5<x.mηπ<1.54, data)\n",
    "    Np = length(LMs)\n",
    "    fold(x) = x[1:Np] + 1im .* x[(Np+1):2Np]\n",
    "    init_pars = fold(2rand(2Np).-1) # get random starting parameters\n",
    "    init_pars .*= sqrt(length(data)/sum(abs2,init_pars)) # normalize\n",
    "    #\n",
    "    f(x) = ellh(;data=ldata, pars=x)\n",
    "    f′(x) = fold(ForwardDiff.gradient(p->f(fold(p)), vcat(real(x), imag(x))))\n",
    "    f′!(stor,x) = copyto!(stor,f′(x))\n",
    "    #\n",
    "    Optim.optimize(f, f′!, init_pars, BFGS(),\n",
    "                Optim.Options(show_trace = true))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if the constraint is fulfilled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(abs2,Optim.minimizer(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mass production:\n",
    " - create a settings file several fit attempts from random starting points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function fit_data!(settings)\n",
    "    ldata = settings[\"data\"]\n",
    "    _Natt = settings[\"Natt\"]\n",
    "    #\n",
    "    Np = length(LMs)\n",
    "    fold(x) = x[1:Np]+x[Np+1:2Np]\n",
    "    #\n",
    "    f(x) = ellh(;data=ldata, pars=x)\n",
    "    f′(x) = fold(ForwardDiff.gradient(p->f(fold(p)), vcat(real(x), imag(x))))\n",
    "    f′!(stor,x) = copyto!(stor,f′(x))\n",
    "    # \n",
    "    frs = [let\n",
    "        init_pars = fold(2rand(2Np).-1) # get random starting parameters\n",
    "        init_pars .*= sqrt(length(data)/sum(abs2,init_pars)) # normalize\n",
    "        Optim.optimize(f, f′!, init_pars, BFGS(),\n",
    "                    Optim.Options(show_trace = settings[\"show_trace\"]))\n",
    "    end for e in 1:_Natt]\n",
    "    settings[\"fit_results\"] = frs\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = Dict(\n",
    "    \"data\"=>filter(x->1.5<x.mηπ<1.54, data),\n",
    "    \"Natt\"=>100,\n",
    "    \"show_trace\"=>false)\n",
    "fit_data!(settings);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of the fit results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Optim.minimizer(settings[\"fit_results\"][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tfr = Table(\n",
    "    [(st = Optim.converged(fr), min = Optim.minimum(fr), pars = Optim.minimizer(fr))\n",
    "        for fr in settings[\"fit_results\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "converged = filter(x->x.st, tfr);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram(converged.min, bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestpars = converged[findmin(converged.min)[2]].pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let pars = bestpars\n",
    "    cosθv = range(-1,1,length=100)\n",
    "    ϕv = range(-π,π,length=100)\n",
    "    calv = [model(cosθ,ϕ; pars=pars) for ϕ in ϕv, cosθ in cosθv]\n",
    "    heatmap(cosθv, ϕv, abs2.(calv), lab=\"\",\n",
    "        color=cgrad(:viridis, scale=:exp))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram2d(settings[\"data\"].cosθ, settings[\"data\"].ϕ, bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let\n",
    "    Nd = length(settings[\"data\"])\n",
    "    Nb = 50\n",
    "    stephist(settings[\"data\"].cosθ, bins=Nb, lc=:black, lab=\"data\")\n",
    "    Int(cosθ) = quadgk(ϕ->abs2(model(cosθ,ϕ; pars=bestpars)),-π,π)[1] / Nb * (1+1)\n",
    "    plot!(Int, -1, 1, lw=3, lab=\"model\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": "33ef23ec9ab74d97a8e76aa10280cf0c",
   "lastKernelId": "121be517-36e0-4efc-af46-53d086db378b"
  },
  "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
}