commit df661e255c61cce470e307e1dfa5ef414f2dce97 Author: xiangyu.xie Date: Wed Oct 8 08:36:06 2025 +0200 Add MC generation codes diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..95f7af4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.npy +*.pyc \ No newline at end of file diff --git a/GenerateFromGgdPar.ipynb b/GenerateFromGgdPar.ipynb new file mode 100644 index 0000000..9acaee9 --- /dev/null +++ b/GenerateFromGgdPar.ipynb @@ -0,0 +1,179 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Welcome to JupyROOT 6.28/10\n" + ] + } + ], + "source": [ + "import os\n", + "import numpy as np\n", + "from ROOT import *\n", + "from multiprocessing import Pool\n", + "from array import array\n", + "from scipy.stats import gennorm\n", + "import imp" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Processing Moench040 150V 15keV simu2 with 128 bins\n", + "Alpha fitting chi2/NDF = 1.9768728049976632\n", + "Beta fitting chi2/NDF = 1.3047556961613742\n", + " FCN=245.132 FROM MIGRAD STATUS=CONVERGED 96 CALLS 97 TOTAL\n", + " EDM=3.53827e-13 STRATEGY= 1 ERROR MATRIX ACCURATE \n", + " EXT PARAMETER STEP FIRST \n", + " NO. NAME VALUE ERROR SIZE DERIVATIVE \n", + " 1 p0 5.20982e-01 1.41176e-02 3.53269e-05 -1.23898e-06\n", + " 2 p1 2.73449e+00 1.50597e-02 1.01242e-05 -3.50997e-04\n", + " 3 p2 -9.48420e-02 3.32013e-03 1.69238e-06 -3.00452e-03\n", + " 4 p3 4.99907e-04 2.88035e-05 3.33069e-08 -1.76453e-01\n", + " FCN=160.485 FROM MIGRAD STATUS=CONVERGED 412 CALLS 413 TOTAL\n", + " EDM=3.93066e-07 STRATEGY= 1 ERROR MATRIX ACCURATE \n", + " EXT PARAMETER STEP FIRST \n", + " NO. NAME VALUE ERROR SIZE DERIVATIVE \n", + " 1 p0 8.57590e-01 1.20963e-01 5.07321e-05 -4.75170e-02\n", + " 2 p1 -2.47722e+00 1.29467e+00 1.23369e-03 -3.44551e-03\n", + " 3 p2 -3.91474e-01 3.56163e-02 1.94337e-05 -5.53030e-02\n", + " 4 p3 -3.36737e-01 1.31473e-01 3.89495e-04 -1.89786e-03\n", + " 5 p4 -2.81985e+00 1.71232e+00 3.71391e-03 7.06359e-04\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Info in : created default TCanvas with name c1\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABjwAAAMECAIAAACE6hs7AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nOzd3XKySBcG0OaruaN4rrkyzZWJ58k18R302NMCIv7SwFo1NZXXoEFDpH3YvbtqmiYAAAAAQEn+N/UOAAAAAECb0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4gitAAAAACiO0AoAAACA4vwz9Q6UqKqqqXcBAGagaZqpd2F16rqOX+x2uyn34x5GVgAwhpFV15Sh1eFw+Pn52W63h8OhtIGXYwUAhkkiPmy3251Op/yW4/FY2gjqGiMrABhmZNXrE9MDq6o6HA7d2w+HQxzBfH9/p8uGAABrttvtqrN0Y1VVrcQqGEEBAEv3oZ5WPz8/cezVTa+6IzAAgHVqlVPF3Cr+f7vdHo/H4/HYNM3xeNxutyGE7+/vqXYVAODdPhFapXFVOKdXu90uTgmMg7AZFbcDALxJXden0ymFU3EEFQdLTdPUdb3b7eKQabfb1XV9PB5D1uXqAXlVVxye3bW3+d3vui8AwBjVJ1sMdHsxbLfbAsvaq+qjLwsAzJHT5cvFTKqV/lRVtd/vr0VCVVU9dvGvruveKq2RY7PuoG7gvg4VALjJ6bLXh6YHHg6H2Ithu93u9/v9fh9rr06nkxorAIDoY+OimFhtt9vmLNZtnU6nmzVTh8MhJlZxrmJ+3wIvRgIA8/WJ0CquEhiHNXVdHw6Hw+FQ13XTNPv9Xm4FAHDNdru9NlKKCdED46gUS+UZ0263i9nTz8/P8N3jBnmF12632+/3+SMDADzvnw/8jJ+fn2tFbrGzlR6iAABxUNSa7jdQu/T9/R2jonvFx+zeN/3c2K9q4L6hE5aJqwCAl3t7pVVd16kLe6/dbldmZysAgE+Kuc/393c1TsgWaE7N1MePqR4rdY872RuWxWr6Bx4TAKDXJyqtun06793g8+JAsJfWaADAmzRN09vjfKSRd7wWbF2rorr2I+q6TlMUdXsAAF7u7aFVWpg5zgTsbhCvyJU20JFMAQCTeLj8PEVID4sdG8ZMOYwpVQqwul2uAACe94klFWMj9vh1aiYaB1VxrFPa9EArTQLATU6XS5JGazdHZXk1ej6uS+sJdnOrgQL2YQ4wANbDyKrXh16UPLdqKS2xCo4VABjB6XIZ6rpOS+Ls9/ubTalSAtXKp1LhVfeocKgAwE1Ol73e3og9OhwOx+Ox1ZF9u90ej8fSEisAgJVIizjHUdn4NuqpxioxogMAXu4Tjdiju5azAQDgra7VTA3bbren02lg+7qudbYCAF7ic6EVAACFiIlVgV0aAACSD00PBACgEHEa4GOJVbzvwB2VWQEAryK0AgBYl7g8zph0qa7rw+GQ97qK9zqdTq3cKt7eamAKAPAM3el7aNoPADc5Xc5X6mZ1TfrNpgWg8991ujG1Y0+LRPceEg4VALjJ6bKXnlYAANwhFl79/PycTqfT6ZRuN9QGAF5LktdDwAkANzldkqYN7na7gcmGDhUAuMnpspcXpYdjBQC6/s5zyr6aJjhdMppDBQBucrrsZXogAHDD360WSAAA8HJCKwDgqlZc9eUCIAAAnyK06jewqo6CPQBWIk+sxFUAAHyY0KqfZAqANRNXAQAwuf9NvQMAQFkkVgAAlEClFQDwL3EVAADlEFoBABquAwBQHKEVAKyauAoAgDIJrQBgpcRVAACUbAmN2He7XXW22+0Oh0PvZofDId+sruuP7iUAFOOvqlrtqyRWAACUpmrmPEit6/r7+7t7+3a7bWVSu93udDq1Ntvv970JV1XN+2UBgGteWF3ldMlI1eVR1+IoAoBgZHXFvCutYmK13W6bs+PxGEI4nU55GlXXdUys9vt93Gy/34cQfn5+1FsBsB6qq5hKc93UuwYAlGvGSd7hcPj5+QmdC3Sp/CrdHq/vteqqYu1VtyYrCDgBWJxWXPWSx3S6ZCSHCgDc5HTZa8YvSkydeqf4xZTqeDzudrv0z5vZVn73+b4sAJB7X7d1p0tGcqgAwE1Ol73mPT0whBBjqQHXJgCmO5ohCMAi6bYOAMCs/TP1DjzuWtiUbo+xVPzndrv9yE4BwPTeV10FAAAfM/tKq6446S+2Wh8WkyyVVgAshuoqAAAWY1Gh1eFwiO2rttttanQVM6mbswhbqoe89ukAwHjiKgAAFmbG0wNzqat66KwSmDa46wH1PwNgLkwGBABgkZZQabXb7WJitd1uj8djK7EaqLE6nU7DGwBA4VRXAQCwVLOvtEqT8o7Ho/gJgPVoxVUT7gkAALxDNet5cKmD1cDsvzRzsPtM4917b5/1ywLA4pWQWDldMpJDBQBucrrsNePpgXEa4HBiFbLZf63NLBoIwBzlDdfNBwQAYMFmHFr9/PyEcR2pttttOIdcSSy/2u/379g3AHi57vqAE+4MAAC824zLz6rLxZK68qeWJhIeDoe6ruu6jl3Ye5++qjwAilLm+oBOl4zkUAGAm5wue834RbkrtEqdra5t0Hrk+b4sACxJmXFV5HTJSA4VALjJ6bLXul6UWGMVQtjtdgPzCh0rAJSg8MmATpeM5FABgJucLnt5UXo4VgCYVuFxVeR0yUg3D5VW+bzjCoAVMrLq9c/UOwAAXJhFYgUvFMfoBusAQIvQCgBKIa4CAIBEaAUA0yu54ToAAExCaAUAk2m18vkVV7FQA4s+mxIIAFyjd0APLRUA+IBUXbWZZ1zldMlIIw8VRxQAa+Y82EulVT/XAwF4n7/OWWZ2iRUAALyb0KqfZAqAd+jvXXX9SgkAAKyW0AoAPmFMq/VWna8rKAAArJnQCgDeaziuikFV/H9MqXQ0AACAILQCgPcZU10Vpbjq7fsEAAAzIbQCgLfIE6trcVUsqkpZlegKAAASoRUAvMzVvCmb/dfauDs9EAAACEIrAHih3xBCCJvzF1+dHlVarUPkbwEAuEloBQAvMLJ9VT4B0Kd01syyAwDATUIrAHhQq1Tkv+qqR5tSpQeUagEAgNAKAB7UNM1fVXUnA4ZOntXtddUbSEmpAAAgEVoBwCO6iwOm8qjWjKfuBCjrAwIAwE1CKwAYq3c+YDCbDwAA3kBoBQBj9c4HDH2LA4qxAADgSUIrABiluz7gzU5V3QxLqgUAACMJrfoNdBvxAQNgbf46J4VN+K/A6q5ThpMIAACMJLTq50MFAKGvuiqEi/mAAADAmwitAOA/11qtf3W6VgUdrAAA4J2EVgDwn1ar9ZDFVb0bf2avAABghYRWAKzXtbqq6N9W65fbNH3914Fh/moAgAcIrQBYr9RJPRZY5d+KBVb5BvFTd/5/lVYw0rU/FmEWADBAaAUA/zVc/8rCqfxjdqX5OgAAfJbQCoAlqzqT+1ob5AVWm/DvyoB5aZWsCgAAJiG0AmDJ8vl94XoTq00ITSedyu8iugIAgA8TWgGwIimH+u27PWSpVsqqdF4HAIBJCK0AWJHebuv5+oAqqgAAoBBCKwBm7GbLqqQ/rhr9+CYJAgDAhwmtAJixVsuqXlU3rrqMosL1TEpKBQAAUxFaAbBkf1X1G5cF7KuukkkBAECxhi5Nr9bwFXsASpO/bw83TR9TmcVIXkZGGjhUWn+8jigAVst5sJdKq34Dn3kcRgATutnEqmmav6rahJDWB8wLrPSoAgCAuRBa9fMxBqBMrVKpVgj1e7mx+YAAADBfQisAZikvlaqq6rezwSYEARXMxfDEXgBgncyZ7GEqKUDhupVW0W8IX01j6t9nOF0y0pieVv5sAVg5I6teKq0AmJk8qMoLrNJkwOacWwEAAPMltAJglv6yWKrbuwqYkXhtWdYMALT8b+odAIA7dAunW4lVVVXxo2/6AgAAmCNzJnuYSgowoVbS1HQyqYHv8klOl4x0s6dV/v8P7xsAFMJ5sNdyKq3quq6qqq7raxscDofqbLfbHQ6Hz+0cAKM1TRNP2OmL6K+q/qrqN4T4X+u7wEwZowMA1yynp9VwCNW6Mn86nU6nU13XAyEXAIX4u3wP18EKAADWYCGVVofD4XQ6XfvubreLXxyPx3hl/ng8hhBOp5N6K4CS/VVVVVWbENJ/EisAAFiJeddjHw6Hn5+f/Jbj8ZgiqiSWWbW+Vdf19/d36OuHokwdYHJVVf1m/9yEEHSwKozTJSMNL4mQjiJHFABr5jzYaznTA29qhVndbAuAQuTzATfxM62zOMxZ/vebD8ot8QkADJj39MDD4dCc3dy41b5KNyuAMuWJ1VfTxDXFJtwfAABgEvMOrUbabrchhDgZMEpzA/f7/WS7BcCluD7g5jwZcHOuwlBjBQAAK7SK6YF1Xe92u9Pp1LpWv9/vNWIHKEFeXfWr2zoAALCS0CqEEEOr8ds/NhVFLQDAGCO7MgMAAGu2iumBu90uLjJ4PB5jA6zj8RhC+Pn5udaOvXnIJ58UwHyl98zf7L/0rdbGVVXFkCt9ASyDv2gAYNjyK60Oh0Ossco/CO12u9jZ93Q6xcmDk+0fwPr8XX5S/Xcy4JWPry4JAADAOi2/0iquEtjbcD02aLeMIMDHxFbr6Z9fTfPVNKmESj4FAAAkyw+tBsQCK6EVwAd046pwnvFnhjUQmTAIAOSWH1oNJFPxRnMDAd6tt7oqZM2tAFqqS1PvDgAwgbWEVrF3VX57Xdex19XhcJhgtwAWrffTZoyrpt0xYC5SqC3dBoDVWn4j9t1ut91uT6fT9/f3drtNhVcxsertdQXAk5rz7L+0LGAeV6UYq1U9obMVAACQLL/SKoRQ13UMp06n08/Pz8/PT0qslFkBvEPewaq3wKq57uM7C+u12+1SReRutxs5LjocDtUVRlYAwAstp9Jq+HPO4XA4HA51Xac+VlpZAbxQq2aqt8AqbSmZgsnVdf39/Z3fcjqdYjuFz6xRE98KNKsCAAb45NDDByqAB/xV1SaLqzYhhL4rCt5jF8OvctZiWrTdblNElWKsm6Xou93udDodj8eRlwBbh0qKq9L/811qHVQOMwBWwimvlxelh2MF4C5pJmAMrVJclT6CtoopvMcug9PlfB0Oh5+fn9D5Y0y51fBv9t72c9dCq/h2IbQCgOCUd8UqeloB8CZ576poc7nUVx5dpVs+vJNAS6yu6i5HkyqnPjNDEABg2HJ6WgHwWsPlUd3v5uFUa33AbnQFTE5/TwCgcCqtAOiXCqa6i/r9VdXvuX1V/G6rnrm7IKC4CspR13XTNN3QKhVYDeRZ+TYvWTfwTy92AOAKoVW/aws5W+MGWK38nTB2rfq67EST/u+tEmYqNWIf2CaFVqfTKYSw3W7jP39+fu762/+9vQkAsHZCq37NdVPvGsA0mqaJHzJjjVWeWHmrhLk7HA5pPcGRNVPH47Fpmli0dTwe0+P0bty6/teKt/Kw25VCACDRnb6Hpv0AUfrQmGoibi71Nfxo3l2XxC90GdKKgSGE/X4/JrGq67o7f/DaioThyuqBcVbgJgvBrR4IwJo55fXyovRwrABE8TNkPotnc9lVfeQb5nBPd2bK6XIBdrtdmuV3OByebM1+Lcjuhlbh/MYitAKAyCmvl9UDAdbr5vqAvyFszv/8dz5gVT1wQnUChgKld4Dj8fjClQR767CG/VXVl3cJAKBDktdDwAksxpgSp/SmN9xBJm3Tu7G3zRVyupy11MEqNVYfKW7fG0s9XGkVQvjKSjgH7gsAS+WU10ulFcCStZKmkRvn8wGvlT84p8J8xcZVDyRW4bzCYLc464GHSjYhhKwRu7cXACASWgHQ07sqdOKqfG2vILSCOYsd08dM4qvrOqZRqUH7drs9nU7f39+tN4EYZu33+wf2J29r5b0FAEiMDHoYMAELM/C2lhKov3MmtRFIMY7T5XwNzwUO2ZtA75qA6e4poorbhFtzkNM/8wVJ0/TA7pa9twDAIjnl9VJpBbBw8ePl8FkwJVZfTRNufZoFVu54PMa6qpRVhYcmG+ZvOHqxAwBdkrweAk5gMeIb2kB35PyfdzXAAofKyqWZg3Ga4cBkw2uVVrH5eqr0TP8cuC8ALJVTXi8vSg/HCjBfAzlU683tL9vy6/LzpPdAxnCoMJLQCgBucsrrZXogwKK0UqpW9/T882GSEiut1gEAgHIIrQCWqfdaTR5jhU4sJaUCAADKIbQCWKZWAhWrq+KUnE3fBgDT+rMKBABw6X9T70Chquum3jWAG7rvVPlHwc1ndwbgJusGAgC9VFr1U4AALEO323qjyyNQgE0I4dxrL5zrQAEAciqtABYlD6SurQ8IMLnm0tS7AwCUSKUVwPy0JgB2V4gPWfGCuAoAAJgjoRXA/MSUKhVVtTKsfJbNV1+elb5Q3QAUxSRBACAntAKYvRRdDcRV+ZYApflqGqsHAgAtQiuA2Wt90jMfEAAAWAChFUARhttUXSOuAmah6lRRpVtUgAIA1witAIrQalN1U09cZWYNUKr8nU1PPQBgpP9NvQMA3C1PrL6aRoEVMBc3G1fpbAUAJCqtAOYkfpzbnP+5CSGMLs4CKERv1D6mF/tjM6kBgJlSaQVQivhhLP0/FzdIH+d+Q2gyaZt8Y4DliW96+RcAwIKptAIoSwyh0hep11Vrm4F/AgAALIDQCmAyvVVReWIVb/mrqt8QwnlWoIgKAABYA9MDASaTT29pLoVzOJWmBOq2DqyEXuwAQKTSCqBErc9sseF6UGYFLNqYXuwAwHoIrfoNdDL2iRF4ufyNpfWBLRZYxW9rsg4AAKyH0KqfZAp4h24/9fyWPK7aXN7FmxKwAPdWUUnqAWDl9LQC+JxWv6q0OGC4/Cz31bRVVRU/vKUvABbMJEEAIKi0AphEjKvyKCqE8Hu927pKK2AlNiH8Tr0PAEAhhFYAE0j1Vq0Cq+n2CKBQeVEqALAqQiuAV+p2rbq2gbgKAABggNAK4JW6zap6SawAAACGLacRe13XVVXVdT2wwW63q84Oh8Pndg4ghJCVWW1C2ITw1TQSK4BeerEDAMvpEbDb7U6n0/F43O123e/Wdf39/d29vffpa50APKn3baSqqry7sLiKuXO6ZKR0qMQcajO4uETcOG0ZtbZ37AGwPM5uvRZSaXU4HE6n08AGMbHa7/dp/fh4e2/CBfC8tCxgkr6lwApYIZVTAMC95h1aHQ6H+FHw5+dnYLOYTG2323xK4PF4DCEMR10AD4j5VNM08YtYXfV7/kJcBazZ5vYmAAD/WkUj9phMtZpY7Xa7/X4/zQ4By5UqqvLEKhJXAdzld3TIdXPlVgBgjpYzZzIOVro9rVI3q/HP1FRSoOWuj0Ot3lXhVgMXmCmnS0bqtqm62dMqZA2wuttfO/YckwDMl7NYr+VXWsX1BLfbbfxnqreyeiAwXjx/jDyR9FRX6eQCrFsrzc+vBBigAwDXLD+0Sg6HQ9766ufnZ7/fi66AF/q77LY+4Z4AFCVVWqV/TrgzAMBcLD+0ipVWp9PpdDptt9vdbrfb7eJqgzHD6s2tqofKIozAYM0kVgADVJ4CAPdafmiV5HVVdV3HXlc/Pz8xxmptLH4CWmKWHWcItnLtfM5Lq31Vqy+79xaAkTaX764AwAr9b+odeLsUSHVXD4xfxFIsgDFSYhXjp6ZpBhKruEHuo/sKAAAwZ8sPraLUiL17o9AKuKlbMBW/SFMCv5rGlEAAAIAXWn5oFSuqTqfT8AYALVUmXNZM/V0ug7U5bxwe7YgHAABAy/J7Wg1kUjHJElrBmrUypnwGX/y61Yvq73L7zfUtAWh5rM2fzlYAsFrLr7QK52mArXAqtbgSWsGapVZTN3tO/VVVd33AeK+8GuvN+wswY/n7be8bZipZjV+8Y9p1denljw8AvNDyK61CCHVdV1V1Op2qqtrv9/GWWGZ1PB6n3jugdL+XBVZfsaFVNhkwr7QC4GH5YhdvelNN79guMwBA+VYRWoUQjsfj9/d3COHn5ye/UZkVkHSnCrYmA6Zr/j7qAAAAvNu6rjLVdR3XCtztdgNxlYtvsCopq8ovv1fnVuubEH7HTVHx1sHaOOYZKT9Ubh42eaVV/m6cvw9fe5Dxx6SjF4DSODf18qL0cKzAqnSnoqS4KjfwtjDQzR0WzOmSkYRWAHCTc1MvL0oPxwqsRzdvSh+Qwnt6AMNiOF0yUjpUxhwz13paNUIrABbNuanXWnpaAUQDn4Ly9QHFVQCTi2/Rf5a5AIC1+t/UOwDwUfk6663EKn0tsQIAAJicSitg7VJclZpYpckpE+0RAG1/VeWKAgCsjUqrftV1U+8a8BrdWSfNpel2DYD/yKoAYLVUWvXzeRUWLKbP5gMCAACUTGgFLFy3QLK1KJXECuDd0rvuyKWRbm5miSUAWAPn+x6GQbA8Y9ZQB+7idMlI8VAZ2S4wvUv3rpURLzNce6jxx6SjF4DSODf18qL0cKzA8sS/a1MC4YWcLhlpuCVo6yjqDa2qqvoNIVyGVt3HEVoBMF/OTb1MDwQWqPfzTOtCPQAfM77SqlsSO7BlMMQHgEWzeiCwQGn5v/j/3+xbEiuAWbOUMwCsh0orYMliddUmhPR/kRXAfP1JrABgTVRaAXNVXZdv9hvC77n2yhQSgJlSJwsAK6TSCpirbjeT9LWG6wAAAHMntAIWJa0wFUmsAAAAZkpoBcxePh8wJlaby2+ZFQiwGL/ZmzwAsGx6WgGz0W1clTKppmnyAqvm0kT7C8Ar5cWzTdNYRhAAFk+lFTAbeROrmFvFL/6qKl11t0QgwCKJqABghVRaAfOTGq7HzzAxqLJEIMCCxff23/M0cABgDVRaAUVrXVqPH1rSR5dEw3WABTMZEADWSaVVv+q6qXcN1iVVTuWlVSGEv+yPUWIFsBLe8AFgVVRa9TO9CEqTTwlMiZVPLwCL9MBIzJVFAFgeoRUwD/nckM35lgn3BwAAgLcSWgHz8FdVqYmVAiuAlftTVwUAKyC0AorWqq76lVgBrNtm6h0AAD5GaAVMr3eJwHBZXRXEVQDccu2EAgDMkdAKmF7qsJ5/urA+IAAt8UwRk6m812Eu3i6uAoAFEFoBxcnjqo3r5AD0yU8WeZgFACzG/6beAYALEiuApfKWDgDcRaUVUIR4eTwlVv/OB3TNHGAR4pv8tdl8H98dAGAehFbA9Kqq+j0vCLU53zTh/gDwWvpMAQAPMD0QmEB1KVwmVulTjY83AAz4vb0JADBjKq2AN7q29HhaLjD/vJGWCEzLQn1kHwGYn6+m+VOTCwBLJ7QC3iiFU90EKv+w8XX5XStAATCeiYcAsFSmBwIf0p0SGH35pAHAaKmn+2bqPQEA3k2lVb+BKg+X8uAB6TL4tSmBADBGGok9OT1QfRYAlE9o1c8gBl4lNajK46pN57ut7Xu/BQCJtlYAsHhCK+CN4nXsqqr+ssTqazCIklIBAAAQ9LQC3ie1HQkhbLLqKgB4oXzWuaU8AGBJVFoBr9T6qJB/kEi5lXl/ALzEV9MEERUALJfQCnil1G09XCZWw1MCAWCYEioAWCGhFfAWeQerKfcDAACAedLTCnhWdSlfy0liBcALXSvbtYwgACzSckKruq6rqqrreszGu91u/MbAsKZp4qeI38sCK4kVAAUyzRAA5mI5odXhcBi5ZV3Xp9PpnfsC65V3W/epAKBw8TJetNvtxo+mWrznAwDvsJCeVofDYXwO9f39/dadgbXJJ2Xotw4wC3Vdt0ZEp9PpdDrVdX1vKfput3vhjgEAJPOutDocDvHK3s/Pz8i7uAwIr6WDFcAcxcRqu902Z8fjMYRwOp3uqre668LhW2lrBQDLM+/Q6l5xELbf76feEZi3vO36xpRAgLlJsVReVLXb7WJuNf5a4L0bv8ljl0yapmktJPLyHQMAnjTv0OpwOKTLgzc3rus6jqse7tcAa9Ya1seG67HzenNp2v0E4KaYVXUv46WJfiNnCMYzQoy6Ziqetpy/AKBMC+lpNUYsgzcigcfEv52qqn6zG80HBJivJ3tRxbvv93s9rQCAN1lLaJXGVSO3f6xEXCIGAJTvWiFVuv1mDpXWYi6kgP2rafS0AoDlWUVoFVuEbrfb8eMq8RO09DRc9/EAYFliWfqYi3z3FrDHy4HjLwrGLe99/Lt+BABQvnn3tBoptrK6d/1moNVwPfo6d64NOq8DLEVckTmEMOYiX6zDuquVVd466qZwf5Oppml+s06L4+8IAJRs+aHVAlqEwlTiZ4Bw/hiQOljpvA6wGHVdV1UVr/Dt9/ubF/lSAXtpray0WQSA5akW84EzhVOtIdSYMpDWi1BVy3lZ4ISOBWoAACAASURBVC6tv5eYWG1C+A1hc7mlvxHA6XIBdrtdbE0VC6zGtLLqnRg4PJsvHirjD5hrW/bent8YZ7J/3fpZaW/v3TEAeB/no16r6GkFjNS7RGAI4cvbJ8DipAsV3Wt+4+/bvXG/30/bnV1HdgBYjOVPD7zWKyGEcDweTW6CYWZbACxS6mDVNE1pE/0e5pwFAAuj0grWrndKYDD0B1iuWAm13W7vXaZmt9tdm7gXzBwHAF5NaAVr1zsl8Ot8Y/ynTyMASxLbro8psKrrOgZb0874AwDWafnTA4HxNueG61VVxUaAVgkEWKqfn5/qirRNXdc/Pz8x5JoFZcIAsCTLqbS66xO1j9+Q5A1r/WkAsB5x6cB3/5TWj3CqBYDxLKnYw0qTrEdKrDaG0cCdnC4ZKR4q4w+Ya1v23t698a+qNiH8Xqm66k2pPnAk+3sBYJgzRa/lVFoBd8njqkjjKgCm9ZJeil9NE67XT7Xisw9UWgEADxNawRr9maoAQHmcjwCAnEbssAq9HXa/mkbDWgCW6k8VFQDMnEorWIV47bqqqt8QQpoSaDQPAABAqYRW/QYaHChcZ47yq81f54NYEysAAACKZXpgv+a6qXcN7pYnVpvLTNYhDcCCjZ8h2DuPHgCYltAKFq61SmDMXo3IAVi2ze1NLqRrky5SAkA5TA+ExeqNpfIpgXIrABbvr6oeXnWkstguAExKpRUs019V/Ybwe/5na5ZrqrRScgXAUj2/Qq7yKwCYlkorWKA0JfCraUJVdcfZRt4AAAAUTqUVLE2eWKmiAqAEeXnv51ueD7RjdxUHAEqm0gqW46+qLvrOVlUjtwKgAPF89PmE6KtpWonVJLsBADxGpRXMw8B16XRjTKx+L78VNK4CYInu6s84UGwFABRLpRXMQ1rvr3t9OKZUmxB+Q9icF/l2GRmAZbt5pmuv/ffOnQEA3kFoBfOWXzreZNnWdHsEANNLExKrqvq9vTkAUCLTA2HGUmKlugqAkpUwXd0MQQCYHZVWMDPxunG3wCp9HkibBTEWAGWY9nzUbcc+oD2p0JkUAKaj0grmIb86XV0mViG7gh1CaDKf308AKNaY6CqdQJ1JAWByKq1gHlJvjr+qSt3W07cm2y0AmIO7iq0AgEKotOpXXTf1rrFqccCdN5SVWAHAkwzwAKBMKq36CQIoTb740VfTNOfmVlPuEwB8xAs7NlpJEABmRKUVzMDfZWI15a4AwMc1lx57ECdQAJgdoRWULu/BYcANAE/S3AoA5sL0QCjXtVW3XzhLAgDWQzt2AJgXlVZQrt+s9UaeTL1klgQArNZf57KQXuwAUCCVVlCW1qA5hlaSKQB4Ru/pFQAonNAKptcdSW/O4+nNJDsEAPOX10/Fyz/xnxIrAJgLoRVML42k82H0Jk4D7ERaAMDz/qrK8iYAUDg9raBQ8fpwTKzSFwDAgDGnTlkVAMyFSiuYwEBnjX9H0pYFBID7tU6dw5d88mIrvdgBoEAqrWACadW/pmlaiZXqKgB4N8VWADALKq1gSnHJ7dRtvVJgBcASDVyGmfas9+f6EAAUTKUVTCMfJTeXJtwrAHiH5rqpdmnyYqvq0rQ7AwBlElrBh1wbm04+aAYAPi/vFeCSFQD0Elr1q66beteYqzQkjU2sNuf/O6gAoASSIwAojdCqX4FF7CxAmhL4e3mMTbtXALBOip0BoHBCK/iQP1MCAaBIn2nHrrYaAO4ltIJPkFgBQIGclAGgZEIreLuYWG3Ofaw0RwOAzxg/Bz9dXnKOBoBy/DP1DsDC5TVW2lcBwLzIsABgQkIreJf2lECjXgAo2Gc6W31eK3dzCQ2AGRFawYu1hoa/+mUAQDF6I5uvpulNrKqqaprm+WKrxx7hVWFTvGN8Lo89AgBMRU8reIEqE0L4DSH9f6OJFQDMx0tSqpdomibGTOkLAFib5YRWdV1XVVXX9bUNDodDniwcDofP7RxLlw8rf883fjVNPso03ASAYimLBoACLSe0Gg6hqqr6+fnJb/n5+RkOueAB+eSCSgt2AJibfBnBdPoupPYKANZmIaHV4XA4nU7Xvrvb7eIXx+OxuSx++f7+/sDusQbttuvn+QWGuQAwC6nY6nd4OwDgU+YdWqUZf60qqpaYZx2Px5RehazyRbEVz6uqahNCCGFzbmIVzv9XYwUA03rsApIzOABMbt6h1Rgpk8oTq2i73QahFU/7q6rfEH7Pna1yIQuwlFwBQFHS2dlpGgDKNO/Q6nA4tKb7XRPzKXhS1SdWV/U2cO1mWABAIXpP0zqyA0A5/pl6B95ut9tdCwvitMFuBRZck4qn4hd/VbUJ4dcAFwAKFrtM3nv16K+qnN8BYFrLD62uSVlVb2il9wHD8rbrmxD84gFgMb6aJnxwtmBr2Gk8CQDJvKcHPqau66qqYpnVtWFBtznRGJ99HpTCrx4ACheLre69199Hoqs0jDSeBICW1VVa7Xa7GFcFWQMPSePXf6cM6NsKALcMBEbrGY89NksRANZsRZVWh8MhFVjt93sjBh4Q266HEDaWHAKA0eZbqP6ZYisAoNdaKq1SgdV+vz8cDlPvDrP0V1W/IQRt1wFghu6dIfjVNBIrAJjWKkKrw+Ew3MEKhuVjVokVAAAAfMAqpgf+/PwEiRWPai0UCACsRLpS9fvo0tIAwDOWX2kVJwNut9upd4RZytuua54KAGvz1TRp0ZU8tzIkAIAPWH6lVV3XIYTT6VRdocUV18TEanNuux5cZQWAtfo9B1Wz6B8PAMuw/Eqr2M0K7pVqrOLAVJkVAKzcY33Zq75Cram09sHYBoDCLSe0unbSdTJmjNYY7vf8hbbrAMAzKwk2WYeByXMrl+IAmJflhFbwjHwMl/exil9MPsQEAMrxTL2SxlgAMJ7QCi6kPlYhhHA5JI2jTONLAFihzbkQ+y8rm3pgVNBY2gUARlt+I3Z4QOqxmrdc1XgVANYslWA/PFUQALiLSisIIavVb9VYuRYKAKyHTu0AFEVo1W+gh5GT9/L8VVVv53WJFQDQa6nFVjq1A1AU0wP7NddNvWu8WD7otFYgADDs5aOFV633Ul16yWMCwLRUWrFqf9mswOayxir/QlgJAMvWOvUP3BhC+MrWGh7zyNcGEq+NlhRJAbA8QitWKh9rbjqxlNEeAKxK69TfvWrVGzD9VdVdhVdx6cBH9xEAVkdoxRrFxGqT3eKyJAAwoJU3fTVNXLbl3twKABhPaMW6dNfEEVcBAABAgTRiZ11+Q4gLBca2+kr0AYAnPbmSoAEJAFwjtGJF/i5bq6qxAgCekSYGPplbAQC9hFasRRpNxvFlWmEnnAOsCfcNAJipGTW0MtoBYHb0tGL52kO08z+VWQHAZwzEJe87Hacf2l0K8B3u7cgeK76XGiR1u4hOtScAzJpKK5Yv9bHKqa4CgI9prvvYD33fDzJJsCu95u9+8QFYNqEVC5cPHz85UAYA1qPwSYIu1AEwU0IrFuuvqlJitZl2VwCAmUi12I8VZT9WbDV8IU3kBMBqWT2th0XlFmCS3hkAq+J0yUjlHCpP7kmrN1b30VqrvvRulh4k9bRKi8P07tu9Dbl6H+faD73rQR5Qzq8eoHzeM3uptGKB/qoq72NlSiAA8AGaWwHAawmt+lXXTb1rXOj+dloXOaVUAMDHvLa5lZEnACsntOo3yRo3PKC1Nk2eWBnnAQAvdNc4cLjYyhwQABhDaMVy9CZW6uMAgE8qeZJg46oeALPyz9Q7AK/xV1X/LRF4Z8tSAIAX+jpXfxeYWwHAjKi0YvbSBcPfLKiSWAEAE0r1Vr+XtxuiAMB4ptP30GVgRnpL3P36AD7A6ZKRSjhUWgOGh/cnfy5jnldrfZh8Z+J904P0Plra7fE7PPA46QfdfMBX/cpK+NUDzIX3zF6mBzJv6erlVzaCnGpnAIAyTfUx4CtbJQYAuJfpgcxY6+pl6rmu+ToAjLHb7aqz3W53OBxG3rGu69Z967p+447O35jo6jOjF5fxAZgR5Wc9VOXNQm+9PQAf43Q5a3Vdf39/d2/fbrc346fD4fDz89O9fb/f98ZeSzpUemf23bxLqzC8d45h76y9azHWwM/t7lW85a6JjaYHAnye98xeKq2Yn3hRdxNC/A8AuFdMrLbbbXN2PB5DCKfT6Wa9VUys9vt9uu9+v0+3L14MgNI/q0u9d0kX2B6bKhg/w6T/j0zKhncJAGZBaMX8/J5bWY0ctwEAuRRL5UVVu90u5lbD2VO873a7zbOt9PX4CYazludWaTQyPCx5Mrd6YA/zqAsAZkpoxczoZgoAT4pZVSyPyu12u3yDgft2w6n4aOvpbPVAGKShAQDcy+qBzElKrMwKBIAnpYjq3ntFrdtjXPXYYy5Yq0FJWkzwr6pkWABwk0ZfPfQ/K9O1xMovC2ASTpfLk7qz3/ubTXc8Ho/d3GpJh8q1dubd1lEDzdpja/beK3C9jdhTG/Uxr2R3DzViB5gF75m9VFoxD62xoD9mAHi5GDx1pw1eczgc6ro+nU7xn/v9fs2VVnctKRhC+L2cMKhjOgB06WlF0XoXvpFYAcBrHQ6HeLZtdVgflidWYXBuYPWQp57SRzy5kzp1AsAwoVW/WY+flqRpmrhQ4G92y2R7AwCLU9d1VVVxxcD9fn9XJ/W6ro/H4/F43G63IYTv7+9rd28e8oKnV6oPLyYIADNlzmQPU0nLEUdysenDyFYOAHyG9+QF2O12sVQqFlg9M7nvcDjE5Ku3hdNiDpVuH6vWrMDWk73W0yremBKrr6bJHy3fMuhpBbAO3jN7eVF6OFYKMbBWoF8QwOScLucu1Y/3dk9/+AGFVuNDq9AZ7XTvey14GrmHQiuAufCe2UsjdoozMAfT3zAAvETqYHXXfMBwPZniMV9N8/IZgj72ALAYelpRnLyP1Xp6WwDAx8RW6w8kVkn3jg8/FF/vHOHEkDH9P10aHLhGqIUrAOUQWlGcvL/DtHsCAIsUm0+NmRJY1/XhcMjXE4w917srDKYg7FU7OWv5Aj7jM6Dfzi2tyGnMz+3e6BIgAPOleLiHmuoJ/VVVq4OV3wVAmZwu5+tm/JF+s73t1dPd9/t9/CJuE66ctZd0qNzsGHWt6dW1B8lvjKHV15XOU3lT9vyOvb+atDPXdqD39jHfvcvDv/qB5wiwVEs6Xb7Qciqt4oLNA6Xph8MhXfXa7XaK2AsUa6xaswKn3ikA4MLxeIxf/JyFELbbbbqdZ9xscZUGSAsbKVWZcA6qFvYcAbjXcpK8uGbzteVv0orOuf1+3y1uDwLOieSr53j9AcrndLlydV2nS4C73W5gsuGSDpW3Vlo1WVP2r85aga1Hu/Yg8YvnK60GNhjvgV/9XascAiyJN71eC3lRUu16b2hV1/X393fIUqrh7R0rn6cIHGB2nC4ZaUmHyrtDqzRJMHQu4/WGVteaWOXfFVoBzII3vV7znh6YZvylTgq9WolVvOO1NqJ83l9V/WbNR/2hAgCrlXpadfuydy11qiAARPMOre7SyqfiP7tzBvmwfFZgkFgBAKuXcqub/a0AYNnmHVodDoeb7bqvNVxPswJ1ZJ9KrJLbhLDRxwoAmKdrc/SetIDc6h0vCwBrM+/QaoyYScXJgBSlaZpY9x7/H0c2ackYAIA1SIOf1hBoc/5ivrkVADxp+aHVgJhkqbSaSr4+Tsu0OwYA8DEDQ6CH663eVP/1gEJ2A4CZWn5oFTOpgWWYe1UPecf+L9XfeWJgUF0FAHDF+L7slp0CYHmWH1pF95ZTdUp/RnnPvi9QvFr4O3hpEQBghbpXQ7+yYdIDJVf3/nSjMgDKsfzQaqDGKi4deG8RFk/KZwVOuycAAKWJmVErOfpqmve1Zr933oCpBgB8zD9T7wArEqcE/selPACA+702t4rjsVhjNSaEyrd/4W58WOuZzvq5ACzYWiqtYlHVwAa8W5oSGDsymBIIACzDO4Y01yIhheqvksaiBqUAJVtLaBU6ba0sGvhJ+ZTAzfCmAMDiWMfmhUbOE5xFEONgAGDY8kOrEMJ2uw0hHA6H/Mbv7+8Qwn6/n2SXViVPrNJwxNAEANbDOjavledWL29x9UnKnQAYtoqeVnVdV1V1Op12u93hcKjrOpVZtZIsXq7bdt2IBACYnfzCWxg9nrl32DOyqxQArMQqQqsQwvF4/P7+Pp1OscAqkp68W0ys/p0PaAQGAMzWzXHj5xuTxyuCcbg1st4qT8Tm3kkdgDVY17kq1VjtdruB/utO4S/RqrHyqgIsjDd2RlreodJdeC5fei892WtPfPgFSQ/e3SZ//PxBqqqKC918nffk5g/qfis9he4TGfkIvY9283G6z2XM/r/K8g5OYL68I/XyovRwrLxEDK1SYhWUtgEsi9MlI63hUBkZWnXTrt6Huvbdm6FVCGEzIjW7tm9CK4AJeUfqtYpG7Hze32XxefpCmwYAYKlufth4XxP6fHXmMVMF39E8yzqAALyc0IrX+6uqTQjxv3SFzSJBAMAavDWyGQibvrJR1iRLCuZLAX7+pwOwSEIrXiwOkn6z8YqBCwDABzyZW8U4rLVOIgBMaC2rB/IZrebrQWIFAKxMc9kN/WNSxhSnCv4ObHpdvttCKwAmJ7TqN3CSlsJcc/GiGeUAAHxQa4waLyXmC+PkWkVV7xvfThXhAbAMpgf2a66betcK9VdVv33X9DTjBADW5h1tzu998OGpgnkbB+NbAIoltOKVvs7dNyV9AACTS9HVX1VN0p0dAJ4htOIF4hhok5Waq64CANasnIt2k68qyOdVl6benWl4EWAZhFY86++yG4LqKgCAt+oOtIaHXl9Nk5dcvWu3KEYaja95WO5FgGUQWvGUfLlAXTYBAIr14amCrV7vj90dgJUTWvG4PLGadk8AANbmgVgnH7N11895lXgh8+Hq+ziZy6VQAEII/0y9A8xVK7EytgAAeJUxw6rH1iiMI7eSJwlOPqRsvaqT78/K+XXAygmtuFt7eHT+p9wKAOBJd0VRzfna4b0/5atpYm71V1VK5lsaV2RL4tcBK2d6IHdrmibWk/9e3uhEAgAwF5vzF59pcTXGMw2wrBMHsEgqrbhbPjEwxlQGBwAAczffC5CKccphNh/wWkIr7vNXVem6XH5Oil87LQEAzEje4ir+v3e2oCSCkQSIwGsJrRjLYAUAeMxAUbYRRQm+miZ1Ke3tcvWBJCI+uPr9FiNwYOX0tGKsvJWV8yUAMF5z3dS7xn9SVrWALleLkf5M/L0A6yS0Yqy8ldW0ewIAwDvkw7zfge1YIv3sgQKZHsgo3VZWLvUAACzPyC5Xr6UF0khvnS2oHRVQIJVW/arrpt61ycSrbYr5AQAWLw+qnpwqqFPVC5ktCKyNSqt+zgG5NFLxsgAArMQmhKZp8pIrAPgwoRU3pDHKJgSRFQDAMuQFUL2VUHGDlFvd9cgP7My9dxmYKGfFPYDFEFoxRC03AMA6xXFgGg3GThExwNpcvdPnDDRg0pvpkx6LCAWLwEhCK4bE0UkclziXAACsRyrF+m8QmAVYm76rm2OGi6tKK9bwZB+LCAWLwEgasXNVXgrudAIAMJW0HNA71gW6KzjIG7T/hvB7HiWO7wu+qlbiq3qyzJRlxyicSiv6tWq/29fZAADISmne+mHv3jFYa69eO4T7ulwNUI92mLVW1dsaygOZF5VW9GjVWCUT7hIAQIGaS1Pvzr8+s1d5Z6uHoyuVHVAU5YGURqUVbYYOAAAr8WRToX9nC1ZVUHIFwBuotKLt99x/XbgOAEDLzZxrLumVK7UA5VNpxYW5DDIAAIry1jZS5YslV2kkGb/4WtmLwBpo+fRhXnCEVvxH83UAgMfMdMj02t3eXP5zFtHVkxMkZ+3DccAy0odW23LezQuO0Ip/tZqvT7gnAACUr2ma3k+SvVVXlaFmeT4cB0gfgAfoaUUI2aii8EthAADMwlfTpIHlwy1TYy72hr0DYB5UWhFCXsttWAAAwIvodQXAM4RW/QYu6SyvnPWvquK1rxhdLe8JAgDwEo/N7drEEeZldJUebWG1VMvo3ARQCNMD+zXXTb1rL9aaGLi8JwgA8GFVVVXnRk4LS2Qek5KpvMYqXTddmPSRYZGfHQA+TKUVIYTwtbhrXABAOVZVwx4W+qRuGlk29W9u1am6aj1U64vSKKcC+AyVVqv21xkQuB4IALzcemrYGa9VddVaybrwI0Q5FcBnqLRaLxeIAABmp1WFNOsh3LVeVwAQCa3WS/N1AIDZKWrY9qo26nmrigWEcRTCDBJYAKHVSrUKsCfcEwAA5uJ948b4yH+XVVdfhqmPktcAy7Cunla73a462+12dV1PvUfTSKOBzbT7AQDAE969TOFbL23Gxu2tG7+a5pkVBl/+IlSXXvvgS+WFAl5oLaFVXddVVZ1Op3TL6XT6/v4+HA7T7dTEUmLlHAwAMEfldLXvTaAe1o2ubra7etNoVsN1gGmtJbT6/v4OIWy323RS3+/3IYSfn5+11VvNa2UWAABK8PlrnA9EV+N7bI1J2VzWfcBr40uAVYRWKZbK86nD4bDdbls3Lt5fVW1C2JzLrJyJAQAo2ebyv97oKg9KXtUbvjTmRjzM6waztqLQKkZUud1uF1YWWoUQfs/tCRRYAQAwxrTlM2nsmvpbjam6ekaBMUccuhceXSmzAl5uFaFVlDe0itYWV7311A4AwJPe3Vh97roTBuP/W0HJUoutwqyuOr+jKf5rHxCYhVWEVqnbeiytSjfGGGslvdjzFQMNhgAAClROY/Vcq3ymlax9fj9b0VV4f+HVejzz6UCZFfAO/0y9Ax8Sr7ecTqfWG/HxeMyTrKVqPWunEwAAojRQjF/cHCj2bvC+4eW1sql/c6vLb8Xo6us8k27845c8PBYkjfeOS/I+ScG01hJa3eux97ti38Li/P/YfL3YnQQA4PMeS6k+8HNH+q/q6jyAz6uuBhKHVC/2kt1YgJhwPfyClPApI+5DeiIvmSiaP+YLdhG40yqmB4bz2Wi73R6Px1huvd/vQwjf39+9na2ah3z4SY2Un7aL3UkAANbm5T28utMGB9YgujmGn6o3VnXpsUe49ryef3CAT1pFpVWcALjdbvN8Kray+vn5+f7+XnCUExOrzfmfI6u+AQDghXoDoHcMSmNulV+1Tf3av2YyBm6V9rw2WrpWN5Tf0vvd7oM8qajIzBzAUNhvBJJVVFpda7iebln8MoL59aV1vgUDALAeMZ/aXN6owuhh6UPEUj9NvOoJzrqQrfwpRKzTKiqtojU0XG9Jl5jMwQYAoBCf/DCf92tPbV5/P/bj3/Zk5z68Tw2nejtPDdd5lUwDLHi5VVRaRQPlVIvMs1JitVlrgSsAAEtSdYy/71fTxP/SLX9V9TefNOSt3hSyjP9lLb6SC3jYKkKr7XYbQvj+/m7dvsisKvqrqk0Im3NR9IyuTgAAsGDP5COvijZa0dXIwqtnZn4tezR+7dld+30pRALGW8X0wLqu07okcdHAEMLPz0/84ng8TrZn7/Qbwtd5qdep9wUAgIWrssYU4UqlfzlD0zgXIVVapS/S+Pla2/i7nsK1hxrv+UcY4x2/lGbOn0Te1Jddu3d4wFzfRx6w2+1iR/bc8Xjs1lvN9+01yicGBu+GALzH3E+XfMzNaUEf2xNaxvwVDy8qN/5xwmWYdW/0E7JDpbWs3rXdu/bjUgyUbumdJLjJYpdrL8LAs8jvOJDiDT/r3h3oPvGB+17bbOBZPJA5du8yvPP5jw7Xf4Otxx9+yt39fODFH97V8cViz3x3KvceYx9W5ov2cit5mvdaRaVVVJ/Ff+52u0VOD2wlVqG8dxwAYG2MQ4gKORJanwxTB9hWehVnDj5c5eTz5wdMUru07MmeUBrvpD3me4JRcQrAx8z3dMmHOVSK9flKqzE1Sr13DA9VWl0rlQp9ZUH5P38vFxncnDtvjH8WrTqv8OZKq97JjLOotMo3eLjSavhXP3DHMY/T+/XwY86u0uraEVXOfhb4or3DSp7mvVZUabUGaR1fxzoAANekFKC0j6aF+GqaUFVfWdOrcGUWYSFaoc/4UqAPFw197Mf58A+LIbRajpLPowAAlMPn+ZG++pqgjxx1d1/kz8yKqM4rUKV/LuPX3S3aGpa/2ot5EWCFhFaLoo8VAABrkE/++oCvNK7u9L26+O6gvB6qqqp85187bn9sGmZ+9wWkPE++CNd8+MADhFYLEc+Xv6NPmQAAwAO+Ol3b09fXhuLdGGi4i1MJ3pepAYwntJq9dpnxVPsBAABvUGZhy+z6Xt1rGfVW7+PFgc8QWvUbOC+W9t7UNM1fVWm+DgAAL3czMkvVVX9PzBzsNRCLvCRR+vBMt2s/TjrGsNSmzUGyTkKrfjP6e1jS9RwAAHihSfpedWuv/qrq9zM78XEDyddLXvzWSpcPeGEoJjeBzxNazZ7m6wAA8IxWIJIvwHdtdD2QocTxeWvmYHCx+SF55FTgLFHg3YRW85Yu2ui/DgDAOrWKce69iJsKgh4oyRne/t8h+qJbX7U8XxhF4dotlX0O5c2EVjN2c5kSAABYvFl8bP6vBe3rWl8VuMBfazfyDKuQPWSkaxlu+QtfsjBCq7m6SLi9ZQAA8DZPljK93Oc/MA//xDFVRXkEECOq3tqru9KrZ7o1LTt0eMdTW/YrBsX639Q7wIN+Q4gTA5um8e4JADxgt9tVZ7vd7nA4fOa+zE5zaerdmV4rv0gvy12vz1cWYCV/VRX/Sz8oNdha8IS7T/bL562qS1PvDksgLe4xixDdXGIApjWL0yXX1HX9/f3dvX273dZ1/dr7OlRmbcyvb+SvON/smdZRebXXwOO0unf33j3fOFwvXOr+xHzJvO6TGri99SNaPy6vvdqcL1F3a69a+zDw1K69Pr2vYffzRevu116cgYKv7gt77ekPFPH13rd3g2v7+cAhOrx7N5v0X3tJx+/MJ107QE1fMQAAIABJREFUcsb8fT3w3cd2L6zgM2+Bx0YJVFrNUuy//ptd9Zp6jwCAOYmp03a7TWOJ4/EYQjidTjdrpp65Lzzs4aqNa015CvTVNPG//MZW7dX7tOrplMk877ESPCAnyetReMCp/zoAJSj8dMmAw+Hw8/MTOh/dUwnVwG/2gfs6VGatnEqrkBWqhHsqrcb86JuVVq3Kpt5ap+crrbo/tJVVbS63eW2l1cA2A482ptIq9L0a3a9D2ZVWrW9d29ub+1PgW6JKq0IUeGyUQKXVXEmsAIDHxEl8+/2+dftut8s3ePl9gXu1aq9+L/+7WX41UDDlszGRoITCWT1wTi5OOd5cAIAnpJjpw/eFxwzXfeTSmHm4HOZhY2pMXj63LuVW3aDqAzMHH/bal2LMMfCOWY1mSsKEBB89Sg6b/6pq48IIAAUo+XTJY8ZMD3zgvg6VWStnemBz2f68ZUyQMTDFbMz0wNCZptQ7PXB4DuC1mYPXnvLwk4o/6Pd8S+rdHrKQa/gZDfygMdMDu7+Lgdf25o8eM/tvzNPp7saT0wN7X4reJzvmeZX2ltg7E9D0wEmUdmwUQqXVnJR8FQUAmLuYOnWn/j1/38fqFIzdp/XueqUHpMTqgZ0ZWaTz2qf51tcwvRr/tQ3J/tDSB4c4i/CFrUVaf86PPak3vTLPl7ld64fVG88tMmJo/WpgcnpazUassYrNF6uq8iYCALzK4XCIQ4vtdnvvCoBj7ts85NlnxXPG/DrSoHQWo9Mn44zHfuKHD+nNlda3aQnC56+Cp+cy/hl1tyz2j/3zB0kJ8r/iMn8vrJlKq34Db1UT/vX+6r8OALxUmtYXQtjv93clVs/cl2V4YGA8bQHXSyKJ8j/Mf13O8/q9/G7Mrf5aVVrz96q6p6ZvRchlu/Yci6omK2pn+CShVb9y/hja9bdT7QcAsDi73e50OoVzkdRd7dWfuS9r9uQwe9pCmN6fXv5n6TyZalVa5QFWGMywVlh/9IzCD4kZedPKBsyI0Kp0KeNXZgUAvFD6DHA8Hu+NnJ65L/A+N4OSfz9QVNVX08xrIcJhrw2Jmie6pwGvpafVDMz35AEAlCl1oWqa5rHE6rH7AmN8pq7kq2lihpW+yMVJhX/nbkdTxTcf+9EvaTA/lerS1LsDr6TSqnTpTWcTQlBlCgA8LTaf2m63dV1/8r6weLOYFNZdHq5bVZRfNf+rqt9sImGZkz8+9rLfXFxvksxo2R24Up/4qXeEaQitSpf6JpZ5egAAZufn5yeEMKZIqq7rGE6lJuvj7wsz8rGmOUV15xkIONIswu63Up5173SQaXvwv0q+28OLd/VGSOX89mckvYxevXVaZhb7pKIi6uWt6wHAMhR1uuQuN8f96Td7OBxiSjX+M0PvhzSHyso9P/EqP4oeOKLG3H0gSbl299Yac/kjDOzkkzszsNm1lKR14/jMKD2vvM1utMkurodbHdxH/r7G7H/vsdTdpvuU730xw4jf1EAsdW0JwrcGdo+92Q7fa+AQfeBeD3jyb39eFv8EH6PSqmi6WQEAMHc+hs1X61P0Jv02+1YhTN56xb385RqBF/LX3qOcN0FlVgAUq5zTJYVzqPC8wiut0t3D4iqtQqfQslVc1vStQpiLNVljPtH07lhvFdVw7dL4qqjhn/7ySqvuzr+WSqu5W/wTfIxKq3IpswIAgI95ycfFMZ2b3jFH7E3tfsYEE61AqvdTTH7jSi7JqwiDlxBaFSq9rW9C8D4HAACzMCakSHU3z/+sD7emHhPEtDOpzh6mTzorSa+6xvfYSoRfrJbQqmhfTdO7ZgcAADBrL4yc8k7wBcrDqbzk6oGZJSU/zWH3Vl3dnGf6Eh+oBZO+8SShVYn+sqJif9UAAFCUMXMA6ZUCrN7E6ndWPVLWMAHw+ef4mfSNBRNaFWpz/sKJEAAAinJtcD7JoH2mccCYTljhMsZ633TCN5WqzbcuDMoxyze4d5v2fT9/U57pGQiANXCSYiSHCs/7zOqBz6z3N37HXrh6YGvLd6weeG3lvocfZHjL1mNeS7I2fYsYDv/EMavdDTzOyMNmYHHDfA3BMDrfvOtof/JPY2Cb/J9jKqd6j7rHTgRjfrmLsexn9zCVVv0GEvHPHEab8G/PQpVWAAAwiUmmAd5b9bPgTwr5VfxWgBX/OZfphHMPI9YwEZJiCa36TfUHmb/nelMAAIAJfWBAPrssYKopb8PNsGYRXU1FK3RmTWhVIosGAgAAL1f4OoNjxOa/X+dJar9925STYZXwas+3Ffoc95mXE1oV5L+3s6nf1wAAAMoXr/df9Ggfty7h+9q6dy0seZk8hmNVhFYFiVcJNot7UwMAgMdM0lWKWRgowxmzNGF+41x6Y8EKCa1K8VdVm/PXTskAwJJMvsQN8+UIKdaMfjXdZdmvRVQPRFerLTu6dwFBeIzQqiC/yqwAgCUyvAGKsunccrM3Vvq6VcaVr6y32gCrHLrOL4/Qqiz+qAAAgJKV0Fz8Sa3e5FXsijWiN1b4f3v3sq0okgUAFHrVH1XOxS8Tv0ycZ34TPYi+0SQvUUEjYO9BrSwvYhzkcTxEBGuMK2wHPb9YS76zzjPlP99uwEc1TVNVVfmjaZpvt+h/jKAGAABWNywwhf8NP4i+1KitxKCeii6+ZfjGf9v210+3rF9j/bOiP2X5pywTmRur/Nuba+vtQg/rQXGum/3tYHzFgXpaVVV1v9+7r5zP59Pp9N3SVb/74rfaAQAAu5PpPO7dMsrqbU55I7zZBWmLNz4s0DwcVxj+9/MPK5zakrnXkrrtT3lnZi1H6WlV13WoWN1ut7Zt27a93W5FUdzv9+8Wrdq2Dae53w45AABYVfu3bzdnqRzb/GHpjP8K1ajYJ2uqW1boihWKWX8mSl2fEferRDbgs2LjM20/z0rlUN9aKMfebreqquKLdV1fr9dhZ6tPngFNFAdAptL5wUDi7Cqs6809aqb7yZLVziz2whqGfameja63toehPbX+3sKvrefZJWc+tPti0dluSz5iaj0zfbuWf6HFbB+xXm+sXxOTvheDKd6Xt2F04SX/+2yM3TcWE79eRzfIa4ft8F0Lm53pdSfTZm/tEMMDY02qW7EqiqKu64+3pS+esD7ZXxQAAI4mzaGC5uQuBl/N/vTmd/+3bUenvvrwfFhb73j2bVZxiKJVKE5dLpepP32dihUAAGwq5R/PiRdrtq739XoArbvyNP0abMapilV4vTfY8Fs7syIUn3eIolWYzSp0s2qaJnS8qqqq1/Hq81J4tAQAAGQhza5Sq4h9UtIs2by2qXf8fW1h2I+h+2sxDNCJQwt7PySP0wci2WOE7RyiaBV19+/r9VoMZrn6iuOcYgAA4GXfrXrsoATz4W4ySW2ip4aqJfJd/9u2Re+xg7O9sbr/+MBvzCWb1AhB3neIHahbqzqdTrHLVXyeYK9u9XLt9qmN+ckTCgCsThrKQnYVsvD+jvraROzdPlYv/MJ/aiL2d1a+/E9vWjgR+wuz1788q/fD1RaPJmKP27/7+gttHv2I+LuyN8X78H+HwzCfquJ1Y5zZdKPT5K81EfvClmd63cm02Vs7Vk+rXn2qqqr7/X4+n4d7xtb7yl8nLLsmAADAkxLpErXQRlN3/b8PRFmGytT89FibWnGUqwoOwYGKVrGPVdQ0zbcGxMYxyY5DAACAoSWjz7b+iA+vZxX/K2N1hhaOlqt+Pzk31msxmoWKNx2oaDUzd1XTNB+e2So+/SGLewIAAACs5cOTPS2cHmuqK9bvn8GGf3orge0domh1Op3C9FXpCD2tHPAAAHAQ241ly2uUHInojftZOHiwO+l77Kvlhy3bOUTRqq7r8/ncNE3v9fjKh7tZfWAsMQAAkJTtakkzaw5lrGfHZ+2vCtbbDjuIqOfhfPAPPSg8lWVvZveu+At3OORw3guT03M0CY283VQ4GEYnYj+dTr161tYdNVWjAdiBpObvIGV2FbLwsacHDt+y1tMDZ5Yp1h6P9pmnBz58Ht+zay7eeILe/Gpn1vNOFMufkzj1gL+ZVc20eWb93bWF/w6rVMPy1vxszkueEflUyzO97mTa7K0doqdVURSXy+V6vZ7P5zgd+/V6DX8a9sDalG5WAABwcL2uTPHffrLSk85eMVNC6i/2zIRZ3b/q2MFQKgfAB9R1HQtVwbCPVbDpecHRCMA+pJNGkzi7Cln4Sk+r7nvDP17rq7Vk5bn0tOqufPc9rbr/u6RX0Xd7Wg3XNvysqQ6DZVlODS0sxnpmRcMfznpaHc3hNkpd1+EfVVVNTWW10b5ivC4AeyK1YiG7CllQtFr9c99feRZFq4frmY/iYR2qeLJoNbXmLYpWw8+aKVo9HB64cEzSr6J42PJMrzuZNntrNsqI7faVPz/T1+lpBUDupFYsND8DtL2Ir1vr1rKi1VoUrbqrLRIrWo2OB5zqeDX16VMr77041f1qpmdWXCDHi4vMatRR5rRKQTy8fxVFYXcEAA5D2kPK7J+kbN1q4/uG1aiFb3whimFXj+UzRA+X1HEkU4pWnxOLwY4WAADgHb2p3NMpauRlre2W+PYPO0mvbJR4m0f927blkkmin58JPq7/1aaxFUWrT3MYAADAnnylfrTFQEL2of2p7HzeB46FhdH1fncv7KKlf1aCFK0AAABepxjEqKcqODNzq705QrA74dTWtcu1HsW4uiW1p9HClv5ZX6doNW7maHntOFw++BYAAIDcPfXLcX7u+ZcNhwSm3OeuN5H8+yt86jELL0+hpX/WphStxm10GNt3AQAAmK/LrFK4Sbk+9RlvlgI37Z+1cP0oWm2rX9n9VjsAAADYo2GFa6bg9cUZr5bLqNy2yiMOVa9mKFptK1Z2f9sRAQAAeN5rlaZc6j47s/CHvxmEFlK02lY8s/wqiiKfajEAAPAxm840lEXPms878jZZcZdYvqoPT6eVfl8tnVoWUrTaVtu2uvwBAACwY+kUiY5cjtyl/3y7AUehYgUAAOSuLMtQFIj/yEUiJRXgKXpabcs4VQAAYDc+U/qJ5bCN+u+s8my+D8uuwaP6Typ7+8tdfYWkRk+rT9DNCgAAYKFYejhmDSKdoXara9s2hBb/8ZpYq1prhSRL0QoAAIBUbN2laMclIaJ9dEyjULTalLGBAAAAAK9RtNqcsYEAAMB36Vu0b5nOjr863ej2x0Tsm/jrTOGwAQAASEDuRY2p9mcRVBiyt0pTV1wVidPTahNt2/4uiqIofmdy+gAAAEhc7Em0jy5FfirCQ3pabSKeQH8VeloBAACswA+roa93OEq201OyDeMpilbjZsr2S3b63z//MKEVAABAItIvYaizQJei1bhVThMqVgAAxdu3AwGAY1K0AgBgWypTwAuOfOqItf7wj91vihf6lw1H/+mktkuKVuv7k/+MgAAAAHzLfPFlB5PQR6OxqD0ReXrgVowNBAAAYEWxe1Eo66Rf3Fn4wMcQ1+qfvtFq+SQ9rQAAAOA79j2obRiaKhJP0dNqZcYGAgAAsFDivYEWdpV6eeWxqvXCdkh807EKPa1W0ztadlsqBwAAYD3dOcVT63WVWns4Gj2tVhPHFf92YAMAAJCATbtKreudnlP7HmV5ZHparSYeXb+KonDAAAAAsMx2I92S/WX6mTKTIYS509NqNW3b/v7pZpXseQEAACBxx/w9dcyoe2wEehStAAAAIA/drkPGxOlItXuKVqvx3EAAAACAtZjTamX/HrvODQAAkLhuT6XijSFpa63nUELfqI36iPlG9kfRah26WQEAAGRhrVpGmjWRYxZuQpgHCfZQFK3WpJsVAAAAXzRfuNlfScu0VvumaDVuZqffwVENAADAASXye/a1Zuyv4sZDilbjntr7jQ0EAACATalSHZCnB67G2EAAAACAtehp9Za/RhFu8/gDAAAAjmw4LM4sThyEohUAAACkS/cIDkvR6i1t24YJrYwNBAAAXmBuaYApilYAAABfo0oFMEXR6i2eGwgA8NDM3Ct+rgO8zCmU3Ttu0WrFzrfGBgIAzPCzCgB4wX++3YDvqKrq200AAAAAYNIRi1Z1Xd/v9/fXY2wgAAAAn6cHKwdRHm1fb5rmfD7H/x0NvywXbRbPDQTgyBZeLsGuAl+X42G4RZsT3A69Jj3VwnTCebYlcfnejIfxxWfjmn9LOhtqXi7t/LDD9bQKFavb7fbthgAAAGyrLMtQF4j/OCbbIU3t377dHFJ0rInYw1RWl8vFnFYAAMDuKQQEtsNGYgVwxQedQdeBup+FgYGn06lpmmL2oFrSKy9OaGV4IADHpBM7C9lVgBfs/tTx/si4fDfR6qP5DA/csQNtlF6Var5o9dpHHGdjAoDUioXsKsALjnnqULR6+NfPrPArcmnnhx1lTqswHnD5VFbtI7/HFtuq9QAAAAAHc4g5req6vt/vp9NpramsDNwFAAAA2NT+u5+FqayKQWnpnTmtTGgFADqxs5BdBXjBMU8dhgc+/OtnVvgVubTzw44yPLD4ebhp1H2xrusXVqhiBQAAALCRQwwPBAAAAFJmHh6G9t/Tqqqq0ZnUw1/Dv1/raQUAsBtVVcUO6VVVvZAdNU1TlmXTNOs3DmDv4nig7sCgo/GgM4b0tHran6OeQQCAXYoTgEb3+/1+vzdN81QFyl1AgJep0WxB160d2H9Pq42Y0AoA2IdQsTqdTvHm9u12K4rifr8vr0OFhzVv10gAeJauWzugaAUAcFyxLNXtVFVVVahbXa/X+beHIYFlWT5cEgDgWcctWim1AgCEWtXlcum9XlVVd4H5twMAbOG4RavXmNAKANifWKJ6Vl3Xhl0A6zIlORCZiP0VJrQCAPZhqqtUfP3lehbAaxTBd2z1mdFNtb57ilYAAPSF2dmHwwZf81pfCb89AHZm9RO7K8XuGR4IAMD/1XUdakyn02n50wPntS9Z5aMBgHzpaTVu9H7g78+3AwDgU5qmCR2siqK4XC5rVawAAF6jp9W4mdt9JrQCAPanqqpQsTqdTrfbTcUKAPg6Pa2W+n/fKxO8AQD7EvOc2+1m5nUAIBGKVgAAhxZnsJp6kiAAwFcoWi3Vtu2fsiwMDwQAdiQMA1SxAgASpGgFAHBc1+u1KIolQwKbpgmFLdNdAQCfUZqbaagsRzbLn5+5HvS0AoBi4nJJdkafmNwVv+W6rkOFa+p7D6sazoplVwHYgrPrzvhCR3l64HNUrAAAAAA+QCVvxExPK0UrAAjcD2QhuwrAFpxdd8YXOsqcVgAAAJCNOLI7/EOlgx1TtAIAAIBsqFJxHOa0AgAAACA5ilaL/Hn0YB0AAAAAVqRo9QSzsAMAAAB8hqIVAAAAAMlRtAIAAAAgOYpW48q/dV/8bsMAAAAAjkDRalz7t+6L320YAAAAwBEoWj3m0YEAAAAAH/bPtxuQDY8OBAB4zcwEC7qxAwBTFK0AANiWyhQA8ALDAwEAAABIjp5Wc/7qyl6WbhICAAAAfIaeVgAAAAAkR0+rOW3bxkcHmogdAAAA4GP0tFpExQoAAADgkxStAAAAAEiOohUAAAAAyVG0AgAAACA5ilYAAAAAJEfRak58dCAAAAAAn6RoNa4sy/KnYvXr53+7LwIAAACwHUWrcW3btm3b/Xf03YYBAAAAHIGiFQAAAADJUbQCAAAAIDmKVgAAAAAkR9FqkkcHAgAAAHyLotUDv77dAAAAAIADUrR6wOMCAQAAAD7vWEWrqqrKH1VV1XX97RYBAAAAMKI8SE+ipmnO5/Pw9dPp1DRN78WyLNu2DXNa/XuM7QMAzwqXy2+3ggzYVQDgIZfLUUfpaRUqVqfTqf1xu92Korjf71P9rf5t2x1PaFUeY5r5I4R5hBiLY4Qpxt04QphHiJF1ldO+3bR17CaQeUcI8wgxFscIU4y7cYQwjxDjyw5Ryavr+nq9FoMJqmL3q97rscC540rnjkPrOkKYR4ixOEaYYtyNI4QZYjxCpKziCLvKEWIsjhHmEWIsjhGmGHfjCGHKrGYcoqdVGAB4uVx6r1dV1V0AAAAAgEQcomgVxBIVAAAAAIk7dPczwwO/3YrNHSHMI8RYHCNMMe7GEcLUiZ2nHGFXOUKMxTHCPEKMxTHCFONuHCFMmdWMf77dgG8KFavhsMGiMxHaUzOi2cMAAAAAVnGg4YFddV2HatTpdBp9emB4wmD8x0IfjgIAAABgrw5XtGqapizL8DDBy+UyMwW7vnkAAAAA33Ks4YFVVd3v9+Kng5Wp2QEAAADSdKCiVZyd6na7rViuerZD1lPLb7fws9Jp9nZhptPsRGJ8dvlEYty0JZnurs8u76t8Z+FnpbMH6l9MItLZzzM99jO9VKXzvT8lnWbbXd/kq0x24WelswfKrNZylKJVnMFqZjwgAAAAAIk4xJxWYap1FSsAAACAXByiaBWmXX9qSGDomVWWZRxUCAAAAMDHHKJoFVyv13LCt5sGAAAAwF+OMqfVs0yZBgAAAPBFhyhaqUABAAAA5OVAwwMBAAAAyIWiFQAAAADJKQ2dGzI1OwAsIYtgCZkVACwhsxpStAIAAAAgOYYHAgAAAJAcRSsAAAAAkqNoBQAAAEByFK0AAAAASI6iFQAAAADJUbQCAAAAIDmKVv9X13X5o6qqbzdnHU3TlGXZNM3UAr2oZ5ZMU1VV3fbXdT26WL5hNk3Ti3Gq8fnGOBSiGP1TpmF2m90z3GkzjTHo7bE7212nvsSpbzPTMIMDnnlYncwqx4NCZhXlG+NQiGL0T5mGKbPqyTTGqS9x6tvMNMzggGee1bS0bdu2oxvndrt9u13vOp1OM4GEv/ZcLpePNvFVt9tt9Fs7nU69JfMN83K5jMY4bHy+MQ7FWGb+lF2Yoy0fbX++MbYTjd/TITn1JY5GkW+YU2fXfZ95WN3oXiSzSpnMqivfGIdkVpnG2MqsZFb5xLgpRau27ewcMQWZOblnJF6YR1OrePDEw2B++dQMz9rDiEZfzCjMYTix8d3Fso6xp5tN9v6UdZjzP3KirGMcNjWG02181jHO6O20WYc5PLuONj7rGNmazGrJ8qmRWUVZx9gjs8o3RpmVzCqjGDeVd+qwltHzeL77x/Am0mgUw8t2+3MBGNbvUzOaYbSdQz2+km+YIcZhI4cR5RtjT+8uRO+vWYc5GtHUYlnH2DvbDHfjrGOcMkyd8w1zeBYNQuN3eeZhCzKrKJeDQma1y/ObzGoHMcqsgnzDlFm9T9Fq8iI9dWFL31OpVe/FqYMqNcODPOqFnG+YU3eQhmeufGPsid/dTGrVezGXMJ9KrXovZhHjTCMvl8vwerz87ekLje9dLPINc+raN5Ul9xbLIka2JrPqyuWgkFnt8vwms8o3RpmVzCrIIsatHTr4YOoivY/9Yyq1moluJhtLx0yX4G77sw6zd0GKents1jF2deMaRpR7mEtOJlnHOPNrpyvrGKfsbHedKjf0ftRlHSNbk1ktfEtSZFYyqyiLMGVWQdYxTtnZ7iqzep+nBxb3+70oiuFDbXbzmJtR4TEEozO9ZaFpmrZth99RfLxC+FPWYdZ1PXz6SdM0vT026xijENfpdJp6SlHWYXZ3y5nnoWQdY3e3bJom7L3Dx51kHeOo8CX2emFkHWbcLbsn2Lquw1cc/5p1jGxNZpUjmZXMKiMyqyjrGEfJrD7auEwoWtEXDpVMH655Pp+LwWluVF5h1nUdrsoxwCWpf0Yxhrhea2r6Yca2hYtTvBpdr9ep50/3pB9jFPbS6/V6vV7P5/PMs5l7Moqx63q9Fp2E46EswmzbtiiK+/0efwaEMG+3287OPPAZWR8UMquujGKUWc1LP8ZIZjUvizBlVm9StDqosNPv6Z5nXdfhKtW9p7SbMONtwKAb0Q5iDI3vzRXas4Mwi59uveFudox3T7trSJFPp1MYghEusefzOV5ldxBj19Suu7MwRx0hRnjK/g4KmdWn27QemVX4xw5ilFkFOwtz1BFifJmi1aHto17bNE0sV18ul6l+s1lrmuZ2u91ut+G1Ki7wlYa9L3SOPZ1OS87RmYZZ13VIqroxVlUVblyHXTfKNMbodrt1O7HHPba7TO4xBsMRJcMFPtmeFcWfqXH2hLCv7unMAxvZx0EhsypyjlFmJbPKlMyqyDnGTSlaHdTMZWz+ZJGaqqriLYjb7dbrSrqbMIuiqKqqqqqmaeI5Lr4+9Zb0Y2yaJiQWD0/QWYdZTDRvf7vrMEXufbM7iDEanXMhyDrM0LbT6dQ0TWxnXdd7OvPAFnZzUMisso5RZvVwmSCLGAuZ1Y+sw5RZvU/R6n92eRNp98qyDMdwuAVxkCN5+RjvjJR/6764y3h7dnO2mTkGdxNj1Js+czem4uqNtoCHZFY5klnthszq201Yh8xqB2RW71O0mpzYbN97TzgDdgfzjy6QstjNcvRhN0HWYXaTjBlZx7hc1mE2TTN/PgmNzzrGhc86yTrGrpBnTEW9gzAftnAHMbIdmdXMAimTWQVZx7hc1mHKrKKsY+ySWe0gxg21hxdneuu9Ho6Zy+XyjUatJoQWR88+/NPU1khN6E4Z8qp5+YY59d0NG59vjFNGW55vmEf4KkMjh4fknmLseniByDfMh1eN+SWziJGtyay6L+ZyUMis9n1+k1nNL5YgmVVPvmHKrN536OCjsB90D5Ld7BwzB0k4NfROhcNNkabl7cw3zNGWtz+N776eb4xTRo++fMOMd416rw8bn2+M7cTZZhhR1jFGDy8Q+YY5tbsOX883Rj5AZtVbPv2DQma1v0tVl8xqfrE0yay68g1TZvW+7FOHVcT53i6Xy+12i/+75HZT4mZSq7ZzkQ5RTx1RCSoeGS6cdZiXH6MBtjnHOGqq8fmGeYSvsnvm3GuMwcLf3vmGObO79q4m+cbI1mRW2R0UxSPDhbMOc8eX41FTjc83zCN8lTKrnnzDnNldZVZLHD3+6DJ4TsEO8qr2UWoVTxBdn23gi4bNnom1FwPMAAAJvklEQVQi3zBHWz66Z+Yb46ip9ucb5kG+yoUn0qxjbBcPosk6zJgkdc0MxMgxRrYms8rroBg2eyaKfMM8yOV4aKr9+YZ5kK9SZtWVdZgyq3eU7YKr1HHEOfzDQ3C/2ZQPinMZ7jvqfMPsTjY53/h8Y3xKvmEub3m+MRaLT6RZx7hcvmE687AKmdWOo843TOe3nnzDlFl1ZR3jcvmG6czzMkUrAAAAAJLzn283AAAAAAD6FK0AAAAASI6iFQAAAADJUbQCAAAAIDn/fLsBAJN2+eyMENQ+Iqrr+mEg8ZE3Pc9ugantFtb/1NqmmjTawnS+r26zV2/PpisHIBEyq8TJrD5JZpWNFiBVp9Npfyery+VSFMXpdPp2Q951u92WBDJ/DbpcLu98XFzPU9tz4fUxrHNhmB/w2nZLYeUAJEJmlTKZ1YfJrHKhpxWQrvv9Hv8d7sl8ry3ZK8uyKIp2cWKxrtPp1P36mqYJX+71eo13fZ8V33W73eLKQ5jdV4ZCdjtsTO/1NPe3kOqt3raw2rquuwcdADsjs1qRzCqSWU2tVma1jm9XzQDG9a5z7lG8afVz/lP3A2+32/BP8X7v6F8fGr2z+sIKQyDpXxA/0MiwSR1rALsks1qXzGqKzCqSWa3CROxAoq7Xa1EUl8slnO7D/7InTdOE7Op8Pr+8kjTv2gFAamRWuyezYpcUrYAUxf7JdV3HiQyHXZ27PdurqirLsizLqqqG80GGeR/DCmcWC/NfxjWXZdn70PlPCW8ftjOuNq65+94YY/zQsiy7UXc/dLS/dwxqtGHdTx/+tbv+7uf29JoxuswLYjjdf8Qwwyd2N05sXlVVIdu+Xq/dST3DX6e+iNdauPr3tXCbLxE/4uEX1PtQKSnAocisZFaFzGoBmVWKvt3VC2BEr39yuGs07C8dXu91dw+6C8fu1t35R6Nuh+f4yvCv3Ren3h5bMhPLsO93eMtoFKMvdlc+1aruYjN/erjputt5tG3vdGLvrjx2nI79qOOHjk7bOWzzMJaFXdnnO7Gv+321i7f50HBtsW0v7ye9TaQTO8BeyaxmXuyuXGYVl5FZzXyuzOqTFK2AFPXO+6MpS/v3Jb93be6+0r2u9FKc3mp7l8lh1tW9BMZPHzZp+NFxgalL9cMGt52L+kyr4nu7jR/NHuKLw8/trnDh9pwyegnv6uXN3SzhdDr1Utthw3p5wMOPG3ottXrt+1q4zUfNfIMP94rhMqM/V6RWAHvVuz7KrLpLyqy6S8qsluwVw2VkVttRtAKSM3qpG71kTs032VvDzNWr9/Z4oeotNnVDcnjF6q1z2LyZ+0vDlRSztxbbscttO5tzzMQ+teRo855NCGZynV5QMbUa/U7TSa16nzuaZ099Xw+3+ajhMsNEbfRz5xPr4RulVgA7I7OaaonMSmY12nKZVVLMaQUkJwxE73WcDv87NUa9N4y8O83BcM1d4VrSe33YJTg8rXY4ij4s2Xs9vHg+n0MzTn8/k3hKb5nR/vnD9YRT+cOVD8U2D9fZ3dTdKTB6LRnt2f6y4fOAE58aYLR589/Xwm3+guG3U4xt0t4yL+88AORFZlXIrGRWz5BZJeWfbzcAoC9cEobZ0v1+H71ajF7jT6fT/X7vTiE5XGfxM+tkb7XLL+qjb6/rumma2NqFk1aukkmEWTOXPA8otOp0Os03Ly42/FP4Rl5q5ojRTDplL3xfC7f56mJTy7K8XC5VVSWetgKwLpnVy2RWHyOzYoqiFZCWeNW5Xq+jKUIvWyomLnLhwv/wGrbkGjOzkqm3N01TlmXx98yjG6nr+uWnVt/v96mHInejntrCqzwtOz7Oprfy99ecpoXbfF232y18aPewulwu7zxeB4AsyKyeJbPKi8xq9xStgLTEc/3wjlC4+xTutn22UU+LmcH5fN60n3BI4IqfrvLB8mRrpoN99/Wppws/19YJK95UzMLCbb6u0GW9+RG2+fV6Df+70YcCkAKZ1VNkVtmRWe2eohWQlnDSv91uw8tMyBiGV+KZC//Da9WS68rMSqY+OjQydKSv63qjuy5xte9kb/NtW+um35SpmR327VvBdruvN01zPp+HAz0A2BmZ1XIyq0zJrPbNROxAQuIlZ/RcH//auzKN3lAanb5hmAmNTk06ZeHbQ2/h2+0Wlg93XZas/1lhtcPJTZeopmeUrKqqLMtum0cXWyXlCtvqtRCy89Q2X1Fd1+Euce8Tw37rfiDAjsmsniKzyovM6iAUrYCEzExOGUxdDKbG7fde7414jzfuHt6fiY+tefj28IlhUsapN26q+ZkudP6SGa+svdhjUKH9cbHelnzzjlbTNHVdxx74B7kZOLUxe9t8C/f7fZgKb/2hAHydzOp9MqtkyawOQtEKSMjDRGfqjsr9fg93PMIFO/aEH66hLMuwWFVVIeNZ8uDk2KTw9pAZDN9e13UvhKl7mKsI67xer6FJ3VYVRTE6W2p3pH18e9h0oTPz8AZdWOx+v3c33bM3A8/nc9lxPp/jGtadUTXmi2ne5oobc36bb/Ghxd97b8xrpVYAOyazeorMakhmNfWhhczqY1qANMRLy/xiYZnL5dK2bbiOXi6X4S3E2+0W3xIv3sOr1+l0Gq68+96u4aeEZjx8e/f10Jju545GHUPrvjh877BJ4a/Dk/zomX80rel96Ohip9Np2JhRw/XPfFD7sxtMtaH7caNLdr/iqe9xKrqZv671fbWLt/nQ8HOnWj58fSp/7b1xauMDkCOZ1fCzZFYzscusZFZpKtstH74AsKmqqu73e3i+bLzT1Z0TMQgzIxY/02rGJV+4RzfzKd8y2qTRGKduTi4MKsHYR73z/X7MCxsz3MF756rdvRs8+rl1XV+vVw9sBjgsmVUhsxqQWT380KnPlVmtQtEKyFg3tZpZrJdaQY7eT60ekloBHJzMiuOQWeXCnFYAAAAAJOefbzcAAFgq9DwPk7auvtrRh0YDAOyVzCp9ilYAkI2Q/YSH46y+WgCAQ5FZpc+cVkDeFl5j4tSMGzcHACBjMisgKYpWAAAAACTHROwAAAAAJEfRCgAAAIDkKFoBAAAAkBxFKwAAAACSo2gFAAAAQHIUrQAAAABIjqIVAAAAAMlRtAIAAAAgOYpWAAAAACRH0QoAAACA5ChaAQAAAJAcRSsAAAAAkqNoBQAAAEByFK0AAAAASI6iFQAAAADJUbQCAAAAIDmKVgAAAAAkR9EKAAAAgOQoWgEAAACQHEUrAAAAAJKjaAUAAABAchStAAAAAEiOohUAAAAAyVG0AgAAACA5ilYAAAAAJEfRCgAAAIDkKFoBAAAAkBxFKwAAAACS81+JaVXXKlo4RQAAAABJRU5ErkJggg==", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import McGenerator\n", + "imp.reload(McGenerator)\n", + "\n", + "attenuationLengthDict = { ### https://henke.lbl.gov/optical_constants/atten2.html\n", + " '5keV': 18.0101, '10keV':133.706, '15keV':442.607, '20keV':1037.80, '25keV':1991.45,\n", + "}\n", + "element = '15keV'\n", + "mcConfig = {\n", + " 'element': element,\n", + " 'energy': 15000, ### eV\n", + " 'attenuationLength': attenuationLengthDict[element],\n", + " 'sensorCode': 'Moench040',\n", + " 'nTotalIncident': 10000000,\n", + " 'nThread': 16,\n", + " 'Roi': [50, 350, 50, 350], \n", + " 'clusterWidth': 5, ### pixel\n", + " 'pixelSize': 25, ### um\n", + " 'noiseEneFrame': np.ones((400, 400)) * 0, \n", + " 'shotNoiseFactor': 0.0,\n", + " 'calibrationNoise': 0.0,\n", + " 'sensorThickness': 650, ### moench040: 650 um, 34.7 V\n", + " 'T': 273 + 20.0 + 14.8,\n", + " 'zBins': 128,\n", + " 'hv': 150,\n", + " 'depletionVoltage': 34.7,\n", + " 'simulationMethod': 'simu2', ### 'simu1' (cpu) or 'simu2' (gpu)\n", + " 'resultsPath': '/home/xie_x1/MLED/McSimulation/BatchResults',\n", + " 'sampleOutputPath': '/home/xie_x1/MLXID/McGeneration/Samples'\n", + "}\n", + "\n", + "McGenerator.init(mcConfig)\n", + "McGenerator.parameterization().Draw()\n", + "McGenerator.generateFromParameters()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(625000, 4)\n", + "2.96790016863627 2.287150706089626 15088.160000000002\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAGiCAYAAADTMXDkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAVyklEQVR4nO3df4yVhbng8ecwwEH5MTIowiw/SrVq5Ye7grVjrFWx5HINkTS9qb3WsDV3E7poMSxpi95EydKOae6mMaElYhurtzXYpqJmtxLprUBNQxewrCy1Fi1bpnsRrnqZgUEOMLz7R6+zUkDnDDPPmRk+n+SNnuP7vufJkeGb9z3vvKdUFEURANDLBtV6AADODYIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0CKqoLz4IMPRqlUOmkZN25cb80GwAAyuNoNpk6dGj//+c87H9fV1fXoQAAMTFUHZ/DgwY5qAKha1cHZtWtXNDY2RrlcjmuvvTa++c1vxkc/+tEzrl+pVKJSqXQ+PnHiRLzzzjsxZsyYKJVK3ZsagJooiiIOHjwYjY2NMWhQdZcBlKr5eoLnn38+Dh8+HJdddlns27cvVqxYEb/73e9i586dMWbMmNNu8+CDD8by5curGgqAvq2lpSUmTJhQ1TZVBecvtbe3xyWXXBJf/epXY8mSJadd5y+PcFpbW2PSpElxffx1DI4h3X1pAGrgeByLl+JnceDAgaivr69q26pPqb3f8OHDY/r06bFr164zrlMul6NcLp/mhYfE4JLgAPQr/3aI0p2PRM7q93AqlUq8+uqrMX78+LPZDQDngKqCs3Tp0ti4cWPs3r07fv3rX8fnPve5aGtriwULFvTWfAAMEFWdUvvTn/4UX/jCF+Ktt96Kiy66KD75yU/G5s2bY/Lkyb01HwADRFXBWbNmTW/NAcAA515qAKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJDirILT3NwcpVIp7r333h4aB4CBqtvB2bJlS6xevTpmzJjRk/MAMEB1KziHDh2KO+64Ix599NEYPXp0T88EwAA0uDsbLVq0KG699da45ZZbYsWKFR+4bqVSiUql0vm4ra2tOy8JZ1Q39fJaj9A/HD1W6wn6hY5df6j1CANW1cFZs2ZNvPzyy7Fly5Yurd/c3BzLly+vejAABpaqTqm1tLTE4sWL44c//GEMGzasS9ssW7YsWltbO5eWlpZuDQpA/1bVEc62bdti//79MXPmzM7nOjo6YtOmTbFy5cqoVCpRV1d30jblcjnK5XLPTAtAv1VVcGbPnh07duw46bkvfelLccUVV8TXvva1U2IDAO+pKjgjR46MadOmnfTc8OHDY8yYMac8DwDv504DAKTo1mXR77dhw4YeGAOAgc4RDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkqCo4q1atihkzZsSoUaNi1KhR0dTUFM8//3xvzQbAAFJVcCZMmBAPPfRQbN26NbZu3Ro333xz3HbbbbFz587emg+AAWJwNSvPmzfvpMff+MY3YtWqVbF58+aYOnVqjw4GwMBSVXDer6OjI37yk59Ee3t7NDU1nXG9SqUSlUql83FbW1t3X/KcU3fhmFqP0C/8bP1TtR6hX3j16OFaj9AvLJn6mVqP0KcNKo5GtHdz22o32LFjR4wYMSLK5XIsXLgw1q5dG1deeeUZ129ubo76+vrOZeLEid2bFIB+rergXH755bF9+/bYvHlzfPnLX44FCxbEb3/72zOuv2zZsmhtbe1cWlpazmpgAPqnqk+pDR06NC699NKIiJg1a1Zs2bIlHn744XjkkUdOu365XI5yuXx2UwLQ75317+EURXHSZzQAcDpVHeHcd999MXfu3Jg4cWIcPHgw1qxZExs2bIh169b11nwADBBVBWffvn1x5513xt69e6O+vj5mzJgR69ati898xlUdAHywqoLz/e9/v7fmAGCAcy81AFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQbXegCghx0rIn5/NGLv8YhSRIwfHHHZ0FpPBYIDA0JRRLx4OEqPt0ZsfDdKleLk/zysFBM/XY537jw/2q8vR5RKNRqUc5ngQH/XcixKS/dHadO7UUwbGsXXG6K4eljEpCERRUT88VjEb47E0B+3xke++E4cvKkc//zQBXF8XF2tJ+ccIzjQn219N0pf3BsxfFCc+NH4iJvOP/XoZfzgiE+eF2/cNTRG/FMlGu87EJf81b/EH59oiCMznGojj4sGoL967WiU/nZvxBVDo/jFxIibh3/wqbJSKQ7dMizeeGFsHP1IXUz+4tsxdPfxvHk55wkO9EfHiyh9ZV/EuLooftgYUd/102MdFwyKPz4+JjouGBSNSw9EdBQfug30BMGB/uhHbRH/uxLFwxdHjKj+x/hE/aD453+4IIZvPRoXPP1uLwwIpxIc6G+KIkqPtUb81fCI/zCs27s5/IlyHLyxHA1PtPfgcHBmggP9ze+PRem1o1H87aiz3tW/3n5+nPfKsRiyx2c59D7Bgf7mfx358z+v6f7RzXsOX/Pnq9TO23HsrPcFH0ZwoJ8p/Z9jUYyvixh19r9H03FhXRwfPcjVaqQQHOhvjkfE4J67U0AxOKLU0WO7gzMSHOhnioZBEW919MjlzKVKEXUHTsTx0f4qoPf5Uwb9zfRylN4tInYdPetdlV87FoOORRyZNqQHBoMPJjjQ3/z7YVGcX4rSfz901ruq/x9HomNUKY5cKTj0PsGB/mb4oIjPjYz4x7aI9hPd3s2gthNxwVOH41//5vwohrl7NL1PcKAfKv7z6IiDJ6L0jbe7vY9x/7UtSseKePvvRvTgZHBmggP90eQhUdw/5s93HPhRa9WbN/ygPUb/+HC8+fej4nijrykgh68ngP7qrvooXj8Wg5b+SxR/PB7Ff2mIKH/wqbHSkSLG/kNbXPhoe7z1d8PjwO3nJw0LggP9V6kUxTcvjGJcXZT+2ztRWncoioWjI+aPiDj/L05etJ+I0T9ujzGPHooh/7cj3rx/VLz9nz7k6wyghwkO9GelUsTihig+MzxKzW9Haen+iK/tj7h86Mnf+Pn7ozG+iDg4e1i0rG6IymWuSiOf4MBAcGU5in9sjNhzLOLFw1F6pRKx93hEKSJmDoviP9bH658qxbF/50ee2vGnDwaSSUMiFtTH6e5BcOzo4fRx4P1cpQZACsEBIIXgAJBCcABIITgApBAcAFJUFZzm5ua45pprYuTIkTF27NiYP39+vPbaa701GwADSFXB2bhxYyxatCg2b94c69evj+PHj8ecOXOivb29t+YDYICo6hc/161bd9Ljxx57LMaOHRvbtm2LG2644bTbVCqVqFQqnY/b2tq6MSYA/d1Z3WmgtfXPt0VvaGg44zrNzc2xfPnys3mZc1ZxyJFjV1z+2JdrPUK/UHfYjTq7YlLlf9Z6hD6tKI51e9tuXzRQFEUsWbIkrr/++pg2bdoZ11u2bFm0trZ2Li0tLd19SQD6sW4f4dx9993xyiuvxEsvvfSB65XL5SiXy919GQAGiG4F55577onnnnsuNm3aFBMmTOjpmQAYgKoKTlEUcc8998TatWtjw4YNMWXKlN6aC4ABpqrgLFq0KJ588sl49tlnY+TIkfHmm29GRER9fX2cd955vTIgAANDVRcNrFq1KlpbW+PGG2+M8ePHdy5PPfVUb80HwABR9Sk1AOgO91IDIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBRVB2fTpk0xb968aGxsjFKpFM8880wvjAXAQFN1cNrb2+Oqq66KlStX9sY8AAxQg6vdYO7cuTF37twur1+pVKJSqXQ+bmtrq/YlARgAqg5OtZqbm2P58uW9/TID0okjR2o9Qr/wkb/fXOsRGECKoqj1CH1aURzv9ra9ftHAsmXLorW1tXNpaWnp7ZcEoA/q9SOccrkc5XK5t18GgD7OZdEApBAcAFJUfUrt0KFD8frrr3c+3r17d2zfvj0aGhpi0qRJPTocAANH1cHZunVr3HTTTZ2PlyxZEhERCxYsiB/84Ac9NhgAA0vVwbnxxhtdNghA1XyGA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApBAeAFIIDQArBASCF4ACQQnAASCE4AKQQHABSCA4AKQQHgBSCA0AKwQEgheAAkEJwAEghOACkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgApuhWc7373uzFlypQYNmxYzJw5M375y1/29FwADDBVB+epp56Ke++9N+6///74zW9+E5/61Kdi7ty5sWfPnt6YD4ABolQURVHNBtdee21cffXVsWrVqs7nPv7xj8f8+fOjubn5lPUrlUpUKpXOx62trTFp0qS4Pv46BseQsxgd/k2pVOsJGEiq+yvxnHM8jsVL8bM4cOBA1NfXV7dxUYVKpVLU1dUVTz/99EnPf+UrXyluuOGG027zwAMPFBFhsVgslgG0vPHGG9XkoyiKohgcVXjrrbeio6MjLr744pOev/jii+PNN9887TbLli2LJUuWdD4+cOBATJ48Ofbs2VN9Hc8hbW1tMXHixGhpaYlRo0bVepw+yXvUNd6nrvE+dc17Z6kaGhqq3raq4Lyn9BenMIqiOOW595TL5SiXy6c8X19f739qF4waNcr79CG8R13jfeoa71PXDBpU/TVnVW1x4YUXRl1d3SlHM/v37z/lqAcA3q+q4AwdOjRmzpwZ69evP+n59evXx3XXXdejgwEwsFR9Sm3JkiVx5513xqxZs6KpqSlWr14de/bsiYULF3Zp+3K5HA888MBpT7Px/3mfPpz3qGu8T13jfeqas3mfqr4sOuLPv/j5rW99K/bu3RvTpk2Lb3/723HDDTdU/eIAnDu6FRwAqJZ7qQGQQnAASCE4AKQQHABSpAbH1xp8uE2bNsW8efOisbExSqVSPPPMM7Ueqc9pbm6Oa665JkaOHBljx46N+fPnx2uvvVbrsfqcVatWxYwZMzp/c76pqSmef/75Wo/V5zU3N0epVIp777231qP0KQ8++GCUSqWTlnHjxlW1j7Tg+FqDrmlvb4+rrroqVq5cWetR+qyNGzfGokWLYvPmzbF+/fo4fvx4zJkzJ9rb22s9Wp8yYcKEeOihh2Lr1q2xdevWuPnmm+O2226LnTt31nq0PmvLli2xevXqmDFjRq1H6ZOmTp0ae/fu7Vx27NhR3Q6qvt1nN33iE58oFi5ceNJzV1xxRfH1r389a4R+JyKKtWvX1nqMPm///v1FRBQbN26s9Sh93ujRo4vvfe97tR6jTzp48GDxsY99rFi/fn3x6U9/uli8eHGtR+pTHnjggeKqq646q32kHOEcPXo0tm3bFnPmzDnp+Tlz5sSvfvWrjBEYwFpbWyMiunX32nNFR0dHrFmzJtrb26OpqanW4/RJixYtiltvvTVuueWWWo/SZ+3atSsaGxtjypQpcfvtt8cf/vCHqrbv1t2iq9WdrzWAriiKIpYsWRLXX399TJs2rdbj9Dk7duyIpqamOHLkSIwYMSLWrl0bV155Za3H6nPWrFkTL7/8cmzZsqXWo/RZ1157bTzxxBNx2WWXxb59+2LFihVx3XXXxc6dO2PMmDFd2kdKcN5TzdcaQFfcfffd8corr8RLL71U61H6pMsvvzy2b98eBw4ciJ/+9KexYMGC2Lhxo+i8T0tLSyxevDheeOGFGDZsWK3H6bPmzp3b+e/Tp0+PpqamuOSSS+Lxxx8/6TvPPkhKcHytAb3hnnvuieeeey42bdoUEyZMqPU4fdLQoUPj0ksvjYiIWbNmxZYtW+Lhhx+ORx55pMaT9R3btm2L/fv3x8yZMzuf6+joiE2bNsXKlSujUqlEXV1dDSfsm4YPHx7Tp0+PXbt2dXmblM9wfK0BPakoirj77rvj6aefjl/84hcxZcqUWo/UbxRFEZVKpdZj9CmzZ8+OHTt2xPbt2zuXWbNmxR133BHbt28XmzOoVCrx6quvxvjx47u8TdoptbP9WoNzxaFDh+L111/vfLx79+7Yvn17NDQ0xKRJk2o4Wd+xaNGiePLJJ+PZZ5+NkSNHdh4519fXx3nnnVfj6fqO++67L+bOnRsTJ06MgwcPxpo1a2LDhg2xbt26Wo/Wp4wcOfKUz/+GDx8eY8aM8bng+yxdujTmzZsXkyZNiv3798eKFSuira0tFixY0PWdnP3Fcl33ne98p5g8eXIxdOjQ4uqrr3YZ62m8+OKLRUScsixYsKDWo/UZp3t/IqJ47LHHaj1an3LXXXd1/rxddNFFxezZs4sXXnih1mP1Cy6LPtXnP//5Yvz48cWQIUOKxsbG4rOf/Wyxc+fOqvbh6wkASOFeagCkEBwAUggOACkEB4AUggNACsEBIIXgAJBCcABIITgApBAcAFIIDgAp/h92ujBYogcxHAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "samples = np.load('/home/xie_x1/MLXID/McGeneration/Samples/15keV_Moench040_150V_0_samples.npy')\n", + "labels = np.load('/home/xie_x1/MLXID/McGeneration/Samples/15keV_Moench040_150V_0_labels.npy')\n", + "print(labels.shape)\n", + "import matplotlib.pyplot as plt\n", + "idx = 6\n", + "plt.imshow(samples[idx], origin='lower', extent = (0, samples.shape[1], 0, samples.shape[2]))\n", + "x,y,z,e = labels[idx]\n", + "print(x,y,e)\n", + "plt.scatter(x, y, s=200, facecolors='none', edgecolors='r')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pytorch_nightly", + "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.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/McGenerator.py b/McGenerator.py new file mode 100644 index 0000000..2b702f2 --- /dev/null +++ b/McGenerator.py @@ -0,0 +1,214 @@ +import os +import numpy as np +from ROOT import * +from multiprocessing import Pool +from array import array +from scipy.stats import gennorm + +EnergyPerPair = 3.62 # eV +FanoFactor = 0.13 +pars_alpha, pars_beta = None, None + +_cfg = {} ### global config + +def init(cfg): + global _cfg + _cfg = cfg + +def parameterization(): + sensorCode = _cfg['sensorCode'] + zBins = _cfg['zBins'] + biasVoltage = _cfg['hv'] + element = _cfg['element'] + simulationMethod = _cfg['simulationMethod'] + sensorThickness = _cfg['sensorThickness'] # cm + energy = _cfg['energy'] ### eV + reultsPath = _cfg['resultsPath'] + + print(f'Processing {sensorCode} {biasVoltage}V {element} {simulationMethod} with {zBins} bins') + xsFiles = [f'{reultsPath}/{sensorCode}_{biasVoltage}V_xs_{element}_{i}_of_{zBins}_{simulationMethod}.npy' for i in range(zBins)] + ysFiles = [f'{reultsPath}/{sensorCode}_{biasVoltage}V_ys_{element}_{i}_of_{zBins}_{simulationMethod}.npy' for i in range(zBins)] + ggdParList, ggdParUncertList = [], [] + + for i in range(zBins): + xs = np.load(xsFiles[i]) + ys = np.load(ysFiles[i]) + coords = np.concatenate((xs, ys)) + _rms = np.std(coords) + _binWidth = _rms / 20 + _nBin = int(100 / _binWidth) + _h1 = TH1F('_h1', '', _nBin, -_rms*5, _rms*5) + _h1.FillN(len(coords), array('d', coords), array('d', [1]*len(coords))) + def ggd(x, par): + beta = par[0] + alpha = par[1] + coef = beta / (2 * alpha * TMath.Gamma(1 / beta)) + return coef * TMath.Exp(-TMath.Power((abs(x[0] - 0) / alpha), beta)) + f1 = TF1('f1', ggd, -5 * _rms, 5 * _rms, 2) ### par[0]: beta, par[1]: alpha + f1.SetParLimits(0, 2, 5) + f1.SetParLimits(1, .5, 30) + f1.SetParameters(2., _rms*1.5) + _h1.Scale(_h1.GetNbinsX()/(10*_rms) /_h1.Integral()) + _h1.Fit(f1, 'Q') + ggdParList.append((f1.GetParameter(0), f1.GetParameter(1))) + ggdParUncertList.append((f1.GetParError(0), f1.GetParError(1))) + del _h1 + + z0List = np.linspace(0, sensorThickness, zBins+1) + z0List = (z0List[:-1] + z0List[1:])/2 + + arr_beta = np.array([ggdParList[i][0] for i in range(len(ggdParList))]) + arr_betaUncert = np.array([ggdParUncertList[i][0] for i in range(len(ggdParList))]) + arr_alpha = np.array([ggdParList[i][1] for i in range(len(ggdParList))]) + arr_alphaUncert = np.array([ggdParUncertList[i][1] for i in range(len(ggdParList))]) + arr_z0 = np.array([z0List[i] for i in range(len(ggdParList))]) + arr_t = np.array([getDriftTime_allpix2_8p23(z0) for z0 in z0List]) + + c = TCanvas() + c.SetCanvasSize(1600, 800) + c.Divide(2, 1) + # c.SetCanvasSize(800, 1375//2) + c.SetTopMargin(0.05) + c.SetRightMargin(0.05) + + global gr_AlphaT, gr_BetaT + pad2 = c.cd(1) + gr_AlphaT = TGraphErrors(len(arr_z0), arr_t, arr_alpha, array('d', np.zeros(len(arr_z0))), arr_alphaUncert) + gr_AlphaT.SetTitle(';Approximated Drift Time [ns];#alpha') + func_alpha = TF1('func_alpha', '[0] + [1] * sqrt((x)) + [2] * x + [3] * x^2') + func_alpha.SetLineColor(kRed+1) + gr_AlphaT.Fit(func_alpha, '') + # set fit line color + gr_AlphaT.Draw() + print(f'Alpha fitting chi2/NDF = {func_alpha.GetChisquare()/func_alpha.GetNDF()}') + + c.cd(2) + gr_BetaT = TGraphErrors(len(arr_z0), arr_t, arr_beta, array('d', np.zeros(len(arr_z0))), arr_betaUncert) + func_betaT = TF1('func_betaT', '[0]*(x-[1])^[2]+ [3]*exp([4]*x) + 2') + func_betaT.SetParameters(1.3, -4, -0.4, -0.6, -4.) + func_betaT.SetLineColor(kRed+1) + gr_BetaT.Fit(func_betaT, '') + gr_BetaT.SetTitle(';Approximated Drift Time [ns];#beta') + gr_BetaT.Draw() + print(f'Beta fitting chi2/NDF = {func_betaT.GetChisquare()/func_betaT.GetNDF()}') + + ### save parameters to file + global pars_alpha, pars_beta + pars_alpha = [func_alpha.GetParameter(i) for i in range(func_alpha.GetNpar())] + pars_beta = [func_betaT.GetParameter(i) for i in range(func_betaT.GetNpar())] + np.save(f'ggdPar_{sensorCode}_{biasVoltage}V_{element}.npy', pars_alpha + pars_beta) + + return c.Clone() + +def singleProcess(thredIdx): + energy = _cfg['energy'] ### eV + Roi = _cfg['Roi'] ### [x1, x2, y1, y2], for noise sampling from measured noise map + sensorThickness = _cfg['sensorThickness'] # cm + pixelSize = _cfg['pixelSize'] ### um + clusterWidth = _cfg['clusterWidth'] ### um + attenuationLength = _cfg['attenuationLength'] # cm + global pars_alpha, pars_beta + func_betaT = TF1('func_betaT', '[0]*(x-[1])^[2]+ [3]*exp([4]*x) + 2') + func_betaT.SetParameters(pars_beta[0], pars_beta[1], pars_beta[2], pars_beta[3], pars_beta[4]) + func_alphaT = TF1('func_alpha', '[0] + [1] * sqrt((x)) + [2] * x + [3] * x^2') + func_alphaT.SetParameters(pars_alpha[0], pars_alpha[1], pars_alpha[2], pars_alpha[3]) + + _frameWidth = 9 + nEvents = _cfg['nTotalIncident'] // _cfg['nThread'] + samples = np.zeros((nEvents, clusterWidth, clusterWidth)) + labels = np.zeros((nEvents, 4)) ### x, y, z, energy + ### set random seed + np.random.seed(thredIdx) + for i in range(nEvents): + ### z0 from 0 to sensorThickness + while True: + z0 = np.random.exponential(attenuationLength) + if 0 < z0 < sensorThickness: + break + ### drift time, alpha, beta + t = getDriftTime_allpix2_8p23(z0) + alpha, beta = func_alphaT.Eval(t), func_betaT.Eval(t) + ### nChargeCarrierPairs + nExpectedPair = round(energy / EnergyPerPair) + nPair = np.random.normal(nExpectedPair, np.sqrt(energy * FanoFactor / EnergyPerPair)) + nPair = round(nPair) + ### sample the incident position + x_center = np.random.uniform(pixelSize*_frameWidth//2, pixelSize*_frameWidth//2 + pixelSize) + y_center = np.random.uniform(pixelSize*_frameWidth//2, pixelSize*_frameWidth//2 + pixelSize) + ### sample pair positions + try: + x_pairs = gennorm.rvs(beta = beta, loc = x_center, scale = alpha, size = nPair) + y_pairs = gennorm.rvs(beta = beta, loc = y_center, scale = alpha, size = nPair) + except Exception as e: + print(f't = {t:.3f} beta = {beta:.3f}, alpha = {alpha:.3f}') + continue + x_pairs = x_pairs//pixelSize + y_pairs = y_pairs//pixelSize + ### put the pairs into the frame + _carrierArray = np.zeros((_frameWidth, _frameWidth)) + np.add.at(_carrierArray, (y_pairs.astype(np.int32), x_pairs.astype(np.int32)), 1) + _energyArray = _carrierArray * EnergyPerPair + ### spread noise + x_pedetal = np.random.randint(Roi[0], Roi[1] - _frameWidth) + y_pedetal = np.random.randint(Roi[2], Roi[3] - _frameWidth) + _pedestalNoiseFrame = _cfg['noiseEneFrame'][y_pedetal:y_pedetal + _frameWidth, x_pedetal:x_pedetal + _frameWidth] + + shotNoiseFactor = _cfg['shotNoiseFactor'] + _shotNoiseFrame = np.sqrt(_energyArray * shotNoiseFactor) + _calibrationNoiseFrame = np.ones_like(_energyArray) * _cfg['calibrationNoise'] + _noiseFrame = np.random.normal(0, np.sqrt(_pedestalNoiseFrame**2 + _shotNoiseFrame**2 + _calibrationNoiseFrame**2), size = _energyArray.shape) + _energyArray += _noiseFrame + + ### truncate to desired cluster width + ### odd width: highest pixel in the center + ### even width: highest pixel in the center 2x2 is selected to maximize the cluster energy + highestPixel = (np.argmax(_energyArray)//_frameWidth, np.argmax(_energyArray)%_frameWidth) + if clusterWidth%2 == 1: + pixelArray = _energyArray[highestPixel[0]-clusterWidth//2:highestPixel[0]+clusterWidth//2+1, highestPixel[1]-clusterWidth//2:highestPixel[1]+clusterWidth//2+1] + carrierArray = _carrierArray[highestPixel[0]-clusterWidth//2:highestPixel[0]+clusterWidth//2+1, highestPixel[1]-clusterWidth//2:highestPixel[1]+clusterWidth//2+1] + x_center -= (highestPixel[1]-clusterWidth//2) * pixelSize + y_center -= (highestPixel[0]-clusterWidth//2) * pixelSize + else: + clusterEnergy = 0 + for i in range(2): + for j in range(2): + _pixelArray = _energyArray[highestPixel[0]-clusterWidth//2+i:highestPixel[0]+clusterWidth//2+i, highestPixel[1]-clusterWidth//2+j:highestPixel[1]+clusterWidth//2+j] + _carrierArray = _carrierArray[highestPixel[0]-clusterWidth//2+i:highestPixel[0]+clusterWidth//2+i, highestPixel[1]-clusterWidth//2+j:highestPixel[1]+clusterWidth//2+j] + if np.sum(_pixelArray) > clusterEnergy: + clusterEnergy = np.sum(_pixelArray) + pixelArray = _pixelArray + carrierArray = _carrierArray + x_center = x_center - (highestPixel[1] - _frameWidth//2 + j) * pixelSize + y_center = y_center - (highestPixel[0] - _frameWidth//2 + i) * pixelSize + samples[i] = pixelArray + labels[i] = np.array([x_center/pixelSize, y_center/pixelSize, z0, np.sum(_energyArray)]) + sampleOutputPath = _cfg['sampleOutputPath'] + element = _cfg['element'] + np.save(f'{sampleOutputPath}/{element}_{_cfg["sensorCode"]}_{_cfg["hv"]}V_{thredIdx}_samples.npy', samples) + np.save(f'{sampleOutputPath}/{element}_{_cfg["sensorCode"]}_{_cfg["hv"]}V_{thredIdx}_labels.npy', labels) + ### save samples and labels into one npy file + np.save(f'{sampleOutputPath}/{element}_{_cfg["sensorCode"]}_{_cfg["hv"]}V_{thredIdx}_samples_labels.npy', (samples, labels)) + +def getDriftTime_allpix2_8p23(z0): + ### unit: V, cm, ns, K + sensorThickness = _cfg['sensorThickness'] # cm + depletionVoltage = _cfg['depletionVoltage'] + T = _cfg['T'] + hv = _cfg['hv'] + def getEz(z): + return (hv - depletionVoltage) / sensorThickness + 2 * depletionVoltage / sensorThickness**2 * z + v_m = 1.62e8 * T**(-0.52) # cm/s + E_c = 1.24 * T**1.68 # V/cm + u0 = v_m / E_c # cm^2/V/s + k = depletionVoltage * 2 / (sensorThickness * 1e-4)**2 ### um to cm + t = 1 / u0 * ((sensorThickness - z0)*1e-4/E_c + np.log(getEz(sensorThickness)/getEz(z0)) /k) + t *= 1e9 # convert to ns + return t + +def generateFromParameters(): + global pars_alpha, pars_beta + + with Pool(processes = _cfg['nThread']) as pool: + results = pool.map(singleProcess, range(_cfg['nThread'])) + + return \ No newline at end of file diff --git a/ggdPar_Moench040_150V_15keV.npy b/ggdPar_Moench040_150V_15keV.npy new file mode 100644 index 0000000..32ee86c Binary files /dev/null and b/ggdPar_Moench040_150V_15keV.npy differ