Touch up dipole examples, preliminary paper plots

Harvard-Neutrino · Feb 2, 2024 · 9195449 · 9195449
1 parent df5cccc
commit 9195449
Show file tree

Hide file tree

Showing 3 changed files with 153 additions and 193 deletions.
diff --git a/resources/Examples/Example2/DipolePortal_CCM.py b/resources/Examples/Example2/DipolePortal_CCM.py
@@ -1,5 +1,4 @@
 import os
-import sys
 import numpy as np
 
 import leptoninjector as LI
@@ -61,7 +60,7 @@
 primary_physical_distributions["flux_units"] = flux_units
 
 # direction distribution: cone from lower W target
-opening_angle = np.arctan(12 / 23.0)
+opening_angle = np.arctan(5 / 23.0)
 # slightly larger than CCM
 lower_target_origin = LI.math.Vector3D(0, 0, -0.241)
 detector_origin = LI.math.Vector3D(23, 0, -0.65)

diff --git a/resources/Examples/Example2/PaperPlots.ipynb b/resources/Examples/Example2/PaperPlots.ipynb
@@ -0,0 +1,152 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15a15725-d033-4b56-b3f4-55c2e949025e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import awkward as awk\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9500a367-6e8a-404e-8e64-fe7442a282ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f = \"output/CCM_Dipole_M2.35e-02_mu3.00e-07_example.parquet\"\n",
+    "data = awk.from_parquet(f)\n",
+    "POT = 2.25e22"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1e42dae5-0015-41c6-9f53-ddc1ce7c67ca",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Kinematic distributions\n",
+    "\n",
+    "# initial nu\n",
+    "nu_flag = data[\"primary_type\"]==\"ParticleType.NuMu\"\n",
+    "nu_momenta = np.squeeze(data[\"primary_momentum\"][nu_flag])\n",
+    "\n",
+    "# N\n",
+    "N_flag = data[\"primary_type\"]==\"ParticleType.N4\"\n",
+    "N_momenta = np.squeeze(data[\"primary_momentum\"][N_flag])\n",
+    "\n",
+    "# gamma\n",
+    "gamma_flag = data[\"secondary_types\"]=='ParticleType.Gamma'\n",
+    "gamma_momenta = data[\"secondary_momenta\"][gamma_flag]\n",
+    "# mask out entries that are not gamma\n",
+    "gamma_momenta = awk.mask(gamma_momenta, awk.num(gamma_momenta,axis=2)>0)\n",
+    "gamma_momenta = np.squeeze(gamma_momenta[~awk.is_none(gamma_momenta,axis=1)])\n",
+    "\n",
+    "# nu out\n",
+    "nuout_flag = data[\"secondary_types\"]=='ParticleType.NuLight'\n",
+    "nuout_momenta = data[\"secondary_momenta\"][nuout_flag]\n",
+    "# mask out entries that are not nuout\n",
+    "nuout_momenta = awk.mask(nuout_momenta, awk.num(nuout_momenta,axis=2)>0)\n",
+    "nuout_momenta = np.squeeze(nuout_momenta[~awk.is_none(nuout_momenta,axis=1)])\n",
+    "\n",
+    "kwargs = {\"bins\":np.linspace(0,0.031,50),\n",
+    "          \"weights\":data[\"event_weight\"]*POT,\n",
+    "          \"alpha\":0.5}\n",
+    "\n",
+    "# Energy\n",
+    "plt.hist(nu_momenta[:,0],**kwargs,label=r\"Initial $\\nu$\")\n",
+    "plt.hist(N_momenta[:,0],**kwargs,label=r\"Upscattered $\\mathcal{N}$\")\n",
+    "plt.hist(gamma_momenta[:,0],**kwargs,label=r\"Outgoing $\\gamma$\")\n",
+    "plt.hist(nuout_momenta[:,0],**kwargs,label=r\"Outgoing $\\nu$\")\n",
+    "plt.legend()\n",
+    "plt.semilogy()\n",
+    "plt.xlabel(\"Energy [GeV]\")\n",
+    "plt.ylabel(\"Event Rate\")\n",
+    "plt.show()\n",
+    "\n",
+    "# Angle\n",
+    "kwargs[\"bins\"] = np.linspace(-1,1,50)\n",
+    "def CosTheta(momenta,axis=1):\n",
+    "    return momenta[:,axis]/np.linalg.norm(momenta[:,1:],axis=-1)\n",
+    "    \n",
+    "plt.hist(CosTheta(nu_momenta),**kwargs,label=r\"Initial $\\nu$\")\n",
+    "plt.hist(CosTheta(N_momenta),**kwargs,label=r\"Upscattered $\\mathcal{N}$\")\n",
+    "plt.hist(CosTheta(gamma_momenta),**kwargs,label=r\"Outgoing $\\gamma$\")\n",
+    "plt.hist(CosTheta(nuout_momenta),**kwargs,label=r\"Outgoing $\\nu$\")\n",
+    "plt.legend()\n",
+    "plt.semilogy()\n",
+    "plt.xlabel(\"Energy [GeV]\")\n",
+    "plt.ylabel(\"Event Rate\")\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a167ce8f-adef-42e9-87cb-5448eaf6f04d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dec_flag = data[\"primary_type\"]=='ParticleType.N4'\n",
+    "gamma_flag = data[\"secondary_types\"]=='ParticleType.Gamma'\n",
+    "fid_flag = data[\"event_weight\"]<3e-24\n",
+    "\n",
+    "plt.hist2d(np.array(data[\"vertex\"][~dec_flag][:,0,0]),\n",
+    "           np.array(data[\"vertex\"][~dec_flag][:,0,1]),\n",
+    "           bins=(np.linspace(-5,30,50),\n",
+    "                 np.linspace(-2,2,50)))\n",
+    "plt.show()\n",
+    "\n",
+    "plt.hist2d(np.array(data[\"vertex\"][dec_flag])[:,0,0],\n",
+    "           np.array(data[\"vertex\"][dec_flag])[:,0,1],\n",
+    "           bins=(np.linspace(-5,30,50),\n",
+    "                 np.linspace(-2,2,50)))\n",
+    "plt.show()\n",
+    "\n",
+    "plt.hist2d(np.array(data[\"vertex\"][dec_flag])[fid_flag][:,0,0],\n",
+    "           np.array(data[\"vertex\"][dec_flag])[fid_flag][:,0,1],\n",
+    "           bins=(np.linspace(20,30,50),\n",
+    "                 np.linspace(-2,2,50)))\n",
+    "plt.show()\n",
+    "plt.hist(data[\"event_weight\"])\n",
+    "plt.loglog()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0cff9723-3d17-4de3-bbe4-3cda1e0d3a74",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}