From e140b902b8e3b58a6c99c8226a5be64b9335aba2 Mon Sep 17 00:00:00 2001 From: JulienDoerner Date: Tue, 9 Jan 2024 08:10:45 +0100 Subject: [PATCH] remove photon propagation with DINT and ELECA from examples --- .../secondaries/photons.ipynb | 288 ------------------ 1 file changed, 288 deletions(-) diff --git a/doc/pages/example_notebooks/secondaries/photons.ipynb b/doc/pages/example_notebooks/secondaries/photons.ipynb index 59814e3e4..529d4b85a 100644 --- a/doc/pages/example_notebooks/secondaries/photons.ipynb +++ b/doc/pages/example_notebooks/secondaries/photons.ipynb @@ -123,294 +123,6 @@ "ylabel(\"Number of Particles\")\n", "show()" ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Photon Propagation outside of CRPropa with EleCa and DINT\n", - "\n", - "There are two main ways to propagate electromagnetic particles (EM particles: photons, electrons, positrons) in CRPropa:\n", - "\n", - "1) propagation as part of the CRPropa simulation chain\n", - "\n", - "2) propagation outside of the CRPropa simulation chain\n", - "\n", - "The following describes option 2, for which CRPropa provides three functions.\n", - "EM particles can either be propagated individually using the external EleCa code (suitable for high energies), or their spectra can be propagated with the transport code DINT (suitable for low energies).\n", - "Alternatively, a combined option is available that processes high energy photons with Eleca and then calculates the resulting spectra with DINT down to low energies.\n", - "\n", - "All three functions take as input a plain-text file with EM particles in the format given in the \"Photons from Proton Propagation\" example below.\n", - "In the following examples the input file \"photon_monoenergetic_source.dat\" contains 1000 photons with E = 50 EeV from a photon source at 4 Mpc distance.\n", - "\n", - "The last example \"Photons from Proton Propagation\" shows how to obtain secondary EM particles from a simulation of hadronic cosmic rays.\n", - "\n", - "Note that the differing results in EleCa (and correspondingly the high energy part of the combined option) are due to an incorrect sampling of the background photon energies in EleCa. The EleCa support will be removed in the near future.\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Photons from Proton Propagation\n", - "\n", - "The generation of photons has to be enabled for the individual energy-loss processes in the module chain. Also, separate photon output can be added:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "crpropa::ModuleList: Number of Threads: 8\n" - ] - } - ], - "source": [ - "from crpropa import *\n", - "\n", - "# source setup\n", - "source = Source()\n", - "source.add(SourceParticleType(nucleusId(1, 1)))\n", - "source.add(SourcePowerLawSpectrum(10 * EeV, 100 * EeV, -2))\n", - "source.add(SourceUniform1D(3 * Mpc, 100.00001 * Mpc))\n", - "\n", - "# setup module list for proton propagation\n", - "m = ModuleList()\n", - "m.add(SimplePropagation(0, 10 * Mpc))\n", - "m.add(MinimumEnergy(1 * EeV))\n", - "\n", - "# observer\n", - "obs1 = Observer() # proton output\n", - "obs1.add( Observer1D() )\n", - "obs1.add( ObserverPhotonVeto() ) # we don't want photons here\n", - "obs1.onDetection( TextOutput('proton_output.txt', Output.Event1D) )\n", - "m.add(obs1)\n", - "obs2 = Observer() # photon output\n", - "obs2.add( ObserverDetectAll() ) # stores the photons at creation without propagating them\n", - "obs2.add( ObserverNucleusVeto() ) # we don't want hadrons here\n", - "out2 = TextOutput('photon_output.txt', Output.Event1D)\n", - "out2.enable(Output.CreatedIdColumn) # enables the necessary columns to be compatible with the DINT and EleCa propagation\n", - "out2.enable(Output.CreatedEnergyColumn)\n", - "out2.enable(Output.CreatedPositionColumn)\n", - "out2.disable(Output.CandidateTagColumn)\n", - "obs2.onDetection( out2 )\n", - "m.add(obs2)\n", - "\n", - "# secondary electrons are disabled here\n", - "m.add(ElectronPairProduction(CMB(), False))\n", - "# enable secondary photons\n", - "m.add(PhotoPionProduction(CMB(), True))\n", - "\n", - "# run simulation\n", - "m.run(source, 10000, True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The file 'photon_output.txt' will contain approximately 300 photons and can be processed as the photon example below." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Propagation with EleCa\n" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2023-02-02 13:46:17 [WARNING] EleCa propagation is deprecated and is no longer supported. Please use the EM* (EMPairProduction, EMInverseComptonScattering, ...) modules instead.\n", - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Run ElecaPropagation\n", - " Started Thu Feb 2 13:46:17 2023 : [\u001b[1;32m Finished \u001b[0m] 100% Needed: 00:00:01 - Finished at Thu Feb 2 13:46:18 2023\n", - "\r" - ] - } - ], - "source": [ - "import crpropa\n", - "\n", - "# Signature: ElecaPropagation(inputfile, outputfile, showProgress=True, lowerEnergyThreshold=5*EeV, magneticFieldStrength=1*nG, background=\"ALL\")\n", - "crpropa.ElecaPropagation(\"photon_output.txt\", \"photons_eleca.dat\", True, 0.1*crpropa.EeV, 0.1*crpropa.nG)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Propagation with DINT\n" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2023-02-02 13:46:20 [WARNING] DINT propagation is deprecated and is no longer supported. Please use the EM* (EMPairProduction, EMInverseComptonScattering, ...) modules instead.\n", - "\n" - ] - } - ], - "source": [ - "import crpropa\n", - "\n", - "# Signature: DintPropagation(inputfile, outputfile, IRFlag=4, RadioFlag=4, magneticFieldStrength=1*nG, aCutcascade_Magfield=0)\n", - "crpropa.DintPropagation(\"photon_output.txt\", \"spectrum_dint.dat\", 4, 4, 0.1*crpropa.nG)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Combined Propagation" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Run EleCa propagation\n", - " Started Thu Feb 2 13:46:28 2023 : [\u001b[1;32m Finished \u001b[0m] 100% Needed: 00:00:00 - Finished at Thu Feb 2 13:46:28 2023\n", - "\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2023-02-02 13:46:28 [WARNING] EleCa+DINT propagation is deprecated and is no longer supported. Please use the EM* (EMPairProduction, EMInverseComptonScattering, ...) modules instead.\n", - "\n" - ] - } - ], - "source": [ - "import crpropa\n", - "\n", - "# Signature: DintElecaPropagation(inputfile, outputfile, showProgress=True, crossOverEnergy=0.5*EeV, magneticFieldStrength=1*nG, aCutcascade_Magfield=0)\n", - "crpropa.DintElecaPropagation(\"photon_output.txt\", \"spectrum_dint_eleca.dat\", True, 0.5*crpropa.EeV, 0.1*crpropa.nG)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### (Optional) Plotting of Results" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": { - "jupyter": { - "outputs_hidden": false - } - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAF3CAYAAACc3I0/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAouElEQVR4nO3de5xU9Znn8c9Di6C2NFFpFAkL2IrBC4W2KJtMaONlTYQxmRgvPZsRNXTUGPWVmFmcJKvJ6kQnycR4SZzOJCEX8cYaR1g1GTKAkYkKxPKKmhbIhGAAMWksjXLx2T+qCqqrq+iqrnPqVJ36vl+vflH1O6dPPT8aeHjO+V3M3REREQnLkKgDEBGReFOiERGRUCnRiIhIqJRoREQkVEo0IiISKiUaEREJ1V5RB1ANBx10kI8fPz7qMERE6sqqVatec/dRlV6nIRLN+PHjWblyZdRhiIjUFTP7XRDX0a0zEREJlRKNiIiESolGRERCFetnNGY2C5jV1tbW79j27dtZv349b7/9dvUDi7Hhw4czduxYhg4dGnUoIlIjrBEW1Wxvb/f8wQBr165l//3358ADD8TMIoosXtydLVu28MYbbzBhwoSowxGRCpnZKndvr/Q6DXvr7O2331aSCZiZceCBB6pKFJE+GjbRAEoyIdDvqYjka+hEE7WmpiYSicSurxtvvBGAjo6OAef9bN++nblz53L44Ydz3HHHMX36dB5++OFqhC0iUpZYDwaodfvssw/JZHJQ3/vlL3+ZV199leeee45hw4axceNGli1bFmyAIiIBUEVT437xi18wffp0jjvuOD7xiU+QSqV46623+N73vsett97KsGHDABg9ejTnnHMOAJdeeint7e0cddRRXHvttVGGLyKiigbgt1f9llQyFeg1mxPNHH7z4Xs85y9/+QuJRGLX+2uuuYZzzz131/vXXnuN66+/nsWLF7Pffvtx00038c///M989KMfZdy4cYwYMaLgdW+44QYOOOAAdu7cySmnnMIzzzzDscceG0i/RETKpUQToYFunT3++OO88MILvP/97wdg27ZtTJ8+fcDr3nvvvXR3d7Njxw5effVVXnjhBSUaEYmMEg0MWHlExd057bTTuOuuu/q0v/XWW/zXf/0XW7du7VfVrF27lm984xusWLGC97znPcyePVvDjUUkUnpGU8NOOukkli9fTk9PDwBvvvkmL7/8Mvvuuy8XX3wxV155Jdu2bQNg8+bN3HfffWzdupX99tuPlpYWNm7cqJFoIhK5WCcaM5tlZt29vb1Rh1JQ9hlN9mvu3Ll9jo8aNYp58+Zx/vnnc+yxxzJ9+nRefPFFAK6//npGjRrF5MmTOfroo5k5cyYjRoxgypQpTJ06lSOPPJLOzs5dt91ERKLSsEvQrF69mve9730RRRRv+r0ViQctQSMiInVBiUZEREKlRCMiIqFSohERkVAp0YiISKiUaEREJFRKNBHKbhNw1FFHMWXKFL75zW/y7rvvArB06VJmzpwJwLx58xgyZAjPPPPMru89+uijWbduHSeeeCKJRIJx48YxatSoXXNy1q1bF0WXRET60RI0Ecpd62zTpk10dnaydetWvvKVr/Q7d+zYsdxwww3cc889fdqfeOIJIJ2MVq5cyW233RZ63CIi5VBFUyNaW1vp7u7mtttuo9Ak2pkzZ/L888/z0ksvRRCdiMjgqaIBrroKBrn/WFGJBNx8c3nfM3HiRHbu3MmmTZv6HRsyZAh///d/zz/+4z/yox/9KJAYRUSqQRVNHens7OTxxx9n7dq1UYciIlIyVTSUX3mEZc2aNTQ1NdHa2srq1av7Hd9rr734/Oc/z0033RRBdCIig6OKpkZs3ryZSy65hMsvvxwzK3re7NmzWbx4MZs3b65idCIig6eKJkLZbQK2b9/OXnvtxSc/+Uk+97nP7fF79t57b6644gquvPLKKkUpIlKZWG8TYGazgFltbW1zfvvb3/Y5pqXsw6PfW5ESdXfD/PmFj3V2QldXdePJo20CSuDuC929q6WlJepQRET6mz+/8JDXZLJ4AqpDunUmIhKlRAKWLu3b1tERQSDhiXVFIyIi0VOiERGRUCnRiIhIqJRoREQkVEo0EfvjH//Ieeedx2GHHcbxxx/PRz7yEV5++eWKrjl79mwWLFjQr33lypVcccUVFV07a968eVx++eWBXEtE4k2jziLk7nzsYx/jggsu4O677wbg6aefZuPGjRxxxBGBf157ezvt7RUPiRcRKYsqmggtWbKEoUOHcskll+xqmzJlCh/4wAf4whe+wNFHH80xxxyzaw+apUuXMmPGDM466ywmTpzI3LlzufPOO5k2bRrHHHMMr7zyyq7rLF68mPb2do444ggWLVq06/uzm6ldd911XHTRRXR0dDBx4kRuueWWXd/705/+lGnTppFIJPj0pz/Nzp07AfjhD3/IEUccwbRp01i+fHnovz8ikqO7Oz3sudBXd3eUkQ1IFQ1Etk/Ac889x/HHH9+v/f777yeZTPL000/z2muvccIJJ/DBD34QSFc8q1ev5oADDmDixIl86lOf4sknn+Tb3/42t956KzdnPnPdunU8+eSTvPLKK5x88sn09PT0+5wXX3yRJUuW8MYbbzBp0iQuvfRSenp6uOeee1i+fDlDhw7lsssu48477+S0007j2muvZdWqVbS0tHDyySczderUSn+XRKRU2cmdiUTf9uy/XRGvIrAnSjQ16LHHHuP888+nqamJ0aNHM2PGDFasWMGIESM44YQTOOSQQwA47LDDOP300wE45phjWLJkya5rnHPOOQwZMoTDDz+ciRMn8uKLL/b7nDPPPJNhw4YxbNgwWltb2bhxI7/85S9ZtWoVJ5xwApBej621tZUnnniCjo4ORo0aBcC5555b8bMkESlTnU7uVKKByPYJOOqoowo+tN+TYcOG7Xo9ZMiQXe+HDBnCjh07dh3LXwG60IrQuddqampix44duDsXXHABX/va1/qc+8ADD5QVp4hIlp7RROhDH/oQ77zzDt0591efeeYZRo4cyT333MPOnTvZvHkzjz76KNOmTSvr2vfddx/vvvsur7zyCmvWrGHSpEklfd8pp5zCggULdu3y+frrr/O73/2OE088kWXLlrFlyxa2b9/OfffdV1Y8ItK4VNFEyMz42c9+xlVXXcVNN93E8OHDGT9+PDfffDOpVIopU6ZgZvzTP/0TBx98cMHbX8WMGzeOadOmsXXrVu644w6GDx9e0vdNnjyZ66+/ntNPP513332XoUOHcvvtt3PSSSdx3XXXMX36dEaOHEki/z6xiEgRsd4mIKu9vd1XrlzZp01L2YdHv7ciJco+Xyn23CW3vZxzA6JtAkREpC4o0YiISKiUaEREJFQNnWga4flUten3VETyNWyiGT58OFu2bNE/jAFyd7Zs2VLyCDcRaQwNO7x57NixrF+/ns2bN0cdSqwMHz6csWPHRh2GiNSQhk00Q4cOZcKECVGHISISezWfaMzso8CZwAjg++7+CzPbD/gOsA1Y6u53RhiiiMgedXen18TMd3MSWkfDmKpHVF2RPKMxsx+Y2SYzey6v/Qwze8nMesxsLoC7P+Duc4BLgHMzp/4NsCDT/tdVDV5EpEzZhZfzpVKwaWPVw6m6qCqaecBtwI+zDWbWBNwOnAasB1aY2YPu/kLmlC9ljgOMBZ7NvN5ZjYBFRCpRaOHl5MgIAolAJBWNuz8KvJ7XPA3ocfc17r4NuBs4y9JuAh52999kzl1POtlAA4+cExGpB7X0jOZQ4Pc579cDJwKfBU4FWsyszd3vAO4HbjOzM4GFhS5mZl1AF6QXmBQRqQUbujewcX76ftm7qfQNmac6ngJgdOdoxnTF74lNLSWagtz9FuCWvLY3gQsH+L5uoBvSi2qGFqCISBk2zt9IKpmiOdHcpz2VTAHEMtHU0m2nPwDvzXk/NtMmIhIrzYlmpi6dSlNzE03NTUxdOrVf4omTWko0K4DDzWyCme0NnAc8GHFMIiJSoaiGN98F/BqYZGbrzexid98BXA78HFgN3Ovuz1f4ObPMrLu3t7fyoEVEZFAieUbj7ucXaX8IeCjAz1kILGxvb58T1DVFRKD4JEyAzk7o6qpuPLWslm6dhUfrmYlIwIpNwkwmiyegRlXzo84C8Xr+lB0RkcoVmoSZ3VlZdot1RZN9RrNjx46oQxERaVixrmiyz2imDN9fz2hEpKpyJ2ae9uSvOGPbfbwx8i9MSO2kqbkJOpppSyXpaU5EG2gVxLqiydqugkZEqiw7MRPgjG33ccTO9PKMTc1N7D16bwB6mhMsbu2MLMZqiXVFIyISpezEzOTIJl4hQeLPS/scv6oj/evVVY+suhqiohERkejEuqIxs1nArCOH7Bd1KCJSY4rNg9EcmODFOtFkBwMcZ3vNKTjmUH+iRBpWdh5MIrG7LTsvRv8sBCvWiSbr7aZ9+jfqT5RIw8ufB6M5MOFoiETz+30maVaViEhENBhARERC1RAVjYg0rmIP/fOfzwRl5oZuTt00n8NIz6FppImZxcQ60WRHnQ0Zcly/O2U3J6F1NMRvLzsRyVXooT+k33eWOFeyULIqlqhO3TSftlQSb27b1ZadmFng9IYQ60STHXW2//79twlIZf6zoUQjEn+FFr8sR6FktadE1dOcwBPfAmDq0qkNMzGzmFgnmqxJBcYCJEdGEYmI1KtKk1Uja4hEIyJSqmSy/6DUsJ7nNAolGhGRjGK3wsp5niP9KdGIiGR0dWkOdxg0j0ZEREIV64omO7y5ra1twHNFRAaSu5lZvtGdoxnTpXGshcS6onH3he7e1dLSEnUoIhIDuZuZ5UolU0UTkMS8ohERCVp2M7NcT3U8FVE09SHWFY2IiERPiUZEREKlW2ciIhHKThBNJdODlpo74rcWoxKNiEhEik0CjdtajEo0IiIRyZ0g+lRHD5BehDNuazHqGY2IiIQq1hWNJmyKyGDlTs48cMP/5T2bHmFCaidNzU3Q0Zw+qbNTa9aUINYVzUATNttSyfRTuNyv7u7qBSgiNSt3cuZ7Nj3CPqmXaGpuYu/Re6dPSCYLb90p/cS6otmTxa3pp3CJ3MZkMv2r/ociIuRMzuxoBo5n/9wNafL3EpCiGjbRLBrTxaIxXX03MtIfHBGRwDVsohERqWWpVN//+9bz3JqGTjT5O+nV8w9SROKjdXT/tnqeW9OwiabQRKl6/kGKSHyMOST9lXtrv57n1jRsoim0k149/yBFRGpVrIc3i4hI9Bq2ohERgeK7ZqaSKZoTzRFEFD+qaESkoRXbNbM50czozgJP5aVssa5otASNSOPo7i48UT+ZhERiz99baNfMSqRScFXOsv+lxBBnsa5oBlqCRkTiY/783Yt75Eokii/HH4bW0dCcd8et2jHUmlhXNCLSWBKJvkOCo5Admvw9di/73+hiXdGIiEj0lGhERCRUSjQiIhIqJRoREQmVBgOISN0pNJS5lCHE2jUzGqpoRKTuFBrKXMoQYu2aGQ1VNCJSlwY7lFm7ZlafKhoREQmVEo2IiIRKiUZEREKlRCMiIqFSohERkVBp1JmI1KxKlv6X2hHrisbMZplZd29vb9ShiMgg1MrS/1KZWFc07r4QWNje3j4n6lhEZHBqYel/qUysKxoREYmeEo2IiIRKiUZEREKlRCMiIqEqaTCAmR1QwmnvuvufKwtHRETiptRRZxsyX7aHc5qAcRVHJCINabB7zEjtKzXRrHb3qXs6wcyeCiAeEWlQ2TkzuYlF82XiodREMz2gc0REitKcmXgqaTCAu79d7JiZHTzQOSIi0riCGHX2/QCuISIiMVXxEjTufmYQgYhIY9BCmY2nrERjZv+7ULu7fzWYcEQk7go99Ac9+I+zciuaN3NeDwdmAquDC0dEGoEe+jeWshKNu38z972ZfQP4eaARiYhIrFQ6GGBfYGwQgYiISDyV+4zmWcAzb5uAUUCsns+0pZLQ0dH/QGcndHVVOxwRkbpX7jOamTmvdwAb3X1HgPFEanFrJ6kUNCf7trelkmx9FcYo0YiUTKPLJKvcZzS/y28zs4Pd/Y/BhRSdEVd38aX5/ZPJ9Y910LwRxkQQk0i90ugyyQpiK+fvA7GYS9PVVfjuWHJk1UMRiQWNLhPQhE0RCYBWXpY90cZnIlKx7G2yXLpFJlllVzRm9h7gcNITNgFw90eDDCrv8yYCXwRa3P3sTNtk4DpgC/BLd18Q1ueLSGl0m0yKKauiMbNPAY+SnqT5lcyv15X7oWb2AzPbZGbP5bWfYWYvmVmPmc0FcPc17n5x3iU+DNzq7pcCf1fu54uISPWUe+vsSuAE4HfufjIwFfjzID53HnBGboOZNQG3k04ik4HzM5VLIT8BzjOzrwMHDuLzRUSkSspNNG9n950xs2Hu/iIwqdwPzdxqez2veRrQk6lgtgF3A2cV+f5N7v4ZYC7wWrmfLyIi1VNuollvZiOBB4B/N7N/A/rNrRmkQ4Hf534WcKiZHWhmdwBTzewaADMbb2bdwI+Brxe6mJl1mdlKM1u5efPmgEIUEZFylTth82OZl9eZ2RKgBXgk8Kj6fuYW4JK8tnXAHqfpu3s30A3Q3t7uezpXRETCU1JFY2a/yW9z92Xu/mDmNlfBc8r0B+C9Oe/HZtpERKSOlVrRvM/MntnDcSNd3VRiBXC4mU0gnWDOAzQKX0SkzpWaaI4s4ZydpX6omd0FdAAHmdl64Fp3/76ZXU56yHQT8AN3f77Uaxb5nFnArLa2tkouIyIiFSgp0RRaTLMS7n5+kfaHgIcC/JyFwML29vY5QV1TRETKoyVoREQkVEGs3iwiDUJ7zMhgxLqiMbNZZtbd29sbdSgisVBo8UzQApqyZ7GuaPSMRiR4oS2eWaxc0jbqdS/WFY2I1JFC5VIyWTj5SF0pqaIxs6MqHWosIjKg/HKpoyOiQCRIpVY0P8m+yGwVQM77fQONSEREYqXURGM5ry/LO/argGIREZEYKjXR5C5KaXnHavY5j0adiYhEr9RRZweb2WzgafonmppdGVmjzkQGr9AgMM2XkcEoNdFcBxwPXAiMNbMXgNXAi8BB4YQmIlHKDgLLTSyaLyODUWqiWQ8scvcNAGY2FjgGOBZ4NKTYRCRioc2ZkYZSaqL5GPBVMxtNuop5GkiSXgDzG+GEJiIicVDSg3x3n+Pu7cB3gZeBNcDJwBMEt5Vz4DQYQEQkeuWOGDvX3T/j7t9x94uBvwIeCyGuQLj7QnfvammpdE82EREZrHITzVYzOz77xt1XAUcEG5KIiMRJuYtqXgzcb2YrgFWkBwRsDzwqEakaLf0vYSuronH3l4HjgIeB0aSHOH8khLhEpEq09L+ErdRFNacDj3vaNuDezJeIxICGMUuYSq1o/g5YZWZ3m9lsMzs4zKBERCQ+Sqpo3P1SADM7EvgwMM/MWoAlwCPAcnffGVqUIiJSt0qqaMzsvwG4+4vu/i13PwP4EOmhzZ8gPZ+m5mgejYhI9EoddfawmbWSXhXgGeDZzK+PuftDYQVXKS2qKSISvVJvnU02s2HAZHavcXYWcKyZvePuE0KMUUQCohWZJQolD29293fc/SngZ6Rvlf0ReJv0umciUgcKDWXWMGYJW6nDmycBZwIzgVHAvwN3Al2Z4c4iUic0lFmqrdRnNKuBp4CbgH9z93fCC0lEROKk1ERzKXA08BngNjPbQnpAwLPAs+7+QDjhiYhIvSt1MMC/5L7P2/js48ADgUdWY1Ip6Ojo397ZCV1dVQ9HRKRulPqM5iNAMrvDpruvJ73r5sMhxlYxM5sFzGpra6voOq2jC7dnH6oq0YiIFFfpDptPAy/U6qoAQc2jGXMIjNmYZCkdfdqTwOINnYAyjcRQsWWdVcZLmUq9dTYHwMz+ATiU3TtsdgOvA2PDCrAmFBn72ZZKZl7pL53UlkCW/s+Ohc79BpXxMgjl7kdzrrtPyb4xs+8AXwg2pBrU1VXwL1bPyI7qxyJSgkI5AgYxZyZ/LHShB5UiAyg30Ww1s+MzO2vi7qvMTDtsitQgzZeRWqEdNkVEJFTaYVNEREI1YEVjZqcB5wC3u3sSmO3u3WiHTRERKUEpt84uIr0ywJfM7AAgEWpEIlJdgQxREymulETzhrv/GbjazG4ETgg3JBEpVWjDmEHLOktgSkk0/y/7wt3nmtlnQ4xHRMoQ2jBmkQANmGjc/d/M7Ch3fz7z/tbwwwpGUEvQiNQy5QipdaWOOvtJ9oWZfSr3gJntG2hEAXL3he7e1dLSEnUoIiINq9REYzmvL8s79quAYhERkRgqNdF4zmvLO1bWXBwREWkspa4McLCZzSa9WnN+ovH+p4uIiKSVmmiuA44HLgTGmtkLpFcFeBE4KJzQRCRXoaHMmuoi9aDUbQK6c9/n7bD5aAhxiUieQkOZyx7GrGwlESh3UU2gfnbYFImbiocyB5KtRMozqEQjInVME2+kyjRiTEREQqVEIyIioVKiERGRUCnRiIhIqDQYQKTGaHsYiRtVNCI1JjsCOZ9GIUu9UkUjUoMqHoGsskhqSKwrGjObZWbdvb29UYciUl0qi6SGxLqicfeFwML29vY5UcciUnWamCk1ItYVjYiIRE+JRkREQhXrW2citU6LKUsjUEUjEqFCz+z1vF7iRhWNSMT0zF7iThWNiIiESolGRERCpVtnFUqloKOjb1tnJ3R1RRKONBqtACB1QBVNBVpHQ3Nz37ZksvDfe5FQaAUAqQOqaCow5hAYszHJUjp2tSWBxRs6AZU0sluohYdGE0iNU6KpRIH/MbalkplXSjSyW7bwyE8qKjykESjRVKKrq9/DmJ6RHdHEIjVPhYc0Kj2jERGRUCnRiIhIqJRoREQkVEo0IiISKg0GEAmQ5k+K9KdEIxKgQIYxF8tWoGUnpC4p0YgErOJhzMWyVXYFACUaqTNKNCK1qFC26uhIJ5vcxfV0T07qgBKNSL0odO9NSwtIHaj5RGNmE4EvAi3ufnambRxwC/A68LK73xhhiCKDV85ezgVWohCpB5EMbzazH5jZJjN7Lq/9DDN7ycx6zGwugLuvcfeL8y5xDLDA3S8CplYpbJHgaS9naQBRVTTzgNuAH2cbzKwJuB04DVgPrDCzB939hQLf/ziwwMwuAn4SfrgiIdIiaBJzkSQad3/UzMbnNU8Detx9DYCZ3Q2cBRRKNBcC12auswD4YZjxihRSzl0vkUZWS89oDgV+n/N+PXCimR0I3ABMNbNr3P1rwCPAdWbWCawrdDEz6yKzVv+4cePCjFsaVKFRyEXvemkmpzSwWko0Bbn7FuCSvLbngLMH+L5uoBugvb3dQwtQGlrJd720IY00sFpKNH8A3pvzfmymTSQe9CxGGlQtLaq5AjjczCaY2d7AecCDEcckIiIVimp4813Ar4FJZrbezC529x3A5cDPgdXAve7+fIWfM8vMunt7eysPWkREBiWqUWfnF2l/CHgowM9ZCCxsb2+fE9Q1RUSkPLX0jEakJpU1YEyjy0T6qaVnNCI1qdDkfSgyYKysk0UaQ6wrGjObBcxqa2uLOhSpc2UNGNPoMpE+Yl3RuPtCd+9qaWmJOhQRkYYV64pGJFRag0akJEo0IUil+u5NlaVdeGOmrDVoRBqXEk3AWkdDWyrJzcmOPu2pFDy+QZkmdvQ8RmRASjQBG/P5TpgPibz21GNJmjdBZp1PqVG6GyYSvFgnmkhGnRXZBbFnZEf1YpBB090wkeDFOtFoZQAZDN0NEwlWrIc3i4hI9JRoREQkVLG+dSYSCK1fJlIRJRqRHDM3dHPqpvnQkdO4bFn61xkz+p6sUQIiJYl1otFaZ1JMsSLl+jXzaSNJnwHqM2Zotq1IBWKdaDTqTIopNIwZoLkZto5O0KxhZyKBiXWiESlm5oZubmZ+v4m1kIRD+reKyOBp1Jk0pFM3zactlex/QM9dRAKnikYaVk9zgoRukYmEThWNiIiESolGRERCFetbZxreLACP/s9uRizqO5Z5Qm+StS2JaAISaTCxrmi0lbMAjFg0nwm9yT5ta1sSbJ2ph/4i1RDrikYka21LgsSfl0YdhkhDUqKR+Cgy3b8tlaSnOVH9eEQEiPmtM2kw2en+eXqaEyxu1W0ykaioopF4KbBr2VUd6V+vrnYsIgIo0VTVhN4kybwtnbfO7OSDP9VijSISX0o0VbJ1ZidrF/Vtm9CbzLQp0QRhw6uwaePuCiZL28aIRCvWiaaW5tGkq5a+CSW/upHKbNoIqVT/di1fJhKtWCcabRPQeJqb+z2iEZGIxTrRSExpGLNIXdHwZqk/GsYsUldU0Uh90jBmkbqhikZEREKlikbqjoYxi9QXVTRSdzSMWaS+qKKR2lZghFl2dJmGMYvUB1U0UtsKjDDT6DKR+qKKRmpf3ggzjS4TqS+xrmjMbJaZdff29kYdiohIw4p1otFWzvVvw6vpO2cdHbu/CszVFJEaFutEI/Wv0AgzjS4TqS96RhOxVCr9v/R8nZ3Q1Ui7BwywfplGmInULyWaCLWOTv9DenOyo097KgWPb2iwTJMdXZY34zI7wixR6HtEpC4o0URozOc7YT79/hFNPZakeRM03IZoWr9MJJaUaKLU1VWwaunRhmgiEiMaDCAiIqFSRSM1QQtlisSXKhqpCVooUyS+VNFI9WmhTJGGoopGqk8LZYo0FFU0Eg0tlCnSMJRopOoKPfjXQ3+R+NKtM6k6rV8m0lhU0dSoCb1JknkTN7fO7OSDP43HagHNzf0WARCRmIp1ojGzWcCstra2qEMpy9aZnaxd1LdtQm8y01ZHiWaAhTJFpDHEOtG4+0JgYXt7+5yoYylHumrpm1Dyq5u6oIUyRYSYJxqpAVooU6ThKdFIaLSsjIiARp1JiLSsjIiAKhoJipaVEZEiVNFIMLSsjIgUoYpGArHhVdhEgqtYuqstCSTG6KG/SKNToqkjhSZxQm1M5NRsfxEpRommThSaxAm1NZFTs/1FpBAlmjpRaBInRDCRU7P9RaRMSjRSlg3fnM+IV/onlRQJHtdsfxEpQIlGyrJpI6whwZcSS/sd0/MYESlEiSYGQhkkMMAtMj2LEZFSaR5Nnds6s5O1LYl+7RN6k4xY1D9RlKzAvBjQ3BgRKZ8qmjq3p0ECbakkdHT0PdDZCV155xeoXt55MsnqvfvOiwHNjRGR8inRxNTi1k5SKWhO7m5L9C6DZcv63xJbtiz964wZu5pW753gh9v6Vy6aGyMi5VKiiakRV3fxpfl9K5cjlnXTyXxGJvNObpnB4tZOFuVURkkgMU3zYkSkcko0MdXVVegOWRfXzS9tcIAqFxEJihJNAymUfEREwqZRZyIiEiolGhERCZUSjYiIhKrmn9GY2UTgi0CLu5+dafsr4G9Jxz/Z3f97hCGKiMgeRFLRmNkPzGyTmT2X136Gmb1kZj1mNhfA3de4+8W557n7r9z9EmAR8KPqRS4iIuWK6tbZPOCM3AYzawJuBz4MTAbON7PJA1ynE6hgnRUREQlbJInG3R8FXs9rngb0ZCqYbcDdwFnFrmFm44Bed38jvEhFRKRStfSM5lDg9znv1wMnmtmBwA3AVDO7xt2/ljl+MfDDYhczs9xFwN7Jv003CC1Ab4XnFTo2UFsprw8CXishtj2ptH+ltu/pfVj9C+tnV6i93P7Vws+u2LHibWaFzgmmf7uvvfszzILtn+W09f283Lb+P8vdcezun9nu/vW51HEF2kqOLfO6aGz9r9H/9yiIf1smDRh9Kdw9ki9gPPBczvuzgX/Nef9J4LaAPmtlANforvS8QscGaivxdeT9K7V9T+/D6l9YP7sg+lcLP7tix0ppi3P/SulrXPoX5r8t7l5Tw5v/ALw35/3YTFutWBjAeYWODdRWyusgVNq/Utv39D6s/oX1syvUHqf+ldIW5/6V2tdK1UL/wvy3Bctkraozs/HAInc/OvN+L+Bl4BTSCWYF0OnuzwfwWSvdvb3S69Qq9a9+xblvoP7Vu6D6F9Xw5ruAXwOTzGy9mV3s7juAy4GfA6uBe4NIMhndAV2nVql/9SvOfQP1r94F0r/IKhoREWkMtfSMRkREYkiJRkREQqVEIyIioWqYRGNmE83s+2a2IK99PzNbaWYzo4otCIX6Z2YdZvYrM7vDzDqii65yRfo3xMxuMLNbzeyCKOOrVJH+/VXmZ/evZvafUcZXqSL9G2dmD2TWPpwbZXyVKNK3yWZ2r5l918zOjjK+SpnZR83se2Z2j5mdnmnbz8x+lGn/24GuUdeJptLFOTP+F3BvNeItVwD9cyAFDCe90kJNCaB/Z5Geb7WdGPbPa3zx2AB+fscAC9z9ImBqlcIuSQB9+zBwq7tfCvxdlcIuWZn9e8Dd5wCXAOdmTv0b0j+7OcBfD/iBQcz6jOoL+CDpdR5yVxhoAl4BJgJ7A0+T3koge3xBzuvTgPOA2cDMqPsTQv+GZH4dDdwZdX9C6N9c4NP57bXyVWn/ctruBfaPuj8h/PwOBJYA/wFcGHV/Au5bK+lFgr8OLI+6PwH175vAcZnX1wCJzOv5A31eXVc0XvninB3ASaRXgZ5jZjX1+1Fp/9z93czLPwHDQgt0kAL4+a0n3TeAneFEOXgB9K+mF48NoH8XAte6+4eAM8OLtHwB/N3b5O6fIf2foUrXsgtcOf2ztJuAh939N5lz15O+mwAl3BmrqX9YA1Jocc5DzexAM7uDzOKcAO7+RXe/ivRWA9/L+Ye5lpXcPzP7GzP7F+AnwG3VD3VQSu4fcD/wP8zsVuDRKsc5WOX0DwZYPLYGldO/R4ArMu3rqhvmoJTzd2+8mXUDPyZd1dSDgv0DPgucCpxtZpdkjt0PfNzMvksJS9bU0urNoXL3LaTvMRY6Nq+60QSvUP/c/X7SfyDqXpH+vUX6H+K6V+zPp7tfG0E4gSvy83uO9GK6da1I39axe/X4uubutwC35LW9SboiLUkcK5paX5yzUupffVP/6lec+wYh9i+OiWYFcLiZTTCzvUk/7H8w4piCpP7VN/WvfsW5bxBm/6Ie/VDhyIm7gFfZPbz14kz7R0ivBP0K8MWo41T/1D/1r76+4ty3KPqnRTVFRCRUcbx1JiIiNUSJRkREQqVEIyIioVKiERGRUCnRiIhIqJRoREQkVEo0IjnMbKeZJXO+amKflJy4xuzhnGvN7Gt5bQkzW515vcTMUmbWHna8Irk0j0Ykh5ml3L054Gvu5e47KrzGgHGZ2RHAI+4+MaftRuAtd/9q5v1S4Gp3X1lJPCLlUEUjUgIzW2dmXzGz35jZs2Z2ZKZ9v8wmUk+a2VNmdlamfbaZPWhm/wH80sz2tfSOiy+Y2c/M7Akzazezi8zs5pzPmWNm3yohntPN7NeZeO4zs2Z3fxn4k5mdmHPqOaRngYtERolGpK998m6dnZtz7DV3Pw74LnB1pu2LwH+4+zTgZODrZrZf5thxwNnuPgO4DPiTu08GvgwcnznnXmCWmQ3NvL8Q+MGeAjSzg4AvAadm4lkJfC5z+C7Sa1RhZicBr7v7b8v/bRAJTsNsEyBSor+4e6LIseyWC6tIb2ULcDrw12aWTTzDgXGZ1//u7tnNpT4AfBvSy+Ob2TOZ16lM1TMz8yxlqLs/O0CMJwGTgeVmBundEH+dOXYP8J9m9nnSCUfVjEROiUakdO9kft3J7r87Bnzc3V/KPTFz++rNEq/7r8A/AC9S2iZnRjqJnZ9/wN1/b2ZrgRnAx4HpJcYgEhrdOhOpzM+Bz1qmtDCzqUXOW076eQlmNhk4JnvA3Z8gvQ9IJ6VVII8D7zeztsz19ssMBMi6C/gWsMbd15fXHZHgKdGI9JX/jObGAc7/P8BQ4Bkzez7zvpDvAKPM7AXgeuB5oDfn+L3Acnf/00ABuvtmYDZwV+YW3K+BI3NOuQ84Ct02kxqh4c0iVWBmTaSfv7xtZocBi4FJ7r4tc3wR8C13/2WR7w9k2LWGN0sUVNGIVMe+wGNm9jTwM+Ayd99mZiPN7GXSgxAKJpmMrQNN2ByImS0BJpLe7EqkalTRiIhIqFTRiIhIqJRoREQkVEo0IiISKiUaEREJlRKNiIiESolGRERC9f8BMJRpbySANp4AAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "%matplotlib inline\n", - "import matplotlib.pyplot as plt\n", - "\n", - "plt.figure(figsize=(6,6))\n", - "\n", - "plt.loglog(clip_on=False)\n", - "plt.yscale(\"log\", nonpositive='clip')\n", - "plt.xlabel('Energy [eV]')\n", - "plt.ylabel ('$E^2 dN/dE$ [a.u.]')\n", - "\n", - "# Plot the EleCa spectrum\n", - "elecaPhotons = genfromtxt(\"photons_eleca.dat\")\n", - "binEdges = 10**arange(12, 24, .1)\n", - "logBinCenters = log10(binEdges[:-1]) + 0.5 * (log10(binEdges[1:]) - log10(binEdges[:-1]))\n", - "binWidths = (binEdges[1:] - binEdges[:-1])\n", - "data = histogram(elecaPhotons[:,1] * 1E18, bins=binEdges)\n", - "J = data[0] / binWidths\n", - "E = 10**logBinCenters\n", - "step(E, J * E**2, c='m', label='EleCa')\n", - "\n", - "#Plot the DINT spectrum\n", - "data = genfromtxt(\"spectrum_dint.dat\", names=True)\n", - "lE = data['logE']\n", - "E = 10**lE\n", - "dE = 10**(lE + 0.05) - 10**(lE - 0.05)\n", - "J = data['photons'] / dE\n", - "step(E, J * E**2 , c='b', where='mid', label='DINT')\n", - "\n", - "#Plot the combined DINT+EleCa spectrum\n", - "data = genfromtxt(\"spectrum_dint_eleca.dat\", names=True)\n", - "lE = data['logE']\n", - "E = 10**lE\n", - "dE = 10**(lE + 0.05) - 10**(lE - 0.05)\n", - "J = data['photons'] / dE\n", - "step(E, J * E**2 , c='r', where='mid', label='Combined')\n", - "\n", - "# Nice limits\n", - "xlim(1e14, 1e20)\n", - "ylim(bottom=1e17)\n", - "legend(loc='upper left')\n", - "show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": {