Warning, /snippets/Tracking/check_dd4hep_geo.ipynb is written in an unsupported language. File is not indexed.
0001 {
0002 "cells": [
0003 {
0004 "cell_type": "markdown",
0005 "id": "53ea04e0-3cc3-4fff-8816-40f6ee85252b",
0006 "metadata": {},
0007 "source": [
0008 "# def\n",
0009 "\n",
0010 "* This is an example code to access sensitive detector modules and generate random position on module --> random cell ID. It can be used to inspect the geometry component, and create random sensor noise for a given epic geometry. \n",
0011 "\n",
0012 "* Please start the jupyter server within eic-shell.\n",
0013 "\n",
0014 "Shujie Li, 082025\n"
0015 ]
0016 },
0017 {
0018 "cell_type": "code",
0019 "execution_count": null,
0020 "id": "2528a4cb-7d28-4df2-b76a-70cbe9ca5e36",
0021 "metadata": {},
0022 "outputs": [],
0023 "source": [
0024 "import dd4hep\n",
0025 "import numpy as np\n",
0026 "from dd4hep import DetElement, VolumeManager\n",
0027 "import ROOT\n",
0028 "import DDRec\n",
0029 "import pandas as pd\n",
0030 "import matplotlib.pyplot as plt\n",
0031 "import random\n",
0032 "import math\n",
0033 "import os\n",
0034 "\n",
0035 "# Get DETECTOR_PATH\n",
0036 "detector_path = os.environ.get('DETECTOR_PATH')\n",
0037 "xml_path = detector_path+\"/epic_craterlake_tracking_only.xml\"\n",
0038 "detector = dd4hep.Detector.getInstance()\n",
0039 "detector.fromCompact(xml_path)\n",
0040 "converter = DDRec.CellIDPositionConverter(detector)\n"
0041 ]
0042 },
0043 {
0044 "cell_type": "code",
0045 "execution_count": 3,
0046 "id": "3142d1e6-0410-4338-bd21-f6c55708b125",
0047 "metadata": {},
0048 "outputs": [],
0049 "source": [
0050 "\n",
0051 "def random_point_in_volume(shape):\n",
0052 " \"\"\"\n",
0053 " Generate a random point inside a TGeo shape volume.\n",
0054 " \n",
0055 " Args:\n",
0056 " shape: TGeo shape object with .type() and .dimensions() methods\n",
0057 " \n",
0058 " Returns:\n",
0059 " tuple: (x, y, z) coordinates of random point inside the volume\n",
0060 " \"\"\"\n",
0061 " \n",
0062 " shape_type = shape.type()\n",
0063 " dims = shape.dimensions()\n",
0064 " \n",
0065 " if shape_type == \"TGeoTrd2\":\n",
0066 " # TRD2: [dx1, dx2, dy1, dy2, dz]\n",
0067 " dx1, dx2, dy1, dy2, dz = dims[0], dims[1], dims[2], dims[3], dims[4]\n",
0068 " \n",
0069 " # Random z coordinate\n",
0070 " z = random.uniform(-dz, dz)\n",
0071 " \n",
0072 " # Linear interpolation factor\n",
0073 " t = (z + dz) / (2 * dz)\n",
0074 " \n",
0075 " # Interpolate half-widths at this z\n",
0076 " dx_z = dx1 + t * (dx2 - dx1)\n",
0077 " dy_z = dy1 + t * (dy2 - dy1)\n",
0078 " \n",
0079 " # Random x, y within interpolated bounds\n",
0080 " x = random.uniform(-dx_z, dx_z)\n",
0081 " y = random.uniform(-dy_z, dy_z)\n",
0082 " \n",
0083 " return x, y, z\n",
0084 " \n",
0085 " elif shape_type == \"TGeoBBox\" or shape_type == \"Box\":\n",
0086 " # Box: [dx, dy, dz] (half-widths)\n",
0087 " dx, dy, dz = dims[0], dims[1], dims[2]\n",
0088 " \n",
0089 " x = random.uniform(-dx, dx)\n",
0090 " y = random.uniform(-dy, dy)\n",
0091 " z = random.uniform(-dz, dz)\n",
0092 " \n",
0093 " return x, y, z\n",
0094 " \n",
0095 " elif shape_type == \"TGeoTube\" or shape_type == \"Tube\":\n",
0096 " # Tube: [rmin, rmax, dz]\n",
0097 " rmin, rmax, dz = dims[0], dims[1], dims[2]\n",
0098 " \n",
0099 " # Random cylindrical coordinates\n",
0100 " r = math.sqrt(random.uniform(rmin*rmin, rmax*rmax)) # Uniform in area\n",
0101 " phi = random.uniform(0, 2 * math.pi)\n",
0102 " z = random.uniform(-dz, dz)\n",
0103 " \n",
0104 " # Convert to Cartesian\n",
0105 " x = r * math.cos(phi)\n",
0106 " y = r * math.sin(phi)\n",
0107 " \n",
0108 " return x, y, z\n",
0109 " \n",
0110 " elif shape_type == \"TGeoTubs\" or shape_type == \"Tubs\":\n",
0111 " # Tube segment: [rmin, rmax, dz, phi1, phi2]\n",
0112 " rmin, rmax, dz, phi1, phi2 = dims[0], dims[1], dims[2], dims[3], dims[4]\n",
0113 " \n",
0114 " # Convert angles to radians if needed\n",
0115 " phi1_rad = math.radians(phi1) if phi1 > 2*math.pi else phi1\n",
0116 " phi2_rad = math.radians(phi2) if phi2 > 2*math.pi else phi2\n",
0117 " \n",
0118 " # Random cylindrical coordinates\n",
0119 " r = math.sqrt(random.uniform(rmin*rmin, rmax*rmax))\n",
0120 " phi = random.uniform(phi1_rad, phi2_rad)\n",
0121 " z = random.uniform(-dz, dz)\n",
0122 " \n",
0123 " # Convert to Cartesian\n",
0124 " x = r * math.cos(phi)\n",
0125 " y = r * math.sin(phi)\n",
0126 " \n",
0127 " return x, y, z\n",
0128 " \n",
0129 " elif shape_type == \"TGeoSphere\" or shape_type == \"Sphere\":\n",
0130 " # Sphere: [rmin, rmax, theta1, theta2, phi1, phi2]\n",
0131 " if len(dims) >= 6:\n",
0132 " rmin, rmax, theta1, theta2, phi1, phi2 = dims[:6]\n",
0133 " else:\n",
0134 " rmin, rmax = dims[0], dims[1]\n",
0135 " theta1, theta2 = 0, 180 # Full sphere in theta\n",
0136 " phi1, phi2 = 0, 360 # Full sphere in phi\n",
0137 " \n",
0138 " # Convert angles to radians\n",
0139 " theta1_rad = math.radians(theta1) if theta1 > 2*math.pi else theta1\n",
0140 " theta2_rad = math.radians(theta2) if theta2 > 2*math.pi else theta2\n",
0141 " phi1_rad = math.radians(phi1) if phi1 > 2*math.pi else phi1\n",
0142 " phi2_rad = math.radians(phi2) if phi2 > 2*math.pi else phi2\n",
0143 " \n",
0144 " # Random spherical coordinates\n",
0145 " r = (rmin**3 + random.random() * (rmax**3 - rmin**3))**(1/3) # Uniform in volume\n",
0146 " cos_theta = math.cos(theta1_rad) + random.random() * (math.cos(theta2_rad) - math.cos(theta1_rad))\n",
0147 " theta = math.acos(cos_theta)\n",
0148 " phi = random.uniform(phi1_rad, phi2_rad)\n",
0149 " \n",
0150 " # Convert to Cartesian\n",
0151 " x = r * math.sin(theta) * math.cos(phi)\n",
0152 " y = r * math.sin(theta) * math.sin(phi)\n",
0153 " z = r * math.cos(theta)\n",
0154 " \n",
0155 " return x, y, z\n",
0156 " \n",
0157 " elif shape_type == \"TGeoCone\" or shape_type == \"Cone\":\n",
0158 " # Cone: [dz, rmin1, rmax1, rmin2, rmax2]\n",
0159 " dz, rmin1, rmax1, rmin2, rmax2 = dims[0], dims[1], dims[2], dims[3], dims[4]\n",
0160 " \n",
0161 " # Random z coordinate\n",
0162 " z = random.uniform(-dz, dz)\n",
0163 " \n",
0164 " # Linear interpolation factor\n",
0165 " t = (z + dz) / (2 * dz)\n",
0166 " \n",
0167 " # Interpolate radii at this z\n",
0168 " rmin_z = rmin1 + t * (rmin2 - rmin1)\n",
0169 " rmax_z = rmax1 + t * (rmax2 - rmax1)\n",
0170 " \n",
0171 " # Random cylindrical coordinates\n",
0172 " r = math.sqrt(random.uniform(rmin_z*rmin_z, rmax_z*rmax_z))\n",
0173 " phi = random.uniform(0, 2 * math.pi)\n",
0174 " \n",
0175 " # Convert to Cartesian\n",
0176 " x = r * math.cos(phi)\n",
0177 " y = r * math.sin(phi)\n",
0178 " \n",
0179 " return x, y, z\n",
0180 " \n",
0181 " elif shape_type == \"TGeoTrd1\" or shape_type == \"Trd1\":\n",
0182 " # TRD1: [dx1, dx2, dy, dz] (trapezoid in x only)\n",
0183 " dx1, dx2, dy, dz = dims[0], dims[1], dims[2], dims[3]\n",
0184 " \n",
0185 " # Random z coordinate\n",
0186 " z = random.uniform(-dz, dz)\n",
0187 " \n",
0188 " # Linear interpolation for x\n",
0189 " t = (z + dz) / (2 * dz)\n",
0190 " dx_z = dx1 + t * (dx2 - dx1)\n",
0191 " \n",
0192 " # Random coordinates\n",
0193 " x = random.uniform(-dx_z, dx_z)\n",
0194 " y = random.uniform(-dy, dy)\n",
0195 " \n",
0196 " return x, y, z\n",
0197 " \n",
0198 " elif shape_type == \"TGeoEltu\" or shape_type == \"Eltu\":\n",
0199 " # Elliptical tube: [a, b, dz] where a, b are semi-axes\n",
0200 " a, b, dz = dims[0], dims[1], dims[2]\n",
0201 " \n",
0202 " # Generate point in ellipse using rejection method\n",
0203 " while True:\n",
0204 " x = random.uniform(-a, a)\n",
0205 " y = random.uniform(-b, b)\n",
0206 " if (x/a)**2 + (y/b)**2 <= 1:\n",
0207 " break\n",
0208 " \n",
0209 " z = random.uniform(-dz, dz)\n",
0210 " return x, y, z\n",
0211 " \n",
0212 " elif shape_type == \"TGeoTrap\" or shape_type == \"Trap\":\n",
0213 " # General trapezoid - more complex, simplified version\n",
0214 " # Trap: [dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2]\n",
0215 " # Simplified: treat as box with average dimensions\n",
0216 " if len(dims) >= 11:\n",
0217 " dz = dims[0]\n",
0218 " h1, h2 = dims[3], dims[7]\n",
0219 " bl1, tl1 = dims[4], dims[5]\n",
0220 " bl2, tl2 = dims[8], dims[9]\n",
0221 " \n",
0222 " # Use average dimensions (simplified)\n",
0223 " dx = (bl1 + tl1 + bl2 + tl2) / 4\n",
0224 " dy = (h1 + h2) / 2\n",
0225 " \n",
0226 " x = random.uniform(-dx, dx)\n",
0227 " y = random.uniform(-dy, dy)\n",
0228 " z = random.uniform(-dz, dz)\n",
0229 " \n",
0230 " return x, y, z\n",
0231 " \n",
0232 " else:\n",
0233 " # Fallback for unknown shapes - try to treat as box\n",
0234 " print(f\"Warning: Unknown shape type '{shape_type}', treating as box\")\n",
0235 " if len(dims) >= 3:\n",
0236 " dx, dy, dz = dims[0], dims[1], dims[2]\n",
0237 " x = random.uniform(-dx, dx)\n",
0238 " y = random.uniform(-dy, dy)\n",
0239 " z = random.uniform(-dz, dz)\n",
0240 " return x, y, z\n",
0241 " else:\n",
0242 " raise ValueError(f\"Unsupported shape type: {shape_type} with dimensions: {dims}\")\n",
0243 "\n",
0244 "# Helper function for multiple points\n",
0245 "def random_points_in_volume(shape, n_points):\n",
0246 " \"\"\"Generate multiple random points in a volume\"\"\"\n",
0247 " points = []\n",
0248 " for _ in range(n_points):\n",
0249 " points.append(random_point_in_volume(shape))\n",
0250 " return points\n"
0251 ]
0252 },
0253 {
0254 "cell_type": "markdown",
0255 "id": "9d80062b-6dab-44a2-8eba-df00c3f0e3de",
0256 "metadata": {},
0257 "source": [
0258 "# analysis"
0259 ]
0260 },
0261 {
0262 "cell_type": "markdown",
0263 "id": "b2799827-ae2c-4202-982b-2c2eb58df98a",
0264 "metadata": {},
0265 "source": [
0266 "## get random position and cell ID for a given detector"
0267 ]
0268 },
0269 {
0270 "cell_type": "code",
0271 "execution_count": 3,
0272 "id": "fd32d3a6-7818-4713-ba4e-8d351811962f",
0273 "metadata": {},
0274 "outputs": [
0275 {
0276 "name": "stdout",
0277 "output_type": "stream",
0278 "text": [
0279 "B0ECal\n",
0280 "B0Tracker\n",
0281 "BackwardMPGD\n",
0282 "BackwardsTaggerStation\n",
0283 "BarrelTOF\n",
0284 "ForwardMPGD\n",
0285 "ForwardOffMTracker_station_1\n",
0286 "ForwardOffMTracker_station_2\n",
0287 "ForwardOffMTracker_station_3\n",
0288 "ForwardOffMTracker_station_4\n",
0289 "ForwardRomanPot_Station_1\n",
0290 "ForwardRomanPot_Station_2\n",
0291 "ForwardTOF\n",
0292 "HcalFarForwardZDC_SiPMonTile\n",
0293 "InnerMPGDBarrel\n",
0294 "InnerTrackerEndcapN\n",
0295 "InnerTrackerEndcapP\n",
0296 "LumiDirectPCAL\n",
0297 "LumiSpecCAL\n",
0298 "LumiSpecTracker\n",
0299 "MPGDOuterBarrel\n",
0300 "MiddleTrackerEndcapN\n",
0301 "MiddleTrackerEndcapP\n",
0302 "OuterSiBarrel\n",
0303 "OuterTrackerEndcapN\n",
0304 "OuterTrackerEndcapP\n",
0305 "SagittaSiBarrel\n",
0306 "VertexBarrel\n",
0307 "ZDC_Crystal\n"
0308 ]
0309 }
0310 ],
0311 "source": [
0312 "for i,j in detector.sensitiveDetectors():\n",
0313 "\tprint(i)"
0314 ]
0315 },
0316 {
0317 "cell_type": "code",
0318 "execution_count": null,
0319 "id": "f437947a-262e-4fa5-9cdb-dbda4247e9b6",
0320 "metadata": {},
0321 "outputs": [
0322 {
0323 "name": "stdout",
0324 "output_type": "stream",
0325 "text": [
0326 "VertexBarrel_layer1\n",
0327 "module1\n",
0328 "component0_0 is sensitive\n",
0329 "\tlocal offset [cm]: 0.0 0.0 0.0\n",
0330 "\tshape: TGeoBBox { 0.088375040, 13.500000, 0.0020000000 }\n",
0331 "\tgenerate ramdom hit:\n",
0332 "\t\tlocal coord ( with offset): 0.07891977012456126 3.2334660085878824 0.001675315310801563\n",
0333 "\t\tglobal coord (from local) : 3.6016753153108017 -0.07891977012456107 -3.2334660085878824\n",
0334 "\t\tCell ID: 0b11001010001000000000010011100000001000000000001000100011111\n",
0335 "\t\tglobal coord (from cellID): 3.6 -0.0779999999999998 -3.234\n",
0336 "VertexBarrel_layer2\n",
0337 "module1\n",
0338 "component0_0 is sensitive\n",
0339 "\tlocal offset [cm]: 0.0 0.0 0.0\n",
0340 "\tshape: TGeoBBox { 0.11783339, 13.500000, 0.0020000000 }\n",
0341 "\tgenerate ramdom hit:\n",
0342 "\t\tlocal coord ( with offset): 0.09520108098441747 9.405985024675164 0.0004767279434898437\n",
0343 "\t\tglobal coord (from local) : 4.800476727943489 -0.09520108098441689 -9.405985024675164\n",
0344 "\t\tCell ID: 0b1001001011111000000000011000000000001000000000001001000011111\n",
0345 "\t\tglobal coord (from cellID): 4.8 -0.09599999999999942 -9.406\n",
0346 "VertexBarrel_layer4\n",
0347 "module1\n",
0348 "component0_0 is sensitive\n",
0349 "\tlocal offset [cm]: 0.0 0.0 0.0\n",
0350 "\tshape: TGeoBBox { 0.29458347, 13.500000, 0.0020000000 }\n",
0351 "\tgenerate ramdom hit:\n",
0352 "\t\tlocal coord ( with offset): -0.037435451024441346 -13.085669240674665 0.0013062410037845283\n",
0353 "\t\tglobal coord (from local) : 12.001306241003785 0.03743545102444055 13.085669240674665\n",
0354 "\t\tCell ID: 0b1110011001110001111111111110110100000001000000000001010000011111\n",
0355 "\t\tglobal coord (from cellID): 12.0 0.0379999999999992 13.086\n"
0356 ]
0357 },
0358 {
0359 "name": "stderr",
0360 "output_type": "stream",
0361 "text": [
0362 "Warning in <TGeoMatrix::dtor>: Registered matrix component0_placement was removed\n",
0363 "Warning in <TGeoMatrix::dtor>: Registered matrix component0_placement was removed\n",
0364 "Warning in <TGeoMatrix::dtor>: Registered matrix component0_placement was removed\n"
0365 ]
0366 }
0367 ],
0368 "source": [
0369 "# detector_name = \"MiddleTrackerEndcapP\"\n",
0370 "detector_name = \"VertexBarrel\"\n",
0371 "dd=detector.detector(detector_name)\n",
0372 "\n",
0373 "## detector->layer->detector element (module)->placement->component->sensitive volume\n",
0374 "for dtemp, _ in dd.children():\n",
0375 " print(dtemp)\n",
0376 " ll=dd.child(dtemp) ## layer\n",
0377 " for ltemp,_ in ll.children(): \n",
0378 " print(ltemp)\n",
0379 " de=ll.child(ltemp) ## detector element (module)\n",
0380 " pl=de.placement() ## placement\n",
0381 " n_daughters = pl.ptr().GetNdaughters() ## number of components per de\n",
0382 " for ii in np.arange(n_daughters):\n",
0383 " comp = pl.daughter(int(ii)) ## component\n",
0384 " if comp.volume().isSensitive(): ## assume only one sensitive component per module\n",
0385 " print(comp.name(),\" is sensitive\")\n",
0386 " offset = comp.position() ## xyz offset of the component wrt the module center (local coordinate)\n",
0387 " print(\"\\tlocal offset [cm]: \", offset.x(),offset.y(),offset.z())\n",
0388 " vv=comp.volume().solid() ## component shape\n",
0389 " print(\"\\tshape: \", vv.type(), vv.dimensions())\n",
0390 " ## get a random point within the comp vol, then local --> global --> cellID\n",
0391 " print(\"\\tgenerate ramdom hit:\")\n",
0392 " x,y,z=random_point_in_volume(vv)\n",
0393 " local_pos = ROOT.dd4hep.Position(x,y,z)+offset\n",
0394 " print(\"\\t\\tlocal coord ( with offset): \", local_pos.x(),local_pos.y(),local_pos.z())\n",
0395 " global_pos = de.nominal().localToWorld(local_pos) ## cm\n",
0396 " print(\"\\t\\tglobal coord (from local) : \", global_pos.x(),global_pos.y(),global_pos.z())\n",
0397 " cid=converter.cellID(global_pos)\n",
0398 " print(\"\\t\\tCell ID: \", bin(cid))\n",
0399 " pp=converter.position(cid)\n",
0400 " print(\"\\t\\tglobal coord (from cellID): \", pp.x(),pp.y(),pp.z())\n",
0401 " ## tbd: loop over layers and pick random points wrt their rates \n",
0402 " break\n",
0403 "\n"
0404 ]
0405 },
0406 {
0407 "cell_type": "markdown",
0408 "id": "1abd06d7-f56f-497f-a0fa-d139643bda5c",
0409 "metadata": {},
0410 "source": [
0411 "## get segmentation index range for endcap trd2\n",
0412 "use de and offset from previous section"
0413 ]
0414 },
0415 {
0416 "cell_type": "code",
0417 "execution_count": 311,
0418 "id": "64c7adab-d74d-47e5-8648-4b07d8242769",
0419 "metadata": {},
0420 "outputs": [],
0421 "source": [
0422 "sens = detector.sensitiveDetector(detector_name)\n",
0423 "decoder = sens.readout().idSpec().decoder()\n",
0424 "# de.volumeID()\n",
0425 "# idDesc[\"system\"].value(volID)\n",
0426 "# decoder.fields()[5].name()\n",
0427 "# sens.readout().segmentation().segmentation().gridSizeX()\n"
0428 ]
0429 },
0430 {
0431 "cell_type": "code",
0432 "execution_count": 299,
0433 "id": "7152b9cc-a490-480c-bf3c-fb291fcb67c8",
0434 "metadata": {},
0435 "outputs": [],
0436 "source": [
0437 "yy = 0 ## ignore y dim for now\n",
0438 "lx=[]\n",
0439 "ly=[]\n",
0440 "lz=[]\n",
0441 "lsegx=[]\n",
0442 "lsegz=[]\n",
0443 "lc=[]\n",
0444 "for zz in np.arange(-z,z,0.1):\n",
0445 " # Linear interpolation factor\n",
0446 " dz = zz+z\n",
0447 " dx_z = (x2-x1)/2/z\n",
0448 " xlim = dx_z*dz+x1\n",
0449 " for xx in np.arange(-xlim,xlim,seg_x):\n",
0450 " local_pos = ROOT.dd4hep.Position(xx,yy,zz)+offset\n",
0451 " global_pos = de.nominal().localToWorld(local_pos)\n",
0452 " cid=converter.cellID(global_pos)\n",
0453 " lx.append(xx)\n",
0454 " ly.append(yy)\n",
0455 " lz.append(zz)\n",
0456 " lsegx.append(decoder[\"x\"].value(cid))\n",
0457 " lsegz.append(decoder[\"z\"].value(cid))\n",
0458 " lc.append(cid)\n",
0459 "df=pd.DataFrame({\"x\":lx,\"y\":ly,\"z\":lz,\"segx\":lsegx,\"segz\":lsegz,\"cell_id\":lc})"
0460 ]
0461 },
0462 {
0463 "cell_type": "code",
0464 "execution_count": 310,
0465 "id": "71d96363-54d5-465b-a5f4-ba7e388ade7c",
0466 "metadata": {},
0467 "outputs": [
0468 {
0469 "data": {
0470 "text/plain": [
0471 "Text(0, 0.5, 'seg z')"
0472 ]
0473 },
0474 "execution_count": 310,
0475 "metadata": {},
0476 "output_type": "execute_result"
0477 },
0478 {
0479 "data": {
0480 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGwCAYAAACaW3CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIHElEQVR4nO3dCXgUVdbw8RNIQhIgYQsBAkLAYTMCQ1AIiCBCwu7CgMiI4EAcBlDWVw0gICroKOA7ziAO4o7gOyp+OrIvMiphFZQoiA6ryqYiqIOs/T3nOt2mkxAS6O7qqvr/nqdIp/qmcutWhz59l1MRHo/HIwAAAAioUoE9HAAAAAiyAAAAgoSeLAAAgCAgyAIAAAgCgiwAAIAgIMgCAAAIAoIsAACAIIgMxkFxYefOnZOvv/5aypcvLxERETQZAAA2oOlFf/jhB6lRo4aUKlV0XxVBlkU0wKpVq5ZVvx4AAFyC/fv3S82aNYssQ5BlEe3B8l6k+Ph4q6oBAABK4Pjx46aTxPs+XiSPjaxZs8bTvXt3T/Xq1fVWQJ6FCxf6PX/u3DnPpEmTzPMxMTGedu3aeXJzc/3K/Pzzz57hw4d7Kleu7ImLi/P06NHDs3//fr8y3333nee2227zxMfHm00fHz161K/M3r17TV30GHqsu+66y3Py5Mlin8uxY8fMOehXAABgDyV5/7bVxPeffvpJmjZtKn/9618Lff7Pf/6zzJgxwzy/ceNGqVatmnTq1MmMnXqNHDlSFi5cKAsWLJD3339ffvzxR+nevbucPXvWV6Zfv36ydetWWbJkidn0cf/+/X3Pa9lu3bqZ+ugx9Fivv/66jBkzJsgtAAAAbMNjU/l7srQXq1q1ap5HHnnEr9cqISHBM3v2bPP9999/74mKivIsWLDAV+arr77ylCpVyrNkyRLz/aeffmqOvW7dOl+ZnJwcs2/Hjh3m+0WLFpmf0Z/1mj9/vqdMmTLF7pmiJwsAAPtxbE9WUXbv3i0HDx6UjIwM374yZcpIu3btZO3ateb7zZs3y+nTp/3K6OqA1NRUX5mcnBxJSEiQli1b+sq0atXK7MtbRn9Gf9YrMzNTTp48aX5HYfQ5HcfNuwEAAOdyTJClAZZKSkry26/fe5/Tr9HR0VKxYsUiy1StWrXA8XVf3jL5f48eU4/tLZPftGnTTKDm3VhZCACAszkmyPLKn3NKRxYvlIcqf5nCyl9Mmbyys7Pl2LFjvk1XFQIAAOdyTJClk9xV/p6kw4cP+3qdtMypU6fk6NGjRZY5dOhQgeMfOXLEr0z+36PH1KHI/D1ceYcuNVVD3g0AADiXY4KslJQUE/wsX77ct08DqjVr1kjr1q3N92lpaRIVFeVX5sCBA5Kbm+srk56ebnqaNmzY4Cuzfv16sy9vGf0Z/VmvZcuWmUBKfwcAAICtkpFquoUvvvjCb7K7pleoVKmSXHbZZSY9w9SpU+U3v/mN2fRxXFycScmgdC7UoEGDTKqFypUrm58bO3asXHnlldKxY0dTplGjRtK5c2fJysqSp59+2uy78847TZqHBg0amO914nzjxo1NWofHHntMvvvuO3Mc/Rl6qAAAgOGxkdWrV5tlk/m3AQMG+CUj1VQOmk7h2muv9Wzbts3vGCdOnDDJSCtVquSJjY01CUX37dvnV+bbb7/1/P73v/eUL1/ebPq4sGSk3bp1M8fQY+kxNWVEcZHCAQAA+ynJ+3eE/kO8GXqawkF71nQYMpC9X9/9eEpu/Msq2Xf81+SqAAC4UZnSETKk3eUyrMPlEh1ZKuTv3wRZDgqyrnpouRz58VRAjgUAgJP88doUye7aOKTv346Z+O52BFgAAJzf0//aLdMWfSqhRJDlADpESA8WAABFm/Pebjl15pyECkGWA/T9+y+3+wEAAOd3ziPyUs4eCRWCLAc4/APzsAAAKI693/1HQoUgywGqlo+2ugoAANhC7UpxIftdBFkOsODOXzLRAwCAovVPryOhQpDlAJXKRUvpou+BDQCA63VqWDVg+bKKgyDLIT68P8PqKgAAENZm394ipL+PIMshEuKiJDkhxupqAAAQlh69uYmULhXaYR+CLAf5IPt6q6sAAEDYiYksJbdcXSvkv5cgy2G2T+lsdRUAAAgrOx7qYsnvJchymNjo0lI9npQOAACo3MmZYhWCLAdaNbaD1VUAAMBy5aJLSbmYSMt+P0GWQ3uz6iWGLtkaAADhaN24Tpb+foIsh1o8op3VVQAAwDIVYiMt7cVSBFkOpcnWOjdOsroaAABYIie7o1iNIMvB/nZbmtVVAAAg5OomxpmpM1YjyHIwTbo2s08zq6sBAEBILQmTKTMEWQ53U/NkqVUp1upqAAAQElltU0J6f8KihEctEFTv3dNBorjSAACHu65Boozv1ljCBW+9LvHJFGuy3QIAEAq6kPC5O66WcEKQ5RLaddqpYVWrqwEAQFBsmRh+t5UjyHKR2be3sLoKAAAEXLX46LBYTZgfQZbLVhv+6ZoUq6sBAEBArQ7T28kRZLnM2K6NrK4CAAABkxImObEKQ5Dlwt6s6b2bWl0NAAACYmmY5MQqDEGWC/VKqymV46KtrgYAAI7JiVWY8K0ZgmrzRGvvTA4AwKVoVz+8cmIVhiDLxXInZ1pdBQAALsoLfwivnFiFIchysXIxkVIhNtLqagAAUCIfTcwQOyDIcrmc7I5WVwEAgGKLiYyQhLgosQOCLJfTZa91E+OsrgYAAMWyaYI9erEUQRZkSRgvfwUAwKtcdCkz1cUuHBVk1alTRyIiIgpsw4YNM88PHDiwwHOtWrXyO8bJkyflrrvukipVqkjZsmWlZ8+e8uWXX/qVOXr0qPTv318SEhLMpo+///57sStd/tqidgWrqwEAQJHWjbPXynhHBVkbN26UAwcO+Lbly5eb/b179/aV6dy5s1+ZRYsW+R1j5MiRsnDhQlmwYIG8//778uOPP0r37t3l7NmzvjL9+vWTrVu3ypIlS8ymjzXQsrNXstKtrgIAAOfVMKmcrXqxlL1qewGJiYl+3z/yyCNSr149adfu1+GwMmXKSLVq1Qr9+WPHjsncuXPlpZdeko4df5kQ/vLLL0utWrVkxYoVkpmZKdu3bzeB1bp166Rly5amzJw5cyQ9PV0+++wzadCggdi1N0uTus15b7fVVQEAoIAlo+w3tcVRPVl5nTp1ygRIf/jDH8ywoNe7774rVatWlfr160tWVpYcPnzY99zmzZvl9OnTkpHx66S6GjVqSGpqqqxdu9Z8n5OTY4YIvQGW0iFH3ectUxgdhjx+/LjfFm40qVuHhv6BKgAAVts+pbPYkWODrDfffNPMk9J5WF5dunSRefPmyapVq2T69OlmeLFDhw4mAFIHDx6U6OhoqVixot+xkpKSzHPeMhqk5af7vGUKM23aNN8cLt20dywcPTvwaokv49iXBQDAZuqG8Q2gL8Sx76Y67KdBlfZEed1yyy3SrVs30zPVo0cPWbx4sezcuVPeeeedIo/l8Xj8esPyPj5fmfyys7PNcKR3279/v4Sr9ePtszwWAOBsS2y8At6RQdbevXvNHKrBgwcXWa569epSu3Zt+fzzz833OldLhxl19WBeOqSovVneMocOHSpwrCNHjvjKFEbngsXHx/tt4Uo/MaSQOwsAYLHMRklhfQPoC7FvzYvw3HPPmeE77bUqyrfffmt6lDTYUmlpaRIVFeVblah0BWJubq60bt3afK8T3LUnasOGDb4y69evN/u8ZZxgqY0/OQAAnGFW/zSxM8cFWefOnTNB1oABAyQy8tfFk5qKYezYsWbi+p49e8wEeB0y1HxYN910kymjc6UGDRokY8aMkZUrV8qWLVvktttukyuvvNK32rBRo0YmDYROmtcVhrrpY03zYNeVhYXRTw6dGhacewYAQCjM7NNMSpc6/zQcO3BckKXDhPv27TOrCvMqXbq0bNu2TW644QazslCDMP2qQVf58uV95WbOnCk33nij9OnTR9q0aSNxcXHy9ttvm5/30snzGnjpKkTdmjRpYtI+OM3s21tYXQUAgAvVSIiRm5oni91FeHTGNkJOUzhoz5kOM4bz/Ky3P/pa7pq/xepqAABcZM8jRU/3scv7t+N6shBYPZrWkPb1yZ0FAAiN7TbNiVUYgixc0PN/uJoXCgAg6FJsnBOrMARZKJaPJ2fSUgCAoFrqsJXtBFkoFr0pZ7loXi4AgODo1LCqrXNiFcZZZ4OgWjeuEy0MAAiK2Q5c0U6QhRL1ZlWI/TX3GAAAgTC9d1Pb58QqDEEWSiQn+5ekrAAABEKF2CjplVbTkY1JkIUS0VUf15MJHgAQIFsnZTi2LQmyUGJzB14lieWiaTkAwCXJdfjKdYIsXJR14xg2BABcvHLRpcxcXycjyMJF0QmKw9rVo/UAABdlnQtWrBNk4aKNzmxA6wEASqxCbKTje7EUQRYuqTfrf/s2owUBACWS45KV6gRZuCQ3NEsmdxYAoNjqOez+hEUhyMIlc8snEgDApVvssPsTFoUgC5dMP5G0qVeZlgQAFGlAem3H3Z+wKO45UwTVvKxWtDAA4LzKlYmUB25IdVULEWQhYHY+1IXWBAAUKvcBZyceLQxBFgJGu4BvT69NiwIA/Gyf0lnciCALATXFZV3BAICipVSJdc1qwvwIshBwc3+fRqsCAIylI9uLWxFkIeDaX5FEqwIAJK12gqtWE+bn3jNHUDPBz+xDJngAcLv5Wa3FzQiyEBQ3NU+W5IQYWhcAXOp2l+XEKoy7zx5B9UH29bQwALhQTGQpFkIRZCHY3LpsFwDcbAd5Ew16shBUumy3bmIcrQwALpE72X1JR8+HIAtBt8RFNwMFADerEBsp5WIira5G2CDIQtDpxMfMhlVpaQBwuJzsjlZXIawQZCEkZt3egpYGAAdLSYxzbWb38yHIQshyZ03v3ZTWBgCHWsrUkAIIshAyvdJqShSvOABwHLdndj8fWgQhtWUiq04AwGncntn9fAiyEFK66iQ1OZ5WBwCHyGqbQi+WG4KsyZMnS0REhN9WrVo13/Mej8eUqVGjhsTGxkr79u3lk08+8TvGyZMn5a677pIqVapI2bJlpWfPnvLll1/6lTl69Kj0799fEhISzKaPv//++5Cdp9398662Qq8yANjfdQ0SZXy3xlZXI2w5KshSV1xxhRw4cMC3bdu2zffcn//8Z5kxY4b89a9/lY0bN5oArFOnTvLDDz/4yowcOVIWLlwoCxYskPfff19+/PFH6d69u5w9e9ZXpl+/frJ161ZZsmSJ2fSxBloovm2TyQQPAHYWXUrkuTuutroaYc1xGcMiIyP9eq/y9mI98cQTMn78eLn55pvNvhdeeEGSkpLklVdekT/+8Y9y7NgxmTt3rrz00kvSseMvuT5efvllqVWrlqxYsUIyMzNl+/btJrBat26dtGzZ0pSZM2eOpKeny2effSYNGjQotF7aQ6ab1/Hjx8XNvJngdx35j9VVAQBchA+ZY+u+nqzPP//cDAempKRI3759ZdeuXWb/7t275eDBg5KRkeErW6ZMGWnXrp2sXbvWfL9582Y5ffq0Xxk9Vmpqqq9MTk6OGSL0BliqVatWZp+3TGGmTZvmG17UTQM3tyMTPADYE5ndXRhkaeDz4osvytKlS03vkgZVrVu3lm+//dY8VtpzlZd+731Ov0ZHR0vFihWLLFO1asHs5brPW6Yw2dnZpqfMu+3fv1/czmSCb+R/PQAA4Y/M7i4cLuzSpYvv8ZVXXmmG8OrVq2eGBbW3Selk+PzDiPn35Ze/TGHlL3Qc7TXTDf5m9U+TeuMW0SwAYBNkdndpT1Z+ujpQgy0dQvTO08rf23T48GFf75aWOXXqlFk9WFSZQ4cOFfhdR44cKdBLhgsjEzwA2AuZ3YvP0UGWTjTXierVq1c3c7Q0QFq+fLnveQ2o1qxZY4YUVVpamkRFRfmV0RWKubm5vjLaO6bDfRs2bPCVWb9+vdnnLYOSZ4KvXDaaZgOAMEdOLBcHWWPHjjVBk05y18Dnd7/7nVnFN2DAADOUp+kZpk6dalI0aOA0cOBAiYuLMykZlE5IHzRokIwZM0ZWrlwpW7Zskdtuu830hnlXGzZq1Eg6d+4sWVlZZoWhbvpY0zycb2UhLmzz/Z1oJgAIY82SE8iJ5eY5WZo09NZbb5VvvvlGEhMTzTwsDYJq165tnr/nnnvkxIkTMnToUDMkqBPlly1bJuXLl/cdY+bMmSYNRJ8+fUzZ66+/Xp5//nkpXfrXO4vPmzdP7r77bt8qRE1Yqrm3cGm2T+ksjSYuoRkBIAy9PqyN1VWwnQiPzthGyGkPm/ac6TBjfDy3mfG6bvpq2U3uLAAIKzP7NJObmidbXQ3bvX87argQ9seESgAILzGREQRYF4kgC2GXOyutdoLV1QAA/NemCb8m6EbJEGQh7MzPYpUmAIQDMrtfGoIshGVv1u3pvyxWAABYh8zul4YgC2Fpyg2pEhPJyxMArHJ9w6oSG/3rynqUHO9iCFs7Hvr1NkkAgNBJLBctcwdeRZNfIoIshLXcyZlWVwEAXGfduF8ScOPSEGQhrJWLiTQTLwEAoTG0bV1zX1lcOoIshD0mXgJA6Izp0pDmDhCCLIQ9nXiZkhhndTUAwBWZ3enFChyCLNgCmeABILjiokuR2T3ACLJgC2SCB4Dg2kxm94AjyIJtkAkeAIKjTb3K5MQKAoIs2Ko3K6ttitXVAADHmZfVyuoqOBJBFmxlfLfG0qFhotXVAADH2D6ls9VVcCyCLNjOswOvljhSZwHAJWuWHM8wYRARZMGWNkwgEzwAXKrXh11DIwYRQRZsiUzwAHBppvduSk6sICPIgm2RCR4ALk5UKZFeaTVpviAjyIKtM8E3SY63uhoAYDtbJjLlIhQIsmBrC5lPAAAl0iCpnJlygeAjyIKt6T22ZvVrbnU1AMA2lo5qZ3UVXIMgC7bXtUl1yWyUZHU1ACDskRMrtAiy4Aiz+qdZXQUACGt1E+PIiRViBFlwzLDhzD7NrK4GAIStJSMYJgw1giw4xk3NkyUmMsLqagBA2NEpFXr/V4QWLQ5H2TQhw+oqAEDYYUqFNQiy4ChkggcAf2R2tw5BFhyHTPAA8Asyu1uLIAuOzATf9vIqVlcDACxHZndrEWTBkV4a3FKYAg/AzVKT48nsbjGCLDjWF1O7Wl0FALBEdGSE/POutrS+xQiy4OjcWd2bkAkegPvkTu5sdRXgtCBr2rRpctVVV0n58uWlatWqcuONN8pnn33mV2bgwIESERHht7Vq1cqvzMmTJ+Wuu+6SKlWqSNmyZaVnz57y5Zdf+pU5evSo9O/fXxISEsymj7///vuQnCeKb0Yf7msIwF3SaieQEytMOCrIWrNmjQwbNkzWrVsny5cvlzNnzkhGRob89NNPfuU6d+4sBw4c8G2LFi3ye37kyJGycOFCWbBggbz//vvy448/Svfu3eXs2bO+Mv369ZOtW7fKkiVLzKaPNdBCeNHke9zXEICbzM9qbXUV8F8RHo/HIw515MgR06Olwde1117r68nSHqc333yz0J85duyYJCYmyksvvSS33HKL2ff1119LrVq1TDCWmZkp27dvl8aNG5tgrmXLlqaMPk5PT5cdO3ZIgwYNLli348ePmx4w/X3x8fEBPW/4O3vOI/XG+QfSAOBE+qHy6QEtrK6Gox0vwfu3o3qy8tMGUJUqVfLb/+6775rgq379+pKVlSWHDx/2Pbd582Y5ffq06QHzqlGjhqSmpsratWvN9zk5OaaBvQGW0iFH3ectk58OQeqFybshNLivIQC3ILN7eHFskKUddKNHj5ZrrrnGBEheXbp0kXnz5smqVatk+vTpsnHjRunQoYMJgtTBgwclOjpaKlas6He8pKQk85y3jAZp+ek+b5nC5ot552/ppj1jCO19DWskxNDkAByLzO7hx7FB1vDhw+Xjjz+W+fPn++3XIcBu3bqZwKtHjx6yePFi2blzp7zzzjsXDNp0krxX3sfnK5NXdna26Vnzbvv377/oc8PFWZt9PU0HwJEql42WXmk1ra4G3BBk6crAt956S1avXi01axb9oqtevbrUrl1bPv/8c/N9tWrV5NSpU2b1YF46pKi9Wd4yhw4dKnQOmLdMfmXKlDFjt3k3hN72KSxrBuA8m+/vZHUV4PQgS3uStAfrjTfeMMOBKSkpF/yZb7/91vQqabCl0tLSJCoqyqxO9NIViLm5udK69S8rNnSCu/ZGbdiwwVdm/fr1Zp+3DML3ljspiXFWVwMAAoYPj+HLUUGWpm94+eWX5ZVXXjG5snR+lG4nTpwwz2sqhrFjx5qJ63v27DET4HXIUPNh3XTTTaaMzpcaNGiQjBkzRlauXClbtmyR2267Ta688krp2LGjKdOoUSOTBkInzeuqQt30saZ5KM7KQlhr6Yh2XAIAjqAfGvXDI8KTo4Ksp556yvQmtW/f3vRMebdXX33VPF+6dGnZtm2b3HDDDWZl4YABA8xXDbo0KPOaOXOmSWTap08fadOmjcTFxcnbb79tft5LJ89r4KWrEHVr0qSJSfsAm+TOalhw4QIA2A0fGsObo/NkhTPyZFmL3FkAnJDZ/fU/XWN1NVznOHmygAvnzhrati7NBMC2yOwe/hw1XAiUxJguDWkwALZ0e3pt7k9oAwRZcHVv1uzbuIE0AHuJiSwlU274Nck2whdBFlytc2p1mdmnmdXVAIBi2/FQF1rLJgiy4Hp6y52EGJZAAwh/uZMzra4CSoAgCxCRdePIlgwgvFWIjZRyMZFWVwMlQJAF/DcTfF0ywQMIYznZvyTEhn0QZAH/tYRM8ADCFJnd7YkgC/gvMsEDCFdkdrcngiwgj1m3t6A9AIQVvQ2YfgiE/XDVgHy5sx77XRPaBEDY4MOffRFkAfn0blGLFTwAwsKTt/7WfPiDPRFkAYUgFw0Aq13XIFF6NK1hdTVwCQiygPPYPqUzbQPAEpERIs/dcTWtb3MEWUARubOqx0fTPgBCbuskMrs7AUEWUIRVYzvQPgBCiszuzkGQBVygN6tZcjxtBCBkyOzuHARZwAW8Puwa2ghASJDZ3VkIsoAL0OXTj95M7iwAwUdmd2chyAKK4Zara0kMGZcBBFFW2xQyuzsMQRZQTDse6kJbAQiKZskJMr5bY1rXYQiygBIgSSmAYHh9WBsa1oFKHGT94Q9/kBdeeKHA/uPHj5vnACcrFxNpllcDQKDM7NOMW+c4VImDrOeff16GDh0qd999t5w7d863/8SJE4UGX4DTsLwaQKDERZeSm5on06AOdVHDhe+8844sXrxYMjMz5ejRo4GvFRDGyAQPIFA2T8igMR3sooKsxo0by7p16+T06dNy1VVXyfbt2wNfMyCMkQkewKXS23bphzY4V4mDrIiICPO1cuXKsmLFCmnfvr20atVK3nrrrWDUDwhL+h9j3cQ4q6sBwMb4sOZ8JZ7B6/F4fv3hyEh55plnTM+WztMC3GTJiHZSf8Jiq6sBwIbaXl6FXiwXKHGQtXr1aqlUqZLfvtGjR0uTJk3kgw8+CGTdgLAWHVnKJA+c895uq6sCwEZ0POilwS2trgZCIMKTt2sKIaMpLxISEuTYsWMSH88NiO1s4LMb5N2dR6yuBgCb+PfUrqRscMn7N8lIgUv0/B+uFqauAiiO6b2bEmC5CEEWEAAfTc6kHQEUKaqUSK+0mrSSixBkAQFAJngAF7JlIh/G3IYgCwgQMsEDOB+9HZd+GIO7lPiKf/zxx+fNnxUTEyOXXXaZlClTRtxi1qxZ8thjj8mBAwfkiiuukCeeeELatm1rdbVgUe6slMQ42X3kP7Q/AD98CHOnEgdZzZo18yUkLUxUVJTccsst8vTTT5ugy8leffVVGTlypAm02rRpY865S5cu8umnn5pgE+6zlNxZAPJpU68yObFcqsTDhQsXLpTf/OY38ve//122bt0qW7ZsMY8bNGggr7zyisydO1dWrVolEyZMEKebMWOGDBo0SAYPHiyNGjUyvVi1atWSp556yuqqwcLcWYOuSaH9AfjMy2pFa7hUiXuyHn74Yfnf//1fc3NoL01EWrNmTbn//vtlw4YNUrZsWRkzZow8/vjj4lSnTp2SzZs3y3333ee3PyMjQ9auXVug/MmTJ82WN88GnOn+7o1l5aeHZM93DBsCbrfzoS5WVwF26snatm2b1K5du8B+3afPeYcUdY6Sk33zzTdy9uxZSUpK8tuv3x88eLBA+WnTppnkZd5Ne7zgXCvHtre6CgAsltkoyfRuw71KfPUbNmwojzzyiOnJ8Tp9+rTZp8+pr776qkDw4VT556dpAv3C5qxlZ2eb7LDebf/+/SGsJUKtdKkIk3QQgHvN6p9mdRVgt+HCv/3tb9KzZ08zPKjDhBpQ6IpD7dX55z//acrs2rXL8TeMrlKlipQuXbpAr9Xhw4cLDTB1xaWbVl3il6SD2Qs/llNnuHMV4DZD29Ylszsu7t6FP/74o7z88suyc+dO03OjPVj9+vWT8uXLu6pJW7ZsKWlpaWZ1oVfjxo3lhhtuMMODReHehe5w4tRZaTRxidXVABBi3J/QuUry/n1RmdHKlSsnQ4YMEbcbPXq09O/fX1q0aCHp6elmleW+fftoG/jlzqoeHy0Hjv86vA7A2Wb2aUYvFoyLmpH30ksvyTXXXCM1atSQvXv3mn0zZ86U//f//p+4ieYD07QNU6ZMMZP9//Wvf8miRYsKXRgA91o1toPVVQAQIskJMXJT82TaGxcXZGkOKO3B0aSbR48eNXOxVMWKFU3A4TY692zPnj0mPYOmdLj22mutrhLCsDfr+oZVra4GgBD4IPt62hkXH2Q9+eSTMmfOHBk/frxERv462qhDZt4UDgD8zR14lSTEct8ywMnIiYVLDrJ2794tv/3tbwvs15VzP/30U0kPB7jGxvGdrK4CgCBJq51ATixcepCVkpJibqeT3+LFi83KOgCF06SEmpwQgPPMz2ptdRUQhko8fvE///M/MmzYMPn5559N+ga9jc78+fNNyoJnnnkmOLUEHJScsN64RVZXA0AAkdkdAQuy7rjjDjlz5ozcc8898p///Mfkx0pOTjb3M+zbt29JDwe4LhO8Lu8e9X8Fe4MB2BOZ3RHQZKR579937tw5qVqVlVMlRTJSd2s8cbH859Q5q6sB4BLp7bP07g5wj+MlSEZa4jlZJ06cMD1Y3lvL6PeaumHZsmUXX2PAZTZPyLC6CgAuUVSpX26fBQQsyNJbxrz44ovm8ffffy9XX321TJ8+3ezXHFoAipc7q029yjQVYGNbJmZaXQU4Lcj68MMPpW3btubxa6+9JtWqVTNZ3zXw+stf/hKMOgKONC+rldVVAHCRUpPjpVwMue8Q4CBLhwq9N4LWIcKbb75ZSpUqJa1atfLdYgdA8ZC8ELCf6MgI+eddv3Q2AAENsi6//HJ58803Zf/+/bJ06VLJyPhlbsnhw4cvOAEMQMHcWQPSudclYCe5kztbXQU4NciaOHGijB07VurUqSMtW7aU9PR0X69WYZngARTtgRtSaSLAJsjsjqCncDh48KAcOHBAmjZtaoYKlSYl1Z6shg0blvRwrkQKB+S1cttBGTRvM40C2GCIX3ug4V7HS5DC4aJm7elkd93y0lWGAC5O+yu43Q4Q7sjsjpIiHAfCKBM8gPBFZneUFEEWECZuap4sNRJirK4GgELohyD9MASUBEEWEEbWZl9vdRUA5FOrYqz5EASUFEEWEGa2T2F5OBAutO/qvXs7WF0N2BRBFhCGt9yplxhndTUAiMinfOjBJSDIAsLQ4hHtrK4C4HopiXHmQw9wsQiygDCkeXgyG1a1uhqAqy3lww4uEUEWEKZm3d7C6ioArqUfckg6iktFkAWEKV0uPrRtXaurAbgSH3IQCARZQBgb04XbVAGhph9uyImFQCDIAsKY/kf/v33JBA+EMmUDH24QKARZQJi7oVmypCYXfRNSAIHx1G3N6cVCwBBkATbwz7vaShQryYGg3zqnc2p1WhkBQ5AF2MTHk8gEDwRLQkxpbp2DgCPIAmxCkyLWJRM8EBTrxnWiZRFwBFmAjSwhOSIQcPrhhczuCAaCLMBGyAQPBB4fXhAsBFmAzZAkEQgcMrsjmAiyABvmzpreu6nV1QAcgQ8tCCaCLMCGeqXVlMrloq2uBmBrs/qREwvB5Zgga8+ePTJo0CBJSUmR2NhYqVevnkyaNElOnTrlVy4iIqLANnv2bL8y27Ztk3bt2pnjJCcny5QpU8Tj8fiVWbNmjaSlpUlMTIzUrVu3wDGAYNs8oZPJTg2g5O5oXUe6NiEnFoIrUhxix44dcu7cOXn66afl8ssvl9zcXMnKypKffvpJHn/8cb+yzz33nHTu/GvOoYSEBN/j48ePS6dOneS6666TjRs3ys6dO2XgwIFStmxZGTNmjCmze/du6dq1qzn+yy+/LB988IEMHTpUEhMTpVevXiE8a7jdp1M6S6OJS6yuBmArsVERMqnnFVZXAy7gmCBLg6a8gZP2Ln322Wfy1FNPFQiyKlSoINWqVSv0OPPmzZOff/5Znn/+eSlTpoykpqaaQGvGjBkyevRoX8/XZZddJk888YT5mUaNGsmmTZvM7zlfkHXy5Emz5Q3mgEuly86rx0fLgeP+PbYAzu/D+zNpHoSEY4YLC3Ps2DGpVKlSgf3Dhw+XKlWqyFVXXWUCJu0B88rJyTFDhRpgeWVmZsrXX39thiS9ZTIyMvyOqWU00Dp9+nShdZk2bZrpMfNutWrVCuCZws1Wje1gdRUA29APJeTEQqg4Nsj697//LU8++aQMGTLEb/+DDz4o//jHP2TFihXSt29fMwQ4depU3/MHDx6UpKQkv5/xfq/PFVXmzJkz8s033xRan+zsbBP0ebf9+/cH7FzhbvqG0YwbSAPFwocShFLYB1mTJ08udLJ63k17kPLSXicdOuzdu7cMHjzY77kJEyZIenq6NGvWzARYOqn9scce8yujx8zLO+k97/7ilMlLe8bi4+P9NiBQXh92DY0JXACZ3RFqYT8nS4f2tMepKHXq1PELsHTSugZSf//73y94/FatWpn5UYcOHTK9UTpXy9tj5XX48GHz1dt7db4ykZGRUrly5RKdHxCo3FmP/a6J/M9rH9OgwHmQ2R2hFvZBls6d0q04vvrqKxNgaWoFXUFYqtSFO+q2bNli0jDoZHilwdm4ceNM6ofo6F/yEC1btkxq1KjhC+a0zNtvv+13HC3TokULiYqKuoizBC5d7xa15IG3PpUfT52hOYF8stqmmNtSAaHkmFec9mC1b9/eTCjXVX5HjhwxvU15e5w0MJozZ45J76Bztp555hkZP3683Hnnnb6J7v369TOPNW2Dllu4cKGZs+VdWah0ntfevXvNvu3bt8uzzz4rc+fOlbFjx1p2/oDKncKqKSC/9vUTZXy3xjQMQi7Ckz/Lpk1pyoU77rij0Oe8p7hkyRIzAf2LL74wKwo1zYPO2Ro2bJgZ6subjFT3bdiwQSpWrGiCqokTJ/rNt9JkpKNGjZJPPvnE9HLde++9BSbZF0WHKHWVoU6CZ34WAunHn89I6uSlNCqgQ+m6EOqRbrQFAqYk79+OCbLshiALwXT11OVymNxZgOROzpRyMWE/MwYOff92zHAhgF+9f8/1NAdcr0JsJAEWLEWQBTiQTvBNq/3r7aIAN8rJ7mh1FeByBFmAQ83Pam11FQDLkNkd4YAgC3Bwb1ZmI/87EwBuQWZ3hAOCLMDBZvVPs7oKQMi1vbwK9ydEWCDIAhyeCX72bc2trgYQMppo56XBLWlxhAWCLMDhOqdWl5l9mlldDSAkvpjalZZG2CDIAlzgpubJUjaq8JuXA07xp2tSTO8tEC4IsgCXWD8+w+oqAEE1tmsjWhhhhSALcAnNeq3JGQEnGtq2Lr1YCDsEWYCLkJwRTjWmS0OrqwAUQJAFuEhsdGlpllz0vbYAu9GFHczFQjgiyAJc5vVh11hdBSBgkhNizMIOIBwRZAEuo5/4n7z1t1ZXAwiID7K5GTrCF0EW4EI9mtaQZsncQBr2tn1KZ6urABSJIAtwqdeHtbG6CsBFq5sYx61zEPYIsgAXDxtO793U6moAF2XJiHa0HMIeQRbgYr3SakoU/wvAZjIbVpXoSF64CH+8SgGX2zIx0+oqACUy6/YWtBhsgSALcDkywcNOyOwOOyHIAkAmeNgGmd1hJwRZAMwqrbaXV6YlENbI7A67IcgCYLw0uBUtgbBVq1Ismd1hOwRZAHx2PtSF1kDYKR0h8t49HayuBlBiBFkAfHRZfGajJFoEYSX3ATK7w54IsgD4mdU/jRZB2CCzO+yMIAtAgUzwOsEYCAdkdoedEWQBKOCm5skSExlBy8BSZHaH3RFkASjUpgkZtAwsRWZ32B1BFoDzZoJvmFSO1oEl9OblOnQN2BlBFoDzWjKqHa2DkKtcNtrcvBywO4IsAEUidxZCbfP9nWh0OAJBFoAL585qWJVWQkhsn0JOLDiHo4KsOnXqSEREhN923333+ZXZt2+f9OjRQ8qWLStVqlSRu+++W06dOuVXZtu2bdKuXTuJjY2V5ORkmTJling8Hr8ya9askbS0NImJiZG6devK7NmzQ3KOgBWYgIxQqB4fbe6jCThFpDiMBkRZWVm+78uV+3Xi7tmzZ6Vbt26SmJgo77//vnz77bcyYMAAE0A9+eSTpszx48elU6dOct1118nGjRtl586dMnDgQBOUjRkzxpTZvXu3dO3a1fyel19+WT744AMZOnSoOW6vXr0sOGsguHQC8tC2dWXWe7toagTNqrHcOgfO4rggq3z58lKtWrVCn1u2bJl8+umnsn//fqlRo4bZN336dBNEPfzwwxIfHy/z5s2Tn3/+WZ5//nkpU6aMpKammkBrxowZMnr0aNM7pr1Wl112mTzxxBPmGI0aNZJNmzbJ448/TpAFxxrTpSFBFoKGzO5wIkcNF6pHH31UKleuLM2aNTOBU96hwJycHBM0eQMslZmZKSdPnpTNmzf7yuhQoQZYect8/fXXsmfPHl+ZjAz/HEJaRgOt06dPF1ov/R3aS5Z3A+yETPAIJjK7w4kcFWSNGDFCFixYIKtXr5bhw4ebniYdxvM6ePCgJCX53/y2YsWKEh0dbZ47Xxnv9xcqc+bMGfnmm28Krdu0adMkISHBt9WqVStAZw2ENhN8jYQYmhwBdXt6bbPAAnCasH9VT548ucBk9vyb9iCpUaNGmV6oJk2ayODBg82w3ty5c83cKy8tn5/Oycq7P38Z76T3kpbJKzs7W44dO+bbdMgSsKO12ddbXQU4SExkKZlyQ6rV1QDcOSdLe6T69u17wVWFhWnVqpX5+sUXX5ghRJ2rtX79er8yR48eNUN83p4pLePtsfI6fPiw+XqhMpGRkeb3FEaHH/MOQQJ2z51Vf8Jiq6sBB9jxUBerqwC4N8jSNAu6XYwtW7aYr9WrVzdf09PTzTytAwcO+PbpZHgNfjQdg7fMuHHjzFwuHUb0ltF5XN5gTsu8/fbbfr9Ly7Ro0UKioqIu4WwBe9ChnbTaCbJ57zGrqwIby52caXUVAHcPFxaXTkafOXOmbN261aRY+L//+z/54x//KD179jQrAZVOVm/cuLH079/fBGArV66UsWPHmlQMurJQ9evXzwRduuIwNzdXFi5cKFOnTvWtLFRDhgyRvXv3mn3bt2+XZ5991gxL6rEAt5if1drqKsDGKsRGmvtjAk7mmCBLA6NXX31V2rdvbwKpiRMnmuBp/vz5vjKlS5eWd955xyQQbdOmjfTp00duvPFGk3rBSyelL1++XL788kvTM6UT5zWY0s0rJSVFFi1aJO+++65Zxfjggw/KX/7yF9I3wFXIBI9LkZPdkQaE40V48qcyR0hoCgcN6HQSvLcXDbCbs+c8Um/cIqurAZtJSYyT1WOus7oaQNDfvx3TkwXAukzwQEksHdGOBoMrEGQBuORM8EBx6YIJcmLBLQiyAFwSMsGjJFgwATchyAIQkEzwtSrF0pIoUlbbFHqx4CoEWQAC4r17Okjpwm94AEj7+okyvltjWgKuQpAFIGByH+hMa6KA0iLy/B+upmXgOgRZAAImNrq01E2Mo0Xh5yMyu8OlCLIABNQSlucjDzK7w80IsgAEFJngkReZ3eFmBFkAAm7W7S1oVZjM7jqEDLgVQRaAoOTOmt67KS3rcmR2h9sRZAEIil5pNaVy2Wha16UGXUNOLIAgC0DQbL6/E63rQnUqxcn93cmJBRBkAQiq7VPIneU2K8e2t7oKQFggyAIQVDrxWSdAwx1m9mlm5uQBIMgCEAJMgHaHmMgIcx9LAL+gJwtASHJntahdgZZ2uE0TMqyuAhBWCLIAhMQrWem0tIOR2R0oiCALQMh6s+5oXYfWdigyuwMFEWQBCJlJPa+QCrFRtLjDtL28MpndgUIQZAEIqa2TmLfjNC8NbmV1FYCwRJAFIOTIneUc/57a1eoqAGGLIAuAJbmzqsVzyx27G9q2LjmxgCIQZAGwxOqxHWh5mxvTpaHVVQDCGkEWAMt6s+qSCd62yOwOXBhBFgDLLBnRjta3ITK7A8VDkAXA0txZmQ2rcgVshszuQPEQZAGw1KzbW3AFbKRhUjkpFxNpdTUAWyDIAmCp0qUiZFa/5lwFm1gyiiFeoLgIsgBYrmuT6jIgvbbV1cAF7HyoC20ElABBFoCw8MANqRIZYXUtcD46d07n0AEoPv5iAISNrZMyra4CzoO5c0DJEWQBCBs6obpCLJOqww2Z3YGLQ5AFIKzkZHe0ugrIh8zugMuDrHfffVciIiIK3TZu3OgrV9jzs2fP9jvWtm3bpF27dhIbGyvJyckyZcoU8Xg8fmXWrFkjaWlpEhMTI3Xr1i1wDAAXnwk+hUzwYYPM7sDFc0y/fOvWreXAgQN+++6//35ZsWKFtGjhn4fnueeek86dO/u+T0hI8D0+fvy4dOrUSa677joTnO3cuVMGDhwoZcuWlTFjxpgyu3fvlq5du0pWVpa8/PLL8sEHH8jQoUMlMTFRevXqFfRzBZxu6Yh2Un/CYqur4XrJCTFyU/Nk17cDIG4PsqKjo6VatWq+70+fPi1vvfWWDB8+3PRW5VWhQgW/snnNmzdPfv75Z3n++eelTJkykpqaagKtGTNmyOjRo309X5dddpk88cQT5mcaNWokmzZtkscff5wgCwjE33NkKRl0TYrMfX837WmhD7Kvp/2BS+CY4cL8NMD65ptvTC9Ufhp4ValSRa666ioTMJ07d873XE5Ojhkq1ADLKzMzU77++mvZs2ePr0xGRobfMbWMBloa3BXm5MmTppcs7wbg/O7v3ljqVIqjiSyyfcqvvf0ALo5jg6y5c+eawKdWrVp++x988EH5xz/+YYYR+/bta4YAp06d6nv+4MGDkpSU5Pcz3u/1uaLKnDlzxgR2hZk2bZoZlvRu+esFoKCVY9vTLBaomxhn5sYBcHiQNXny5PNOaPdu2oOU15dffilLly6VQYMGFTjehAkTJD09XZo1a2YCLJ3U/thjj/mVyT+86J30nnd/ccrklZ2dLceOHfNt+/fvL3FbAG685Y5OvEZoLRnBrXMAV8zJ0qE97XEqSp06dQpMbK9cubL07Nnzgsdv1aqVGbo7dOiQ6Y3SuVreHiuvw4cPm6/e3qvzlYmMjDS/tzA6/Jh3CBJA8ejE6+w3PpKfz/iv8EVwkNkdcFGQpXOndCsu7VHSIOv222+XqKioC5bfsmWLScOgk+GV9nKNGzdOTp06ZSbTq2XLlkmNGjV8wZyWefvtt/2Oo2V0FWNxfieAktk0IUNSJy+l2UKAzO6Ai4YLS2rVqlUmxUJhQ4UaGM2ZM0dyc3Pl3//+tzzzzDMyfvx4ufPOO329TP369TOPdcK8llu4cKGZs+VdWaiGDBkie/fuNfu2b98uzz77rJkDNnbs2JCfL+AGZIIPjem9m5ohWgCBEeHJn2XT5jRI0gBIc1flt2TJEjM36osvvjArCjWJ6ODBg2XYsGFmqC9vMlLdt2HDBqlYsaIJqiZOnOg330qTkY4aNUo++eQT08t17733mnLFpUOUOgFe52fFx8cH4MwBZztx6qw0mrjE6mo4VlQpkc+ndrO6GkDYK8n7t+OCLLsgyAJKrv8z6+W9LwpfwYtLkzs50/QYAgjc+7fjhgsBONdLg1vyn1YQpCbHE2ABQUCQBcBWdj3CkFYgRZYS+eddbQN6TAC/IMgCYDsfTfS/4wIu3rbJZHYHgoUgC4DtJMRFSVw0/31dKjK7A8HF/1IAbGnzBHqzLhWZ3YHgIsgCYEt6b71KsayGu1iZjZIkWidkAQga/sIA2Na/7r3e6irY1qz+aVZXAXA8giwAtqV5nRoklbO6GrZDZncgNAiyANja0lHtrK6CrVQuGy290mpaXQ3AFQiyANje9imkISiuzfd3Cuq1APArgiwAjpgEn5IYZ3U1wh7BKBBaBFkAHGHpCIYNi6JBqAajAEKHIAuAI2g6grTaCVZXI2wRhAKhR5AFwDHmZ7W2ugphSYNPcmIBoUeQBcAxNJDQJJvwR/AJWIMgC4CjkGTTH5ndAesQZAFwlNKlImRmn2ZWVyNsEHQC1iHIAuA4NzVPlloVY8XtZt/W3ASdAKxBkAXAkd67t4O4ObzQ3rzOqdWtrgbgagRZABzrU5dmgk+IKW168wBYiyALgGNp8s0myfHiNuvGcescIBwQZAFwtIXDrhE3IbM7ED4IsgA4mk78/tM1KeIWZHYHwgdBFgDHG9u1kbgBmd2B8EKQBcAVvVmP/a6JOB2Z3YHwQpAFwBV6t6gl5cpEilNltU3h/oRAmCHIAuAauQ9kihNd1yBRxndrbHU1AORDkAXAVXInOyvQii4l8twdV1tdDQCFIMgC4CrlYiKlQqxzhg0/nOisoBFwEoIsAK6Tk91RnKBx9XgTNAIITwRZAFyZCb5T46piZ3rf50Uj2lpdDQBFIMgC4Epzbr/KtoFWXKTIrmndrK4GgAugnxmAqwOtE6fOyrg3PpSFWw9LuKseX0beGt5WEuPLWF0VAMUQ4fF4PMUpiMA6fvy4JCQkyLFjxyQ+3n03sAUAwOnv37YZLnz44YeldevWEhcXJxUqVCi0zL59+6RHjx5StmxZqVKlitx9991y6tQpvzLbtm2Tdu3aSWxsrCQnJ8uUKVMkf5y5Zs0aSUtLk5iYGKlbt67Mnj27wO96/fXXpXHjxlKmTBnzdeHChQE+YwAAYGe2CbI0WOrdu7f86U9/KvT5s2fPSrdu3eSnn36S999/XxYsWGACoTFjxvhFn506dZIaNWrIxo0b5cknn5THH39cZsyY4Suze/du6dq1q7Rt21a2bNki48aNM8GaHssrJydHbrnlFunfv7989NFH5mufPn1k/fr1QW4FAABgF7YbLnz++edl5MiR8v333/vtX7x4sXTv3l32799vgiilgdbAgQPl8OHDpkvvqaeekuzsbDl06JDpgVKPPPKICba+/PJLiYiIkHvvvVfeeust2b59u+/YQ4YMMcGUBldKAywN2PR3enXu3FkqVqwo8+fPL7TeJ0+eNJuX/nytWrUYLgQAwEYcOVx4IRoApaam+gIslZmZaQKbzZs3+8roUKE3wPKW+frrr2XPnj2+MhkZGX7H1jKbNm2S06dPF1lm7dq1563ftGnTzEXxbhpgAQAA53JMkHXw4EFJSkry26c9S9HR0ea585Xxfn+hMmfOnJFvvvmmyDLeYxRGe9A06vVu2uMGAACcy9Iga/LkyWaIrqhNe5CKS8vnp6OheffnL+MdLQ1EmcJ+v5f2nmm3Yt4NAAA4l6V5soYPHy59+/YtskydOnWKdaxq1aoVmHh+9OhRM8Tn7XXSMvl7m3S+lrpQmcjISKlcuXKRZfL3bgEAAPeyNMjSNAu6BUJ6erpJ83DgwAGpXr262bds2TLTg6TpGLxldLWgrlTUYURvGZ3H5Q3mtMzbb7/td2wt06JFC4mKivKVWb58uYwaNcqvjKaYAAAAsNWcLM2BtXXrVvNV0zXoY91+/PFH87xORNd8VZpOQVMvrFy5UsaOHStZWVm+obl+/fqZoEtXHObm5prcVlOnTpXRo0f7hvp0JeHevXvNPl1h+Oyzz8rcuXPNsbxGjBhhgqpHH31UduzYYb6uWLHCrHoEAACwVQoHDYxeeOGFAvtXr14t7du3N481ABs6dKisWrXKJBvVoErzYOVdTajJSIcNGyYbNmwwE+M1qJo4caLffCpNRqq9VJ988onp5dK0Dlour9dee00mTJggu3btknr16pletJtvvrnY56OT3zWpqk6AZ34WAAD24E3BpKmkNFuAI4Isp9G8XKRxAADAnrSTpGbNmkWWIciyyLlz50x+rvLlyxe5KtGJ0b/beu/cet5uPne3nrebz92t5+3Gc/d4PPLDDz+Yka5SpUqF78R3N9MLc6EI2KncmsLCreft5nN363m7+dzdet5uO/eECwwT2m7iOwAAgJ0QZAEAAAQBQRZCRld5Tpo0yW+1pxu49bzdfO5uPW83n7tbz9vt534hTHwHAAAIAnqyAAAAgoAgCwAAIAgIsgAAAIKAIAsAACAICLJwyfS+ja1bt5a4uDhzP8bCaFb7/Nvs2bP9yuh9Jdu1a2fuO5mcnCxTpkwxmXXz0vtKpqWlSUxMjNStW7fAMcLtvPV+mj169JCyZctKlSpV5O6775ZTp07Z+rwLU6dOnQLX97777gtKW9jBrFmzJCUlxVwvvW7vvfee2NnkyZMLXN9q1ar5ntdrpGU0A7ZeO72frN77Na+TJ0/KXXfdZa69vgZ69uxpbi8WTv71r3+Z16ieh57jm2++6fd8oM7z6NGj0r9/f5PQUjd9rPfBC9fz1nsH57/+rVq1sv15h4TeuxC4FBMnTvTMmDHDM3r0aE9CQkKhZfSl9txzz3kOHDjg2/7zn//4nj927JgnKSnJ07dvX8+2bds8r7/+uqd8+fKexx9/3Fdm165dnri4OM+IESM8n376qWfOnDmeqKgoz2uvvRaW533mzBlPamqq57rrrvN8+OGHnuXLl3tq1KjhGT58uK3PuzC1a9f2TJkyxe/6/vDDDwFvCztYsGCBuT56nfR66XUrW7asZ+/evR67mjRpkueKK67wu76HDx/2Pf/II4+Ya6XXTK/dLbfc4qlevbrn+PHjvjJDhgzxJCcnm2uvrwF9LTRt2tS8NsLFokWLPOPHjzfnof9nLVy40O/5QJ1n586dzd/D2rVrzaaPu3fv7gnX8x4wYICpc97r/+233/qVseN5hwJBFgJGg6iigqz8f7h5zZo1y/zszz//7Ns3bdo080Z87tw58/0999zjadiwod/P/fGPf/S0atXKE47nrf9xlSpVyvPVV1/59s2fP99TpkwZE1DY/bzzB1kzZ8487/OBags7uPrqq80bTl56/e677z6PnYMsfcMsjF6batWqmQDES6+hXsvZs2eb77///nsTeGoA6qWvBX1NLFmyxBOO8v+fFajz1MBbj71u3TpfmZycHLNvx44dHqudL8i64YYbzvszTjjvYGG4ECEzfPhw05V81VVXmeEuvUm2V05OjhkmypvMLjMz09xEe8+ePb4yGRkZfsfUMps2bZLTp0+H3ZXU+qamppou+Lz11W71zZs3O+68H330UalcubI0a9bMDKXmHQoMVFuEOz1nPZ/810u/X7t2rdjZ559/bq6fDoP27dtXdu3aZfbv3r1bDh486HfOeg31WnrPWdtEX6t5y+ix9DVhl3YJ1Hnq61yHylq2bOkro0Nvui+c2+Ldd9+VqlWrSv369SUrK0sOHz7se87J532pCLIQEg8++KD84x//kBUrVpj/oMeMGSNTp071Pa//eSUlJfn9jPd7fa6oMmfOnJFvvvlGwk1h9a1YsaJER0df8Jy8z9nlvEeMGCELFiyQ1atXm2D6iSeekKFDhwa8LcKdXo+zZ88Weh52OYfC6Bvjiy++KEuXLpU5c+aYc9H5iN9++63vvIo6Z/2q11qv+fnKhLtAnad+1WAlP90Xrm3RpUsXmTdvnqxatUqmT58uGzdulA4dOpgPSU4+70AgyEKxJ7rm37QnpbgmTJgg6enpppdDAyyd0PzYY4/5ldFj5uWd8Jx3f3HKhNN5F1YvrXNJzynY532pbTFq1Cjzib5JkyYyePBg01M5d+5c8yZ8vnPwnofV5xkMhZ2H3c4h/5tsr1695Morr5SOHTvKO++8Y/a/8MILl3TOdmyXQJxncf4Wwsktt9wi3bp1Mz1TOkF+8eLFsnPnTt/rwKnnHQiRATkKHEd7I7TH6UIryi6WdhMfP35cDh06ZD7t6Eql/J9mvN3R3k+O5ysTGRlphqnC7by1vuvXry+wuka71S90TqE+70C3hXfl0RdffGHqGKi2CHc6HF66dOlCz8Mu51AcunpMAy4dQrzxxhvNPj3n6tWrF3rOem11KFWved7eDi2jPWJ24F1NeannqWX0/738jhw5YpvXiJ5/7dq1zfV303lfDHqycN43i4YNGxa56fL0i7Vlyxbz897UB9rLpcuI887jWbZsmRnX976Ra5nly5f7HUfLtGjRQqKiosLuvLW+ubm5cuDAAb/66jwOXdYfTucd6LbQ66u8b0aBaotwp0Mmej75r5d+b5dgojh0mGj79u3m+uocLX0DzXvOeg017Yj3nLVN9LWat4y+FvQ1YZd2CdR56uv82LFjsmHDBl8Z/QCi++zSFtpDvX//ft/ft1vO+6IEbUo9XEOXpm/ZssXzwAMPeMqVK2ce6+Zdwv/WW295/v73v5slz1988YVZ2h4fH++5++67/Van6PL9W2+91ZR74403TJnCUhmMGjXKrFSZO3eupakMLnTe3rQF119/vVnSvGLFCk/NmjX90hbY8bzz06XYmspCz13r+uqrr5oVgT179vSVCVRb2CmFg14nvV4jR440KRz27NnjsasxY8Z43n33XXN9dXWYLrvXVAbec9IVd7rKTq+ZXju9hoWlNtBrrtdeXwMdOnQIuxQO+rfr/TvWt0fv69qbfiNQ56mpDJo0aWJW1+l25ZVXWprKoKjz1uf0+uvf+e7duz2rV6/2pKenm3QNdj/vUCDIwiXT5b36h5l/0z9GtXjxYk+zZs1MIKLBgr7ZPvHEE57Tp0/7Hefjjz/2tG3b1izr16XSkydPLrB0X/+j/+1vf+uJjo721KlTx/PUU0+F7Xkr/U+qW7duntjYWE+lSpVMUJE3RYEdzzu/zZs3e1q2bGnefGJiYjwNGjQwS/5/+uknv3KBags7+Nvf/mbSWuj1at68uWfNmjUeO/Pmg9LgUQPom2++2fPJJ5/4ntdrpNdcr5leu2uvvdYEIXmdOHHCXHO99voa0DfXffv2ecKJ/u0W9jetf+uBPE/NMfX73//eBKq66eOjR496wvG8NZ9hRkaGJzEx0Vz/yy67zOzPf052PO9QiNB/Lq4PDAAAAOfDnCwAAIAgIMgCAAAIAoIsAACAICDIAgAACAKCLAAAgCAgyAIAAAgCgiwAAIAgIMgCAAAIAoIsAACAICDIAgAACAKCLAAAgCAgyAKAfF577TW58sorJTY2VipXriwdO3aUn376yff8c889J40aNZKYmBhp2LChzJo1y+/n165dK82aNTPPt2jRQt58802JiIiQrVu3FtrWO3bskLi4OHnllVd8+9544w3z89u2beP6ADYVaXUFACCcHDhwQG699Vb585//LDfddJP88MMP8t5774nH4zHPz5kzRyZNmiR//etf5be//a1s2bJFsrKypGzZsjJgwABTvkePHtK1a1cTNO3du1dGjhxZ5O/UQO3xxx+XoUOHSps2bSQqKsoc85FHHjHBHgB7ivB4/+cAAMiHH34oaWlpsmfPHqldu3aBFrnsssvk0UcfNYGY10MPPSSLFi0yPVizZ8+WCRMmyJdffml6otQzzzxjgiYNyLSH63y6d+8ux48fl+joaClVqpQsXbrU9IABsCd6sgAgj6ZNm8r1119vepAyMzMlIyNDfve730nFihXlyJEjsn//fhk0aJAJmrzOnDkjCQkJ5vFnn30mTZo08QVY6uqrry5WGz/77LNSv359E2Dl5uYSYAE2R5AFAHmULl1ali9fbnqlli1bJk8++aSMHz9e1q9fb+ZNeYcMW7ZsWeDnlA4O5O99Ku6AwUcffWTmfmmQdfDgQalRowbXBrAxJr4DQD4aJOncqAceeMAM8enw3cKFCyUpKUmSk5Nl165dcvnll/ttKSkpvvlVH3/8sZw8edJ3vE2bNl2wjb/77jsZOHCgCejuuOMO+f3vfy8nTpzg2gA2Rk8WAOShPVYrV640w4RVq1Y13+swoa4mVJMnT5a7775b4uPjpUuXLiaY0iDq6NGjMnr0aOnXr58JlO6880657777ZN++fWZSuypqftWQIUOkVq1aZj7XqVOnpHnz5jJ27Fj529/+xvUB7EonvgMAfvHpp596MjMzPYmJiZ4yZcp46tev73nyySf9mmfevHmeZs2aeaKjoz0VK1b0XHvttZ433njD9/wHH3zgadKkiXk+LS3N88orr+h4oWfHjh2FNvMLL7zgKVu2rGfnzp2+fZs2bTI//84773BpAJtidSEABNm8efPMEOCxY8dM7i0A7sBwIQAE2Isvvih169Y187d0Mvu9994rffr0IcACXIYgCwACTFcGTpw40XytXr269O7dWx5++GHaGXAZhgsBAACCgBQOAAAAQUCQBQAAEAQEWQAAAEFAkAUAABAEBFkAAABBQJAFAAAQBARZAAAAQUCQBQAAIIH3/wEMOBzLj1i64wAAAABJRU5ErkJggg==",
0481 "text/plain": [
0482 "<Figure size 640x480 with 1 Axes>"
0483 ]
0484 },
0485 "metadata": {},
0486 "output_type": "display_data"
0487 }
0488 ],
0489 "source": [
0490 "plt.scatter(df.segx,df.segz)\n",
0491 "plt.xlabel(\"seg x\")\n",
0492 "plt.ylabel(\"seg z\")"
0493 ]
0494 },
0495 {
0496 "cell_type": "markdown",
0497 "id": "9061f942-1a48-47fd-8172-428d49a522f7",
0498 "metadata": {},
0499 "source": [
0500 "## get neighbour cell ID for a given cell"
0501 ]
0502 },
0503 {
0504 "cell_type": "code",
0505 "execution_count": 4,
0506 "id": "b1311639-555a-4663-a035-01bee8fd5f31",
0507 "metadata": {},
0508 "outputs": [
0509 {
0510 "name": "stdout",
0511 "output_type": "stream",
0512 "text": [
0513 "Cell ID 16845714603764658719, Position: 4.276590568660096 2.181567580444739 11.376\n",
0514 "Cell ID 16845996074446402079, Position: 4.275647775186444 2.1833314229734357 11.374\n",
0515 "Cell ID 16845996083036336671, Position: 4.277533362133748 2.1798037379160426 11.374\n",
0516 "Cell ID 16846277553718080031, Position: 4.276590568660096 2.181567580444739 11.372\n"
0517 ]
0518 }
0519 ],
0520 "source": [
0521 "\n",
0522 "sens = detector.sensitiveDetector(detector_name)\n",
0523 "neighbours = ROOT.std.set(\"dd4hep::CellID\")()\n",
0524 "sens.readout().segmentation().neighbours(cellID, neighbours)\n",
0525 "\n",
0526 "converter = DDRec.CellIDPositionConverter(detector)\n",
0527 "\n",
0528 "for nn in neighbours:\n",
0529 "\tpos=converter.position(nn)\n",
0530 "\tprint(f\"Cell ID {nn}, Position:\",pos.x(),pos.y(),pos.z())"
0531 ]
0532 }
0533 ],
0534 "metadata": {
0535 "kernelspec": {
0536 "display_name": "Python 3 (ipykernel)",
0537 "language": "python",
0538 "name": "python3"
0539 },
0540 "language_info": {
0541 "codemirror_mode": {
0542 "name": "ipython",
0543 "version": 3
0544 },
0545 "file_extension": ".py",
0546 "mimetype": "text/x-python",
0547 "name": "python",
0548 "nbconvert_exporter": "python",
0549 "pygments_lexer": "ipython3",
0550 "version": "3.12.5"
0551 }
0552 },
0553 "nbformat": 4,
0554 "nbformat_minor": 5
0555 }
