1.0 Data exploration. Eta pi.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 UpROOT # reading ROOT files\n",
    "using Parameters # pack, unpack\n",
    "using TypedTables # \n",
    "using StaticArrays # fixed size arrays\n",
    "using Plots # plotting\n",
    "using LinearAlgebra # norm for a vector\n",
    "using DelimitedFiles # read and write tabled data\n",
    "using Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theme(:wong, bottom_margin=3Plots.PlotMeasures.mm) # colors of the plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data exploration: $\\eta\\pi$ from COMPASS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download(\"https://apearce.web.cern.ch/apearce/aman-intro/etapi_predicts.root\",\n",
    "    # joinpath(\"data\", \"etapi_predicts.root\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = TFile(joinpath(\"data\",\"etapi_predicts.root\")) # open the file\n",
    "t = f[\"pi03pic\"][:]; # get the tree from the disc to the memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "propertynames(t) # list of branches, the names of columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram(t.pPimZ, bins=100, xlab=\"pz(π⁻) (GeV)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mda, what should I expect??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Apgebra of four vectors\n",
    "Define the particle $(E,\\vec p)$:\n",
    " - invariant mass function, and basic transformations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@with_kw mutable struct Particle\n",
    "    E::Float64\n",
    "    p::MVector{3,Float64}\n",
    "end\n",
    "Particle(E,px,py,pz) = Particle(E=E, p=SVector(px,py,pz))\n",
    "Particle(px,py,pz; msq=error(\"mass needed\")) = Particle(E=sqrt(px^2+py^2+pz^2+msq), p=SVector(px,py,pz))\n",
    "\n",
    "# properties\n",
    "invmasssq(p::Particle) = p.E^2-sum(abs2,p.p)\n",
    "mass(p::Particle) = sqrt(invmasssq(p))\n",
    "# \n",
    "import Base:+,-,*,show\n",
    "+(p1::Particle, p2::Particle) = Particle(E=p1.E+p2.E, p=p1.p+p2.p)\n",
    "*(p::Particle, α::T where T<:Real) = Particle(E=p.E*α, p=p.p*α)\n",
    "*(α::T where T<:Real, p::Particle) = Particle(E=p.E*α, p=p.p*α)\n",
    "-(p1::Particle, p2::Particle) = p1+(-1)*p2\n",
    "# \n",
    "show(io::IO, p::Particle) = print(io, \"P($(p.E),$(p.p))\")\n",
    "\n",
    "# transformations\n",
    "function Ry!(p,θ)\n",
    "    sinθ, cosθ = sin(θ), cos(θ)\n",
    "    p.p[1], p.p[3] =  p.p[1]*cosθ+sinθ*p.p[3],\n",
    "                     -p.p[1]*sinθ+cosθ*p.p[3]\n",
    "end\n",
    "function Rz!(p,ϕ)\n",
    "    sinϕ, cosϕ = sin(ϕ), cos(ϕ)\n",
    "    p.p[1], p.p[2] = p.p[1]*cosϕ-sinϕ*p.p[2],\n",
    "                     p.p[1]*sinϕ+cosϕ*p.p[2]\n",
    "end\n",
    "function Bz!(p,γ)\n",
    "    _γ = abs(γ)\n",
    "    _βγ = sqrt(γ^2-1)*sign(γ)\n",
    "    p.p[3], p.E = _γ *p.p[3]+_βγ*p.E,\n",
    "                  _βγ*p.p[3]+ _γ*p.E\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise: $\\pi^0\\to \\gamma\\gamma$ decay:\n",
    " - get an example of two four vectors of $p_{\\gamma1}$ and $p_{\\gamma2}$ that are decays of the\n",
    "$$\n",
    "\\begin{align}\n",
    "p_{\\pi0} = (E,1,-1,30), \\text{ with } E = \\sqrt{30^2+1^2+1^2+m_{\\pi}^2}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const mπ0 = 0.135\n",
    "pπ0 = Particle(1,1,30; msq = mπ0^2)\n",
    "pγ1 = Particle(10,0,0,10); pγ2 = Particle(20,0,0,10) # is incorrect choise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@test invmasssq(pγ1) == 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@test invmasssq(pγ2) == 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@test invmasssq(pγ1+pγ2) == invmasssq(pπ0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@test sum((pγ1+pγ2).p .== pπ0.p) == 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Presentation of the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Converting the TTree to a table of `Particles`\n",
    "row by row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process(row) = \n",
    "   (pb = Particle(row.pBeamE,   row.pBeamX,   row.pBeamY,   row.pBeamZ),\n",
    "    pr = Particle(row.pRecoilE, row.pRecoilX, row.pRecoilY, row.pRecoilZ),\n",
    "    pπ = Particle(row.pPimE,    row.pPimX,    row.pPimY,    row.pPimZ),\n",
    "    pη = Particle(row.pEtaE,    row.pEtaX,    row.pEtaY,    row.pEtaZ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fv = Table(process.(t)); # create the four vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figs/etapi_d.PNG\" alt=\"etapi\" style=\"width: 300px;\"/>\n",
    "\n",
    "The scattering amplitude dependence on\n",
    "$$\n",
    "\\begin{align}\n",
    "s &= (p_\\eta + p_\\pi)^2\\\\\n",
    "t_1 &= (p_\\eta + p_{b})^2\\\\\n",
    "t &= (p_{t} - p_r)^2\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "The other two:\n",
    " - $s_0 = (p_b+p_t)^2$\n",
    " - $s_2 = (p_r+p_\\pi)^2$\n",
    " Compute them all:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(layout=grid(1,3), size=(900,250),\n",
    "    histogram(mass.(fv.pπ+fv.pη+fv.pr), bins=100, xlab=\"√s₀ (GeV)\"),\n",
    "    histogram(mass.(fv.pπ+fv.pη), xlab=\"m(ηπ) (GeV)\", bins=100),\n",
    "    histogram(-invmasssq.(fv.pb-fv.pη-fv.pπ), xlab=\"transf. momentum, -t (GeV²)\", bins=100, yscale=:log10)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram2d(mass.(fv.pπ+fv.pη),\n",
    "            invmasssq.(fv.pb-fv.pη), xlab=\"m(ηπ) (GeV)\", ylab=\"t₁=(pb-pη)² (GeV²)\",\n",
    "    color=cgrad(:viridis, scale=:exp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angles(p::Particle) = p.p[3]/norm(p.p), atan(p.p[2],p.p[1])\n",
    "boostzfactor(p::Particle) = p.E/mass(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Getting the right angles\n",
    "Introducing the Gottfried-Jackson frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figs/GJ.PNG\" alt=\"GJ angles\" style=\"width: 800px;\"/>\n",
    "\n",
    "\n",
    "Here:\n",
    "$$\n",
    "p_\\eta = p_2+p_3,\\quad p_1 = p_\\pi\n",
    "$$\n",
    "\n",
    "Get the angle $\\theta = \\theta_1$ and $\\phi = \\phi_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting recipe for the particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@recipe function f(p::Particle, proj=:zx)\n",
    "    arrow   := true\n",
    "    linewidth  :=  2\n",
    "    tip = (p.p[3],p.p[1])\n",
    "    tip = (proj == :xy) ? (p.p[1],p.p[2]) : tip\n",
    "    return [(0,0), tip]\n",
    "end\n",
    "function plot_particles(pb, pη, pπ, pr)\n",
    "    plot(layout=grid(1,2, widths=(0.7,0.3)), aspect_ratio=1, size=(850,280), frame=:origin)\n",
    "    plot!(sp=1, xlab=\"z\", ylab=\"x\");    plot!(sp=2, xlab=\"x\", ylab=\"y\")\n",
    "    plot!(sp=1, pb, :zx, lab=\"beam\");   plot!(sp=2, pb, :xy, lab=\"\")\n",
    "    plot!(sp=1, pη, :zx, lab=\"η\");      plot!(sp=2, pη, :xy, lab=\"\")\n",
    "    plot!(sp=1, pπ, :zx, lab=\"π\");      plot!(sp=2, pπ, :xy, lab=\"\")\n",
    "    plot!(sp=1, pr, :zx, lab=\"recoil\"); plot!(sp=2, pr, :xy, lab=\"\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@unpack pb, pη, pπ, pr = fv[6]\n",
    "plot_particles(pb, pη, pπ, pr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@unpack pb, pη, pπ, pr = fv[6]\n",
    "cosθ,ϕ = angles(pπ+pη)\n",
    "γ = boostzfactor(pπ+pη)\n",
    "#\n",
    "Rz!.([pb, pη, pπ, pr], -ϕ)\n",
    "Ry!.([pb, pη, pπ, pr], -acos(cosθ))\n",
    "Bz!.([pb, pη, pπ, pr], -γ)\n",
    "# #\n",
    "cosθ,ϕ = angles(pb)\n",
    "Rz!.([pb, pη, pπ, pr], -ϕ)\n",
    "Ry!.([pb, pη, pπ, pr], -acos(cosθ))\n",
    "# # \n",
    "cosθ,ϕ = angles(pr)\n",
    "Rz!.([pb, pη, pπ, pr], π-ϕ)\n",
    "# \n",
    "plot_particles(pb, pη, pπ, pr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function anglesGJ!(;pb, pη, pπ, pr)\n",
    "    pb, pη, pπ, pr\n",
    "    cosθ,ϕ = angles(pπ+pη)\n",
    "    γ = boostzfactor(pπ+pη)\n",
    "    #\n",
    "    Rz!.([pb, pη, pπ, pr], -ϕ)\n",
    "    Ry!.([pb, pη, pπ, pr], -acos(cosθ))\n",
    "    Bz!.([pb, pη, pπ, pr], -γ)\n",
    "    # #\n",
    "    cosθ,ϕ = angles(pb)\n",
    "    Rz!.([pb, pη, pπ, pr], -ϕ)\n",
    "    Ry!.([pb, pη, pπ, pr], -acos(cosθ))\n",
    "    # # \n",
    "    cosθ,ϕ = angles(pr)\n",
    "    Rz!.([pb, pη, pπ, pr], π-ϕ)\n",
    "    return angles(pη)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kinvars = Table([NamedTuple{(:cosθ,:ϕ,:mηπ)}([anglesGJ!(;v...)...,mass(v.pη+v.pπ)]) for v in fv])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(layout=grid(1,3), size=(900,250),\n",
    "    histogram(kinvars.cosθ, lab=\"\", xlab=\"cosθ\"),\n",
    "    histogram(kinvars.ϕ, lab=\"\", xlab=\"ϕ\"),\n",
    "    histogram(kinvars.mηπ, lab=\"\", xlab=\"mηπ\", bins=100)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore how the angular dependence changes with mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram2d(kinvars.mηπ, kinvars.cosθ, xlab=\"m(ηπ) (GeV)\", ylab=\"cosθ\",\n",
    "        color=cgrad(:viridis, scale=:exp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram2d(kinvars.cosθ, kinvars.ϕ, xlab=\"cosθ\", ylab=\"ϕ\", bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writedlm(joinpath(\"data\",\"data_metapi_costheta_phi.txt\"), [kinvars.mηπ kinvars.cosθ kinvars.ϕ])"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.9.3",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}