diff --git a/resources/Examples/Example1/DIS_ATLAS.py b/resources/Examples/Example1/DIS_ATLAS.py new file mode 100644 index 00000000..19d6718b --- /dev/null +++ b/resources/Examples/Example1/DIS_ATLAS.py @@ -0,0 +1,84 @@ +import os + +import leptoninjector as LI +from leptoninjector import _util +from leptoninjector.LIController import LIController + +resources_dir = _util.resource_package_dir() + +# Number of events to inject +events_to_inject = 10000 + +# Expeirment to run +experiment = "ATLAS" + +# Define the controller +controller = LIController(events_to_inject, experiment) + +# Particle to inject +primary_type = LI.dataclasses.Particle.ParticleType.NuMu + +cross_section_model = "CSMSDISSplines" + +xsfiledir = _util.get_cross_section_model_path(cross_section_model) + +# Cross Section Model +target_type = LI.dataclasses.Particle.ParticleType.Nucleon + +DIS_xs = LI.interactions.DISFromSpline( + os.path.join(xsfiledir, "dsdxdy_nu_CC_iso.fits"), + os.path.join(xsfiledir, "sigma_nu_CC_iso.fits"), + [primary_type], + [target_type], +) + +primary_xs = LI.interactions.InteractionCollection(primary_type, [DIS_xs]) +controller.SetInteractions(primary_xs) + +# Primary distributions +primary_injection_distributions = {} +primary_physical_distributions = {} + +# energy distribution +# HE SN flux from ATLAS paper +flux_file = os.path.join( + resources_dir, + "Fluxes", + "HE_SN_Flux_Tables", + "dN_dE_SNe_2n_D1_0_s20_t100d_NuMu_d10kpc.txt", +) +print(flux_file) +edist = LI.distributions.TabulatedFluxDistribution(100, 1e6, flux_file, True) #bool is whether flux is physical +#edist_inj = LI.distributions.PowerLaw(2, 100, int(1e6)) +primary_injection_distributions["energy"] = edist +primary_physical_distributions["energy"] = edist + +# we need this conversion to make sure the flux is in units of 1/m2 +flux_unit_conv = LI.distributions.NormalizationConstant((100 / 1)**2) +primary_physical_distributions["flux_norm"] = flux_unit_conv + +# direction distribution +# let's just inject upwards +injection_dir = LI.math.Vector3D(0, 0, 1) +injection_dir.normalize() +direction_distribution = LI.distributions.FixedDirection(injection_dir) +primary_injection_distributions["direction"] = direction_distribution +primary_physical_distributions["direction"] = direction_distribution + +# position distribution +muon_range_func = LI.distributions.LeptonDepthFunction() +position_distribution = LI.distributions.ColumnDepthPositionDistribution( + 600, 600.0, muon_range_func, set(controller.GetDetectorModelTargets()[0]) +) +primary_injection_distributions["position"] = position_distribution + +# SetProcesses +controller.SetProcesses( + primary_type, primary_injection_distributions, primary_physical_distributions +) + +controller.Initialize() + +events = controller.GenerateEvents() + +controller.SaveEvents("output/ATLAS_DIS") diff --git a/resources/Examples/Example1/DIS_DUNE.py b/resources/Examples/Example1/DIS_DUNE.py new file mode 100644 index 00000000..49594075 --- /dev/null +++ b/resources/Examples/Example1/DIS_DUNE.py @@ -0,0 +1,68 @@ +import os + +import leptoninjector as LI +from leptoninjector import _util +from leptoninjector.LIController import LIController + +resources_dir = _util.resource_package_dir() + +# Number of events to inject +events_to_inject = 10000 + +# Expeirment to run +experiment = "DUNEFD" + +# Define the controller +controller = LIController(events_to_inject, experiment) + +# Particle to inject +primary_type = LI.dataclasses.Particle.ParticleType.NuMu + +cross_section_model = "CSMSDISSplines" + +xsfiledir = _util.get_cross_section_model_path(cross_section_model) + +# Cross Section Model +target_type = LI.dataclasses.Particle.ParticleType.Nucleon + +DIS_xs = LI.interactions.DISFromSpline( + os.path.join(xsfiledir, "dsdxdy_nu_CC_iso.fits"), + os.path.join(xsfiledir, "sigma_nu_CC_iso.fits"), + [primary_type], + [target_type], +) + +primary_xs = LI.interactions.InteractionCollection(primary_type, [DIS_xs]) +controller.SetInteractions(primary_xs) + +# Primary distributions +primary_injection_distributions = {} +primary_physical_distributions = {} + +# energy distribution +edist = LI.distributions.PowerLaw(2, 1e3, 1e6) +primary_injection_distributions["energy"] = edist +primary_physical_distributions["energy"] = edist + +# direction distribution +direction_distribution = LI.distributions.IsotropicDirection() +primary_injection_distributions["direction"] = direction_distribution +primary_physical_distributions["direction"] = direction_distribution + +# position distribution +muon_range_func = LI.distributions.LeptonDepthFunction() +position_distribution = LI.distributions.ColumnDepthPositionDistribution( + 600, 600.0, muon_range_func, set(controller.GetDetectorModelTargets()[0]) +) +primary_injection_distributions["position"] = position_distribution + +# SetProcesses +controller.SetProcesses( + primary_type, primary_injection_distributions, primary_physical_distributions +) + +controller.Initialize() + +events = controller.GenerateEvents() + +controller.SaveEvents("output/DUNE_DIS") \ No newline at end of file diff --git a/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py index 87981247..2730fb3e 100644 --- a/resources/Examples/Example1/DIS_IceCube.py +++ b/resources/Examples/Example1/DIS_IceCube.py @@ -1,6 +1,4 @@ import os -import sys -import numpy as np import leptoninjector as LI from leptoninjector import _util @@ -9,7 +7,7 @@ resources_dir = _util.resource_package_dir() # Number of events to inject -events_to_inject = 1000 +events_to_inject = 10000 # Expeirment to run experiment = "IceCube" @@ -67,4 +65,4 @@ events = controller.GenerateEvents() -controller.SaveEvents("IceCube_DIS.hdf5") +controller.SaveEvents("output/IceCube_DIS") diff --git a/resources/Examples/Example1/PaperPlots.ipynb b/resources/Examples/Example1/PaperPlots.ipynb new file mode 100644 index 00000000..b38bd4ea --- /dev/null +++ b/resources/Examples/Example1/PaperPlots.ipynb @@ -0,0 +1,101 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "249e815d-3d83-4c81-bce7-bac786f95549", + "metadata": {}, + "outputs": [], + "source": [ + "import awkward as awk\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "id": "c80f20e3-357d-4157-a279-c8f491f0c50d", + "metadata": {}, + "outputs": [], + "source": [ + "f = \"output/IceCube_DIS.parquet\"\n", + "data = awk.from_parquet(f)" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "ef6c4d55-6687-425e-a521-de52d27de2f0", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG1CAYAAADwRl5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3c0lEQVR4nO3dB3RU1b7H8X9C71KkGxBFFJEiRVCQIkVYgKA8sUFAH1wFFeEiwsIHYgHliiKaiKgUuahYkYWXJhcFqQEElUiTSC/SQieSzFv/fdfkzpAyk2QmM7PP97PWWTAzZ2Z2ThLmx97/vXeUy+VyCQAAgIWiQ90AAACAYCHoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsVVAcLi0tTQ4ePCilSpWSqKioUDcHAAD4Qdc7PnPmjFStWlWio7Put3F80NGQc8011/hzTQEAQJjZt2+fVK9ePcvHHRt04uLizHH58uX0C1W6dOlQNwsAAPjh9OnTpqNCR2SyE+X0va70QpUpU0aSk5MJOgAAWPb5TTEyAACwFkEHAABYy7E1OgDgtBmmKSkpoW4G4LdChQpJgQIFJK8IOgBgOQ04SUlJJuwAkeSqq66SypUr52n5F4IOAFhM55scOnTI/M9YZ6hkt94IEE4/t+fPn5ejR4+a21WqVMn1axF0AMBiuoSGfmDoomrFixcPdXMAvxUrVsz8qWGnYsWKuR7GItoDgMVSU1PNn4ULFw51U4Acc4fzv/76S3KLoAMADsAWN3Dqzy1BBwAAWIugAwAArEXQAQAA1mLWFQA40JtLd+Tr+w3tcEPQ36NNmzbSsGFDmTx5ckDOy+35CC8EHQBAWOrXr5+cOnVK5s2b59f5X331lVlN11dAufI82I2gA8BR/OnJyI/eBwReuXLlAnoe7ECNDgAgImgPzdNPPy0jRowwYUW3BnjhhRe8Hn/mmWfSe4N++OEHeeutt8wUZT3++OOPDOepRYsWScuWLc12A+XLl5euXbvK77//7ne7unTpIrGxsem3ly9fLhUqVEhfwwihRdABAESMWbNmSYkSJWTdunUyceJEefHFF2Xp0qUZztOA06JFCxkwYIDZAkMP3QIjM+fOnZNhw4bJhg0bZNmyZWabjJ49e/q9N1i1atXkwIED6bdbt24tFy5ckLVr1+bhK0WgOHboKi4uzhwkbgCIHPXr15exY8eav9euXVveeecdE046dOjgdV6ZMmXMatC6sq72/GTnvvvu87o9ffp0ufrqqyUxMVHq1avnV9BZuXJl+m0NSrp9gXufJoSWY3t0Bg8ebH6IExISQt0UAEAOgo4n3ewxr4Fi586d8uCDD0qtWrWkdOnSUrNmTXP/3r17c9Wjs3nzZlNErT1KCD3H9ugAACLPlbOltPbG3yGmrHTr1k1q1Kgh77//vtn8VF9Pe3JSUlL8Djpnz56V06dPS8mSJWXo0KHy8MMP++xJQv4g6AAArKRDV77KE44fPy7bt283IadVq1bmvh9//DFH76NBR+3fv19mzJghhw8flm+++SYPLUcgEXQAAFbSISgtWtbZVtrTojO1tH7GU9myZc1Mq2nTpplhMB2uGjlyZK6Czt///nfZsWOHrFixwgyBITwQdADAgZywVtDw4cPNtO+6deuaWVBJSUnp9TduGnw+/fRTM21dh6vq1KkjU6ZMMVPQ/aVTyYsUKSJ79uwxU9rdwQfhIcrlcrnEwXRMVavzk5OTSeCAAzhtwcCLFy+aD/hrr71WihYtGurmAAH7+fX389uxs64AAID9CDoAAMBaBB0AAGAtgg4AALAWs64AIBcFy7YVLQO2okcHAABYi6ADAACsRdABAADWIugAAABrUYwMwB7LJ/hx0n350BAA4YKgA8Aaa3Yf931STH60BE6i+2I1bNhQJk+eHOqmIBMEHQBwIr96vwKo7ahcPW3fvn0yduxYWbRokRw7dszsMN6jRw8ZM2aM2XU8HMLIV199JYUKFQr46yIwqNEBAISl3bt3S5MmTWTnzp3yySefyK5du2Tq1KmybNkyadGihZw4cULCQbly5aRUqVKhbgayQNABYI350bt8HogcgwcPlsKFC8uSJUukdevWEhMTI507d5bvvvtODhw4IKNHjzbn1axZM0NPjfbevPDCC+bv/fr1kx9++EHeeustiYqKMscff/xhHjtz5ow8/PDDUqJECdNb9Oabb5ren2eeecY8funSJXn66aelYsWKZvfsli1bSkJCgtd7eZ7vvq3PGTFihAlBlStXTm+Lm6/3zcyPP/5oeo50R283/Tr069mzZ08errTdCDoAgLCjvTWLFy+WQYMGSbFixbwe0+CgIWHu3Lnicrl8vpYGHO0BGjBggBw6dMgc11xzjXls2LBhsmrVKpk/f74sXbpUVq5cKZs2bUp/roaVL7/8UmbNmmXuv/7666VTp04+e5P0fA0x69atk4kTJ8qLL75oXt/N1/tmZvPmzXLTTTeZwOX2008/SdmyZaVGjRo+r4NTUaMDAEHcKoJtInJHh6s0xOgHe2b0/pMnT8qff/7p87XKlCljeoaKFy9uQpJnr4oGko8//ljuuusuc9+MGTOkatWq5u/nzp2Td999V2bOnGl6ktT7779vgsmHH34ozz77bJbvWb9+fVNbpGrXri3vvPOOGXLr0KGDz/fNypYtW6RRo0YZwk+DBg18XgMno0cHABC2/OmxyUsN0F9//SXNmjXzCkV16tQxf//999/N43fccUf64zp0pOf/9ttv2b62Bh1POjx19OhRv943KxpqdEjOk/boXHkfvBF0AABhR4eItPYkq0Ch9+uQzdVXXy3R0dEZApEGiVC6chaWfi1paWm5fr3U1FT59ddfM/To6HAXQcfyoHPq1ClTla/f6Hr16pluRQDISvO903weCD2dOq7DPPHx8XLhwgWvxw4fPixz5syR3r17mwChYUfrbtxOnz4tSUlJXs/RoSsNC55q1aplAolncXFycrLs2PGfIcnrrrvOPE9raTwDlJ5ft27dXH9tvt43M9u3bzdFyJ7DW2vWrDFF2QQdy4OOTulbsWKF6dLToq/x48fL8eN+LBoGAAhrWteis560+Ff/ndc1dXQ9HQ1A1apVk1deecWc165dO5k9e7Yp6P3ll18kNjZWChQo4PVaOjNLPyN0lpKux6O9K/r5oedqrc3y5ctl69at8thjj5keIg1QWkz8xBNPmMf1fRMTE01B8/nz5815ueXrfTOjn3Hq7bffNvVLCxculL59+5r7UlJSct0WJ4j4YmT9YdYCM6W/ENp9GcwxXQChEb85nkvvMFrEu2HDBlPUe//995uZTlpMrAsG6n06dVuNGjXK9OB07drV1Lq89NJLGXp0hg8fbsKF9sRoD5E+ruHnjTfekMcff9w8t3Tp0maWlQYq98ymV1991YSiPn36mCJiHUHQ2WA6bJYXvt43s6CjgU/re2655RbzdYwbN84EsSlTppigh8xFuUKcCjSl/+Mf/5CNGzearsevv/7a/BB7iouLM+dod6VWl2ui9Szi0uErXWNBU66ep2sv+Eu7OPUXQ7sN9YcNQOQGnX0/LfF5Tve0632eszZmoARKqGdd6XCHfqhfe+21WX6I4r90ppX2Fk2aNClPvTaBfl8NOU2bNpWXX37ZUd+ui9n8/Pr7+R3yoSv95mp40TCTGV0nQdcb0PSuRVd6rn7D3dXr6qqrrjLT7vRi6HS9I0eOZPl+2uujF8fzAAA4k85a0lWXdYaVfsbo+jzqnnvuCav31c847clBBA5d6doE7vUJsure0zHR/v37m9u6/Pe3334r06dPl5EjR3qdW6lSJROEdJy2V69emb7ehAkTTHcfgAiTtDLULYClXn/9dVPsq4XHjRs3Np8hFSpUCJv31dEM/Q88QSdCg052tMBKh7R0/NVNi7Xat29vqs2VfvO1RkeLu7T7SofCdMwyK/pa2kPkpj067hUyAdjPn20gKuZLSxAOdLq2fs6E8/tqXRK1p5YGHa2M1+mA2lPjSW9v27bN/F339xg4cGB6EfJTTz2VbeotUqSIOQAAgP3COuj4Q4uS3dPuAAAAwqoYOTs6VqnTx68sLtbbnvuVAAAARFyPjrtASzdCc0851/UM9PaTTz6Zp9fWWV56XLlSJgAEEht/Ag4POmfPnpVdu/5bHKhTxHUoSheCiomJMYXDusiTLtKkw1STJ082U9Lds7ByS9fa0cM9Dx8AANgn5EFHV71s27Zt+m33jCgNNzNnzjR7mfz5558yZswYM8VO9/TQpbivLFAGAAAIu6DTpk0bn9PmdJgqr0NVAOAvfzf2DOQKygAcWIwMAAAQ0T06oUIxMhA+2LATker777835RcnT5402xFFqpkzZ8ozzzxj9o60jWODDsXIAJwsv8PloIaDcvwc3c1b9znUukxdQLZKlSpmBq7WbJYvX97v1/njjz/MppC6v5TWeQbS7bffbjakDvakluwCle7CriFFD2TE0BUAIOzs3r3bzLbduXOn2fxSZ+fqXoe6vEiLFi3kxIkTEi7LoOi6blFRUWKzlJQUiVQEHQBAWPa6a4hYsmSJtG7d2iw3ohtAf/fdd3LgwAEZPXp0+rkaMubNm+f1fO310OEYpb057v2l9FydBKMuX74sTz/9tDlXe4iee+45M+PXvW6bunTpkjmnYsWKUrRoUWnZsqUkJCR49bToa7qHfPQ99fUWL14sN910k5QsWVLuvvtu0+vj5s/75oVuhq1bIZUoUcLs5Tho0CCzlIsnbadeU90rsmfPnnL8+HGvx1944QXT+/XBBx+Y66dfu9q7d6/ZYV2/rtKlS8v999/vtaiv+3mzZ882PU3a0/XAAw/ImTNn0s/54osvTPuKFStmvn7dv1KXjQkWxw5dAUBeNv4M5OafLCroTXtrNCi88sor5sPQk/aePPzwwzJ37lyJj4/3qydl/fr1Zh02DUk333yzCVDqtddekzlz5siMGTNMKHnrrbdMYPJc8mTEiBHy5ZdfyqxZs6RGjRoyceJE6dSpk+lh0vXeMnP+/HmzM7l+2OtG1I888ogMHz7cvJe/75sX+p5TpkwxAWX37t0m6OjXoddLrVu3Th577DGZMGGCCVc6NKhDhFfSr1G/9q+++srsUqAL9rpDzg8//GACmwZSXQZGA5/b77//br6eBQsWmKE2DUOvvvqq+X5q4HvwwQfNddSApQFId20P5qalBB0AQFjR4Sr94NMQkBm9Xz9AdY017Wnx5eqrrzZ/au+B5/ZBb7/9towaNcp84Kp33nlH/vWvf6U/rr0M7777run90N4k9f7778vSpUvlww8/lGeffTbT9/vrr7/MMNt1111nbuvyKC+++KLf75ud6tWrZxqsPHnW6tSsWVNefvllefzxx9ODjgYr7WXS8KNuuOEGWb16tQk8Vw5XffTRR+nXT7/uX375xSzsqz1FSh/X8Ki9XE2bNjX3aSDSa1aqVClzu0+fPmbI0R10NCDde++9Jjiq7DbiDgTHBh1mXQFhVPCatDLYTUEECub/8pOTk82Qi/b0uGmvhW47pB/U7p4JDS133HFH+jmFChUyz/ntt9+yfG0dDnKHHKVF1EePHvX7fbOjvR/uAOHmHopz054r7a3Ztm2bWf1fg8XFixdNINK2advdIctN656uDDoaRNwhR+nzNOC4Q46qW7euGYLTx9xBR8OVZxs9v/4GDRrIXXfdZcKN9ox17NhRevXqJWXLlpVgcWyNjna3JSYmeo21AgBC7/rrrzdDUlmFCb1fPxjdH8J67pWhSANKqGgY8pRZ+3JLh6P0+ngeBQsW9Jph1rVrV6lfv74Zdtq4caP5j31uCoq1xidQX787xGmo056hhQsXmpCkvVt16tQxvUTB4tigAwAITzrE1KFDBzPUcuHCBa/HdCsgrW/RuhB3fY4GHs9iXx368hzOcdfkeG7irEWyupWQ53929fFNmzal39ZeGX3uqlWrvAKUPkc/pHPDn/fNCw02GiomTZokzZs3N8NSBw8ezDD0p3U6ntauXevztfV5OuVfDzftMNBC7JxcD/2+aS/ZuHHjzJR/vcZff/21BItjh64AAOFL61Z0jRod3tAaE+3J2Lp1q6mLqVatmqn3cGvXrp05X4dfNDToLCbPXgWt49GiZh2a0RoXnUGkgeOpp54yQzzaK3LjjTea3gWt/XEHKO3ReOKJJ8x7ujea1iJaDVFazJtbvt43L/Q1NYzpa3br1s2ENK0X8qQzvjRoaMG0Fhdr4feVw1aZ0dlROuSkxeC6wbYOiWmhs86K06UA/KEBS+t1dMhKvy96W2utsqrHCgR6dAAAYad27dpm0+datWqZWTvauzJw4EAzM2nNmjVeM56090LrRlq1aiUPPfSQmeGktShuOrSjs5Dee+89qVq1qvlwVxqIdAZQ3759TUjS2UQarNxTqZXOFrrvvvtMQe2tt95qZiJpMMhLTYk/75tbWgOj08t1Zle9evVM75eGKk/a06NF1VqUrOfrFP7nn3/e52trEPvmm2/M137nnXea4KPfH50B5y+dkr5ixQrp0qWL6W3S99Xvn7vYOxiiXMGs9ooAWqilyV4LxPQbACByi5H3nfIe5gi2itUn5tt7De1wQ66ep0WoWv/guRYKMqdDPtqzoMHqpZdesv59I0F2P7/+fn4zdAUAQdzlnB3Ow9eePXvSFyTUhQF1+Es/VLVXyMb3dSrHBh2mlwPIj4UFA7WoIAJPF9bT9V50qEsHN3SoR6dmB7NeJJTv61SODTps6gkAzqZ1PZ4zqmx/X6eiGBkAAFiLoAMAAKzl2KErAPmArR3ChsMn2CJC+bMthi8EHQCwmC6cp+uf6KJsuoJwIBalA/IjmOuWFfpzq8Xb7tWtc4OgAwAW072FdDXg/fv3m32QgEiiCz/qitQadnKLoAMguIsBIuR05V1daTiUG10CuQnpuqp1XnshHRt0WEcHgJMWFdQPDT0Ap3Fs0GEdHQCR5M2lO4K2TQRgM6aXAwAAaxF0AACAtQg6AADAWgQdAABgLccWIwPIGlPHAdiCoAMgd9jeAUAEYOgKAABYi6ADAACs5dihK1ZGBpAf5kfv8nlO9wCtnuzPooKKhQXhJNFOXhk5MTFREhISQt0UAAAQJI4NOgAAwH4EHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAazl2wUAAcCp/FhZkUUHYgqADICM27Ay71ZMr5ktLAPsQdACHid8cH+omAEC+iXbyXld169aVpk2bhropAAAgSBwbdNjrCgAA+zk26AAAAPsRdAAAgLUIOgAAwFoEHQAAYC2mlwNABGi+d5rPc9bGDMyXtgCRhKADAA4KQ4pABCdh6AoAAFiLHh0AsGSbiO5p1+dLW4BIQo8OAACwFkEHAABYi6ErAEAGby7d4fOqDO1wA1cOYY8eHQAAYC2CDgAAsBZDV4BF4jfHh7oJABBW6NEBAADWIugAAABrOXboKi4uzhypqamhbgqQv5JWcsUBOIZjg87gwYPNcfr0aSlTpkyomwMEBiEGALw4NugAgBO3iVAVg94SIHxQowMAAKxFjw4QAZg2DgC5Q48OAACwFkEHAABYi6ErAHCY5nun+TxnbczAfGkLEGz06AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIvp5QCAXHlz6Q6f5wztcANXFyFFjw4AALAWPTpAiLGPFQAED0EHABxmfvQun+dUzJeWAMFH0AEABA11PAg1gg4QCZJWhroFABCRKEYGAADWIugAAABrMXQFBBEzqgAgtOjRAQAA1qJHBwCQQfO903xelbUxA7lyCHsR36Ozb98+adOmjdStW1fq168vn3/+eaibBAAAwkTE9+gULFhQJk+eLA0bNpTDhw9L48aNpUuXLlKiRIlQNw0AAIRYxAedKlWqmENVrlxZKlSoICdOnCDoAACA3A1dXb58Wb777jt577335MyZM+a+gwcPytmzZ3P8WitWrJBu3bpJ1apVJSoqSubNm5fhnLi4OKlZs6YULVpUbrvtNlm/fn2mr7Vx40ZJTU2Va665JhdfFQAAEKcHnT179sgtt9wi99xzjwwePFj+/PNPc/9rr70mw4cPz3EDzp07Jw0aNDBhJjNz586VYcOGydixY2XTpk3m3E6dOsnRo0e9ztNenL59+8q0adkX0F26dElOnz7tdQAAADvleOhqyJAh0qRJE9myZYuUL18+/f6ePXvKgAEDctyAzp07myMrb7zxhnnd/v37m9tTp06Vb7/9VqZPny4jR45MDy89evQwt2+//fZs32/ChAkybty4HLcTAOCNmVmwskdn5cqV8vzzz0vhwoW97tehpQMHDgSybZKSkmKGo9q3b59+X3R0tLm9Zs0ac9vlckm/fv2kXbt20qdPH5+vOWrUKElOTk4/dNYWAACwU46DTlpamqmDudL+/fulVKlSEkjHjh0z71WpUiWv+/W2zrBSq1atMsNbWtujM6/0+OWXX7J8zSJFikjp0qW9DgAAYKccD1117NjRTOd218JoAbEWIWsNjU7rzm8tW7Y04QuIWOxMDgDhE3QmTZpkioF1gb6LFy/KQw89JDt37jTTuj/55JOANk5fs0CBAnLkyBGv+/W2TiUHAAAI6NBV9erVTSHy6NGjZejQodKoUSN59dVX5aeffpKKFStKIGkdkC4AuGzZsvT7tPdGb7do0SJPr62zvDSsNW3aNAAtBQAAVvTo6Lo3OrPp4YcfNofn2jr62J133pmj19Nhr127dqXfTkpKks2bN0u5cuUkJibGTC2PjY01M72aNWtmhs10Srp7FlZu6dR4PXR6eZkyZfL0WgAAwJKg07ZtWzl06FCG3hudwaSPZVaonJ0NGzaY57lpsFEabmbOnCm9e/c2a/WMGTPGFCBrsfGiRYsyFCgDAADkOejodG4tQL7S8ePHc7Xtgm7Iqa+ZnSeffNIcAAAAQQk69957r/lTQ46uW6PTtN20F+fnn3/2uVgfAABXenPpDr8uytAON3DxELyg465j0d4XXS+nWLFiXkXDzZs3z9XKyKGixch65HSoDQAAWBh0ZsyYkb4Csu5plZthqnBCMTIAAPbLcY2OLgwIAEA4DnExvIU8Bx31xRdfyGeffSZ79+41+1F50h3GAQCRbX70f5f9yEr3tOvzpS1Avi4YOGXKFLOGjU7v1kUCdW0b3cV89+7d2e5CDgAAEPY9OvHx8WafqwcffNCsczNixAipVauWWefmxIkTwWklEIbiN8eHuglA2Gu+9z/7ImZnbczAfGkLnCnHPTo6XOWeRq4zr86cOWP+3qdPn4DvdRVMbAEBAID9chx0dDNNd8+NbtGwdu3a9K0bfC38F26zrhITEyUhISHUTQEAAOESdNq1ayfz5883f9daHd3Ys0OHDmarhp49ewajjQAAAPlTo6P1ObqDuLtXRAuRV69eLd27d5e//e1vuWsFEGaovwEAhwad6Ohoc7g98MAD5lAHDhyQatWqBbaFQCRLWhnqFgCAo+V46Cozuqv4U089JbVr1w7EywEAAORvj87Jkydl0KBBsnTpUrO31ciRI82O4i+88IK8/vrrUr9+/fRtIgBHoLcGyLcp6Ipp6Ahq0NFgo7U4unP54sWLTRHyokWLzDDWv//9b7OpZyRhU08AAOznd9BZuHChWSBQZ11pT44uEtiwYUMZP368RCI29QSAvGGbCFgVdA4ePCg33XRT+g7mRYsWlUceeSSYbQMCjtlUAOAsfhcj62KABQv+NxcVKFDArIwMAAAQ8T06GnTuuuuu9LBz4cIF6datmylM9sTu5QAAIOKCztixY71u33PPPcFoDwAAQOiDDgAAgCMWDAQAALBiCwhbsI4OANjnzaU7fJ4ztMMN+dIWhAfH9ujoOjqJiYmSkJAQ6qYAAIAgcWyPDgAg+FhUEBEXdD766CPp3bu3FClSxOv+lJQU+fTTT6Vv376BbB8AwHL+hCFVMegtgY1yPHTVv39/SU5OznD/mTNnzGMAAAARG3R04cCoqKgM9+/fv1/KlCkTqHYBAADk39BVo0aNTMDRw3OFZJWamipJSUly9913571FAAAA+R10evToYf7cvHmzdOrUSUqWLJn+mG4DoRt93nfffYFqFwAAQP6vjKyBRouRdfdyAAAAq2ZdxcbGps+yOnr0qKSlpXk9HhMTE7jWAQAQYCwq6Cw5Djo7d+6URx99VFavXp1pkbLW60QCVkYGAMB+OQ46/fr1M4XICxYskCpVqmQ6AytSVkbW4/Tp08wWAwDAUjkOOlqMvHHjRrnxxhuD0yIAADLRfO80n9dlbcxArh3yto5O3bp15dixYzl9GgAAQPgHnddee01GjBgh33//vRw/ftwM/XgeAAAAETt01b59e/OnLhoYycXIAAD7MLyFPAed5cuX5/QpAAAAkRF0WrduHZyWANmI3xzv8/oMajiIawgAyFuNjlq5cqU88sgjcvvtt8uBAwfMfbNnz5Yff/wxNy8HAAAQHkHnyy+/NHtdFStWTDZt2iSXLl0y9ycnJ8v48eOD0UYAAID8CTovv/yyTJ06Vd5//30pVKhQ+v133HGHCT4AAAARW6Ozfft2ufPOOzPcX6ZMGTl16lSg2gUEpY5Hklb692LXtuI7ADgY+2E5OOhUrlxZdu3aZXYx96T1ObVq1ZJIwV5XCEggAgDYNXQ1YMAAGTJkiKxbt86sm3Pw4EGZM2eODB8+XJ544gmJFLrPVWJioiQkJIS6KQAAIFx6dEaOHClpaWlmwcDz58+bYawiRYqYoPPUU08Fp5UAAAD5EXS0F2f06NHy7LPPmiGss2fPmv2vSpYsmZv3BwAACJ+hq3/+85+mJ6dw4cIm4DRr1oyQAwAA7Ag6Q4cOlYoVK8pDDz0k//rXv9jbCgAA2BN0Dh06JJ9++qkZwrr//vulSpUqprB39erVwWkhAABAfgWdggULSteuXc1Mq6NHj8qbb74pf/zxh7Rt21auu+663LYDAAAg9MXInooXL262gzh58qTs2bNHfvvtt8C1DAAAIBSbemoxsvbodOnSRapVqyaTJ0+Wnj17ytatW/PaHgAAgND16DzwwAOyYMEC05ujNTr/93//Jy1atAhciwAAAEIVdAoUKCCfffaZGbLSvwMAAFgTdHTICgAAp/Nn4081tMMNQW8LAlCjo/U4ycnJ6bdfffVVr93Kjx8/bhYQBAAAiLgencWLF8ulS5fSb48fP97U6Fx11VXm9uXLl2X79u3BaSUAAAHSfO80n+esjRnI9XZaj47L5cr2NgAAgBXTywEAAKwKOrrlgx5X3gcAABDxNTo6VNWvXz8pUqSIuX3x4kV5/PHHpUSJEua2Z/1OJIiLizNHampqqJuCQElaybUELDY/eldAXqd72vUBeR1YFnRiY2O9bj/yyCMZzunbt69ECt2IVI/Tp09LmTJlQt0cAAAQyqAzY8aMYLw/4B96awBE4Oyt/3g9yC1BdihGBgAA1iLoAAAAa+V4CwgAABDYrSLYJiJ46NEBAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLdXSQK/Gb4wN25QY1HMR3AUBYbQ5aMV9agvxAjw4AALAWQQcAAFiLoSsAAHK9MznCHT06AADAWgQdAABgLYauEHLxXz8Y6iYAQNhjF/TcoUcHAABYi6ADAACsxdAVAACWYHgrI3p0AACAtQg6AADAWlYEnZ49e0rZsmWlV69eoW4KAAAII1bU6AwZMkQeffRRmTVrVqibAgBAUGpr4OAenTZt2kipUqVC3QwAABBmQh50VqxYId26dZOqVatKVFSUzJs3L8M5cXFxUrNmTSlatKjcdtttsn79+pC0FQAARJaQB51z585JgwYNTJjJzNy5c2XYsGEyduxY2bRpkzm3U6dOcvTo0Vy936VLl+T06dNeBwAAsFPIa3Q6d+5sjqy88cYbMmDAAOnfv7+5PXXqVPn2229l+vTpMnLkyBy/34QJE2TcuHF5ajMAwG7zo3f5PKd72vX50hZEeI9OdlJSUmTjxo3Svn379Puio6PN7TVr1uTqNUeNGiXJycnpx759+wLYYgAAEE5C3qOTnWPHjklqaqpUqlTJ6369vW3btvTbGny2bNlihsGqV68un3/+ubRo0SLT1yxSpIg5AACA/cI66Pjru+++C3UTAABAGArroasKFSpIgQIF5MiRI1736+3KlSuHrF0AACAyhHWPTuHChaVx48aybNky6dGjh7kvLS3N3H7yySfz9No6y0sPHRqDt/ivH/R9Sa5txWUDAD803zvN5zlrYwbm23uJvC5OEvKgc/bsWdm167/V7UlJSbJ582YpV66cxMTEmKnlsbGx0qRJE2nWrJlMnjzZ1OK4Z2Hl1uDBg82h08vLlCkTgK8EAACEm5AHnQ0bNkjbtm3Tb2uwURpuZs6cKb1795Y///xTxowZI4cPH5aGDRvKokWLMhQoAwAAhF3Q0e0bXC5XtufoMFVeh6oAAIDzhHUxMgAAQET36ISKU4uR4zfH5+8bJq3M3/cDgDBaPdnfFZT9KyIOvx3Vh3a4QcKdY3t0tBA5MTFREhISQt0UAAAQJI4NOgAAwH4EHQAAYC2CDgAAsBZBBwAAWItZVw6bdQUAgJO2k3Bsjw6zrgAAsJ9jgw4AALAfQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsxvZzp5bnDZp0AELDNP/3Z+DM/N+tUzcUOju3RYXo5AAD2c2zQAQAA9iPoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwlmODTlxcnNStW1eaNm0a6qYAAIAgcWzQYR0dAADs59igAwAA7EfQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYq6A4eGVkPVJTUyUiLJ/g+5y2o/KjJQAARAzH9uiwMjIAAPZzbNABAAD2I+gAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGux11Wk7HXlh/jN8b5PSlqZH00BAOTA/OhdAble3dOu93lO873TxEkc26PDXlcAANjPsUEHAADYj6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsFZBcai4uDhzpKamhrYhyyeE9v0BALCYY3t0Bg8eLImJiZKQkBDqpgAAgCBxbNABAAD2I+gAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYy4qgs2DBAqlTp47Url1bPvjgg1A3BwAAhImCEuEuX74sw4YNk+XLl0uZMmWkcePG0rNnTylfvnyomwYAAEIs4nt01q9fLzfffLNUq1ZNSpYsKZ07d5YlS5aEulkAACAMhDzorFixQrp16yZVq1aVqKgomTdvXoZz4uLipGbNmlK0aFG57bbbTLhxO3jwoAk5bvr3AwcO5Fv7AQBA+Ap50Dl37pw0aNDAhJnMzJ071wxNjR07VjZt2mTO7dSpkxw9ejRX73fp0iU5ffq01wEAAOwU8hodHWrSIytvvPGGDBgwQPr3729uT506Vb799luZPn26jBw50vQEefbg6N+bNWuW5etNmDBBxo0bJ/li+YT8eR8AAPw0P3pXwK7V2qU7fJ4ztMMN4ugeneykpKTIxo0bpX379un3RUdHm9tr1qwxtzXU/PrrrybgnD17VhYuXGh6fLIyatQoSU5OTj/27duXL18LAABwYI9Odo4dOyapqalSqVIlr/v19rZt28zfCxYsKJMmTZK2bdtKWlqajBgxItsZV0WKFDEHAACwX1gHHX91797dHAAAABEzdFWhQgUpUKCAHDlyxOt+vV25cuWQtQsAAESGsA46hQsXNgsALlu2LP0+HZ7S2y1atMjTa+ssr7p160rTpk0D0FIAABCOQj50pQXEu3b9twI8KSlJNm/eLOXKlZOYmBgztTw2NlaaNGliCo8nT55spqS7Z2Hl1uDBg82h08t1RWUAAGCfkAedDRs2mEJiNw02SsPNzJkzpXfv3vLnn3/KmDFj5PDhw9KwYUNZtGhRhgJlAACAsAs6bdq0EZfLle05Tz75pDkAAACsqdEBAADIC8cGHYqRAQCwn2ODjhYiJyYmSkJCQqibAgAAgsSxQQcAANiPoAMAAKxF0AEAANYi6AAAAGuFfB2dUM660uPy5cvmtq6QHHDnLkp+ulCosO+Tzv+VH00BQuLSBX6+4Wzn0i75POdSdOB+Ty6eO+vznKB8vnq8rq+1+KJcvs6w3P79++Waa64JdTMAAEAu7Nu3T6pXr57l444POrpJ6MGDB6VUqVISFRWVfmF0s88rp577c5/7tiZNDVD6DShdurQEWmZtCcRzfJ2T1eM5uTa2XCtf53GtQnutVDB/tmz6uQrHa+Xv84J1ra68j2sVfp+F2k9z5swZqVq1qkRHZ12J49ihKze9OJklwQIFCmT4pvhz35W39e/B+PDOrC2BeI6vc7J6PDfXJtKvla/zuFbhca2C9bNl489VOF0rf58XrGt15X1cq/D8LPRnU26KkbNZUDA392V2TjDk5n38eY6vc7J6PDfXJtKvla/zuFZcK36ucv876O/zgvU7eOV94fzvVbhdq/y8Xv5w/NBVsGh3nSbN5OTkoKRYm3CtuFb8bIUev4dcK1t/rujRCZIiRYrI2LFjzZ/gWvFzFRr8HnKt+LkKrXD4HaRHBwAAWIseHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLohMCCBQukTp06Urt2bfnggw9C0YSI0rNnTylbtqz06tUr1E0Ja7rEeps2baRu3bpSv359+fzzz0PdpLB16tQpadKkiTRs2FDq1asn77//fqibFPbOnz8vNWrUkOHDh4e6KWGtZs2a5vdPf7batm0b6uaEtaSkJHON9N+sW265Rc6dOxeU92F6eT7T3dL1m7p8+XKziFLjxo1l9erVUr58+fxuSsT4/vvvzX4ms2bNki+++CLUzQlbhw4dkiNHjph/YA8fPmx+tnbs2CElSpQIddPCTmpqqly6dEmKFy9u/nHVsLNhwwZ+D7MxevRo2bVrl9m36PXXX8+/b1YEBp1ff/1VSpYsGeqmhL3WrVvLyy+/LK1atZITJ06YBQULFgz8zlT06OSz9evXy8033yzVqlUzvwidO3eWJUuW5HczIor2Uuimq8helSpVTMhRlStXlgoVKph/PJCR7sOjIUdp4NHNAfVA5nbu3Cnbtm0z/14BgbB161YpVKiQCTmqXLlyQQk5iqCTQytWrJBu3bqZ3VJ1t/N58+ZlOCcuLs6k+qJFi8ptt91mwo2b7pSuIcdN/37gwAGxVV6vl5ME8lpt3LjR9Fro/75tFIhrpcNXDRo0MJv6PvvssyYY2igQ10qHqyZMmCC2C8S10udpT4Xu3j1nzhyx1Yo8XisNz/qffX2NW2+9VcaPHx+0thJ0cki7ufUfR/0GZmbu3LkybNgws+T1pk2bzLmdOnWSo0ePihNxvfL/WmkvTt++fWXatGliq0Bcq6uuukq2bNli6gQ+/vhjM+xno7xeq2+++UZuuOEGc9guED9XP/74o/mPxvz5882H988//yw2OpfHa6VlHCtXrpT4+HhZs2aNLF261BxB4UKu6eX7+uuvve5r1qyZa/Dgwem3U1NTXVWrVnVNmDDB3F61apWrR48e6Y8PGTLENWfOHEd8F3JzvdyWL1/uuu+++1xOkdtrdfHiRVerVq1cH330kcsp8vJz5fbEE0+4Pv/8c5ftcnOtRo4c6apevbqrRo0arvLly7tKly7tGjdunMt2gfi5Gj58uGvGjBku20kurtXq1atdHTt2TH984sSJ5ggGenQCKCUlxST59u3bp98XHR1tbmtiVc2aNTOFajpcdfbsWVm4cKFJuU7kz/WC/9dK/73p16+ftGvXTvr06ePYS+fPtdLeGy1wV7qrsnbD60xIp/HnWumQlc7o++OPP0wR8oABA2TMmDHiNP5cK+3lcP9c6b/v//73v01NptOk+HGtdGhPe3dOnjwpaWlp5nfwpptuCkp7glP541DHjh0zdRGVKlXyul9vayGf0mKrSZMmmSl1+s0dMWKEY2d6+HO9lP5y6BCD/iOi9RQ6bbpFixbiJP5cq1WrVpnuYp3a6h4vnz17tpm26ST+XKs9e/bIwIED04uQn3rqKcddp5z8DsK/a6UBWpfDUHquhkL9QHeaY35+FurQ3p133ml+Bzt27Chdu3YNSnsIOiHQvXt3c8A/3333HZfKDy1btjThGb5pz+rmzZu5VDmkPYbIWq1atcx/yuAfncWXHzP5GLoKIJ21odNWryxq1Ns63RdcL362go/fQ64VP1ehVSHMPgsJOgFUuHBhs0jbsmXL0u/T/2HrbacNtfiD68W14ueK38FIwb9XkXutGLrKIS0w09VB3XRqqnaB62JHMTExZjpdbGysWV5eu8cnT55sakv69+8vTsT14lrxc8XvYKTg3ytLr1VQ5nJZTKc562W78oiNjU0/5+2333bFxMS4ChcubKbYrV271uVUXC+uFT9X/A5GCv69svNasdcVAACwFjU6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAHKtNmzYSFRVljlDuZj5z5sz0djzzzDMhawdgI4IOAL/069cv/cPY87j77rsj+goOGDBADh06JPXq1fO6/8svv5R27dpJ2bJlpVixYlKnTh159NFH5aeffvLrdVNSUswuzq+++mqmj7/00ktSqVIl+euvv6R3796mDWz+CwQeQQeA3zTU6Aey5/HJJ58E9QpqYAim4sWLS+XKlaVgwf/ucfzcc8+Z8NGwYUOZP3++bN++XT7++GOpVauWjBo1yu8dnB955BGZMWNGhsdcLpfpxenbt68UKlTIBCltgz4HQGARdAD4rUiRIuYD2fPQHg837eH54IMPpGfPniZA1K5d2wQFT7/++qt07txZSpYsaXo0+vTpI8eOHfMaTnryySfNEI72iHTq1Mncr6+jr1e0aFFp27atzJo1y7zfqVOnzK7IpUuXli+++MLrvebNmyclSpSQM2fO+P01rl27ViZOnChvvPGGOVq1amV2Y27cuLE8//zzsnDhQq/zv/nmG7n11ltNuzQIjRs3Ti5fvmwee+yxx2THjh3y448/ej3nhx9+kN27d5vHAQQXQQdAQOkH/f333y8///yzdOnSRR5++GE5ceKEeUxDiQ4HNWrUSDZs2CCLFi2SI0eOmPM9aYjR3o1Vq1bJ1KlTJSkpSXr16iU9evSQLVu2yN/+9jcZPXp0+vkaZh544IEMvSd6W59XqlQpv9uvPVQawgYNGpTp4xqu3FauXGl6ZYYMGSKJiYny3nvvmZ6aV155xTx+yy23SNOmTWX69OkZ2nX77bfLjTfe6He7AORSSPZMBxBxYmNjXQUKFHCVKFHC63jllVfSz9F/Up5//vn022fPnjX3LVy40Nx+6aWXXB07dvR63X379plztm/fbm63bt3a1ahRI69znnvuOVe9evW87hs9erR53smTJ83tdevWmfYdPHjQ3D5y5IirYMGCru+//z7Lr0nfa8iQIV733X333a769et73Tdp0iSvr/nUqVPm/rvuuss1fvx4r3Nnz57tqlKlSvrtqVOnukqWLOk6c+aMuX369GlX8eLFXR988IFf7QGQN/ToAPCbDhnp7CTP4/HHH/c6p379+l49LTqkdPToUXNbe2OWL19uekzch7tX4/fff09/ng4TedIaGe0Z8dSsWbMMt2+++WbTG6T++c9/So0aNeTOO+/M83dYi5D1a9UeGx0m+0+m+8/X8+KLL3p9Pe7i5vPnz5tzHnzwQUlNTZXPPvvM3J47d65ER0ebGiAAwfff6jsA8EGDy/XXX5/tOVpce+VQT1pamvn72bNnpVu3bvLaa69leF6VKlW83ic3/vd//1fi4uJk5MiRZniof//+XkNN/tA6IK2p0dlQ7q/lqquuMsf+/fu9ztWvR4fq7r333gyvozU7SoOeDp9pezQw6Z86VKehCEDw0aMDIN9o0e7WrVulZs2aJjB5HtmFG53arTU9nhISEjKcp7Oc9uzZI1OmTDE1M7GxsTluo/bAaICJj4/36+vR3qYrvxY9tNfGTYuONTwtWLBAVq9eTREykI8IOgD8dunSJTl8+LDX4TljypfBgwebwmQNExpUdLhq8eLFpudFh3eyosXH27ZtM9O+dRaTDgNp0a/y7LHRGWDau/Lss89Kx44dpXr16jn+7upaNn//+9/NMWzYMBNQNDzpbKwPP/zQvJ87xIwZM0Y++ugj06ujAe63336TTz/91MzO8qTDZxp+tHBZh+q0EBlA/iDoAPCbzpLSISbPo2XLln4/v2rVqmYmlYYaDSI6K0mnkeuwkGcPyJWuvfZaM3X8q6++MjVA7777bvqsK53y7kl7T3TtHR0myq3XX3/drJujiwN27drVDGf9z//8jxmCW7NmjRmOUjr1XXtplixZYmqImjdvLm+++aapDfKk4Ujbc/LkyTy1C0DORWlFci6eBwAhpVO4der5vn37vO6fPXu2DB06VA4ePOhzAT5ds0cXBZw8ebKEg3BrD2ADenQARAStmdHhLl1oT8PMP/7xD68aHJ3lpENhuuWCDnX5u8qwvq4WBv/yyy8SKnPmzDFt0HV5AAQWPToAIoL20ujUbK3x0ZWKdUVl3Y7BvXXDCy+8YHp5tB5GVyv2Z1bTgQMH5MKFC+bv+pqh2oJBV27WhROVDuPpitAAAoOgAwAArMXQFQAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAAAQW/0/46sojNq1CGoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# kinematics\n", + "\n", + "nu_momenta = np.squeeze(data[\"primary_momentum\"])\n", + "\n", + "# muon\n", + "mu_flag = data[\"secondary_types\"]=='ParticleType.MuMinus'\n", + "mu_momenta = np.squeeze(data[\"secondary_momenta\"][mu_flag])\n", + "\n", + "# nu out\n", + "target_flag = data[\"secondary_types\"]=='ParticleType.Hadrons'\n", + "target_momenta = np.squeeze(data[\"secondary_momenta\"][target_flag])\n", + "\n", + "kwargs = {\"bins\":np.logspace(0,6,50),\n", + " #\"weights\":data[\"event_weight\"],\n", + " \"alpha\":0.5}\n", + "\n", + "# Energy\n", + "plt.hist(nu_momenta[:,0],**kwargs,label=r\"Initial $\\nu$\")\n", + "plt.hist(mu_momenta[:,0],**kwargs,label=r\"Outgoing $\\mu$\")\n", + "plt.hist(target_momenta[:,0],**kwargs,label=r\"Outgoing Hadrons\")\n", + "plt.legend()\n", + "plt.loglog()\n", + "plt.xlabel(\"Energy [GeV]\")\n", + "plt.ylabel(\"Event Rate\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eeb54511-c655-478e-9fac-43246279faf4", + "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 +}