From 7e8d8da92fa48b284320a141503047990299f831 Mon Sep 17 00:00:00 2001
From: Zineb Aly <zaly@km3net.de>
Date: Mon, 3 Feb 2020 15:00:56 +0100
Subject: [PATCH] outsource reconstruction data in separate arrays
---
CHANGELOG.rst | 4 +-
README.rst | 84 +++++++++++++
examples/pictures/reco.png | Bin 0 -> 15443 bytes
km3io/offline.py | 156 ++++++++++++++++++++++++-
notebooks/OfflineReader_tutorial.ipynb | 92 +++++++++++----
tests/test_offline.py | 52 +++++++++
6 files changed, 363 insertions(+), 25 deletions(-)
create mode 100644 examples/pictures/reco.png
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 5a77c7a..c881c2e 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,6 @@
Unreleased changes
------------------
-
+* update of reco data from offline files
Version 0
---------
@@ -44,7 +44,7 @@ Version 0
0.3.0 / 2019-11-19
~~~~~~~~~~~~~~~~~~~
* Preliminary Jpp timeslice reader prototype
-* Updated ``AaetReader``
+* Updated ``AanetReader``
* Updated docs
0.2.1 / 2019-11-15
diff --git a/README.rst b/README.rst
index 1653a02..5c74862 100644
--- a/README.rst
+++ b/README.rst
@@ -645,6 +645,90 @@ to get a specific value from track 0 in event 0, let's say for example the likli
>>>r[0].tracks[0].lik
294.6407542676734
+to get the reconstruction parameters, first take a look at the available reconstruction keys:
+
+.. code-block:: python3
+
+ >>>r.best_reco.dtype.names
+ ['JGANDALF_BETA0_RAD',
+ 'JGANDALF_BETA1_RAD',
+ 'JGANDALF_CHI2',
+ 'JGANDALF_NUMBER_OF_HITS',
+ 'JENERGY_ENERGY',
+ 'JENERGY_CHI2',
+ 'JGANDALF_LAMBDA',
+ 'JGANDALF_NUMBER_OF_ITERATIONS',
+ 'JSTART_NPE_MIP',
+ 'JSTART_NPE_MIP_TOTAL',
+ 'JSTART_LENGTH_METRES',
+ 'JVETO_NPE',
+ 'JVETO_NUMBER_OF_HITS',
+ 'JENERGY_MUON_RANGE_METRES',
+ 'JENERGY_NOISE_LIKELIHOOD',
+ 'JENERGY_NDF',
+ 'JENERGY_NUMBER_OF_HITS']
+
+the keys above can also be accessed with a tab completion:
+
+.. image:: https://git.km3net.de/km3py/km3io/raw/master/examples/pictures/reco.png
+
+to get a numpy `recarray <https://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html>`__ of all fit data of the best reconstructed track:
+
+.. code-block:: python3
+
+ >>>r.best_reco
+
+to get an array of a parameter of interest, let's say `'JENERGY_ENERGY'`:
+
+.. code-block:: python3
+
+ >>>r.best_reco['JENERGY_ENERGY']
+ array([1141.87137899, 4708.16378575, 499.7243005 , 103.54680875,
+ 208.6103912 , 1336.52338666, 998.87632267, 1206.54345674,
+ 16.28973662])
+
+**Note**: In km3io, the best fit is defined as the track fit with the maximum reconstruction stages. When "nan" is returned, it means that the reconstruction parameter of interest is not found. for example, in the case of muon simulations: if `[1, 2]` are the reconstruction stages, then only the fit parameters corresponding to the stages `[1, 2]` are found in the Offline files, the remaining fit parameters corresponding to the stages `[3, 4, 5]` are all filled with nan.
+
+to get a numpy recarray of the fit data of tracks with specific reconstruction stages, let's say `[1, 2, 3, 4, 5]` in the case of a muon track reconstruction:
+
+.. code-block:: python3
+
+ >>>r.get_reco_fit([1, 2, 3, 4, 5])
+
+again, to get the reconstruction parameters names:
+
+.. code-block:: python3
+
+ >>>r.get_reco_fit([1, 2, 3, 4, 5]).dtype.names
+ ('JGANDALF_BETA0_RAD',
+ 'JGANDALF_BETA1_RAD',
+ 'JGANDALF_CHI2',
+ 'JGANDALF_NUMBER_OF_HITS',
+ 'JENERGY_ENERGY',
+ 'JENERGY_CHI2',
+ 'JGANDALF_LAMBDA',
+ 'JGANDALF_NUMBER_OF_ITERATIONS',
+ 'JSTART_NPE_MIP',
+ 'JSTART_NPE_MIP_TOTAL',
+ 'JSTART_LENGTH_METRES',
+ 'JVETO_NPE',
+ 'JVETO_NUMBER_OF_HITS',
+ 'JENERGY_MUON_RANGE_METRES',
+ 'JENERGY_NOISE_LIKELIHOOD',
+ 'JENERGY_NDF',
+ 'JENERGY_NUMBER_OF_HITS')
+
+to get the reconstruction data of interest, for example ['JENERGY_ENERGY']:
+
+.. code-block:: python3
+
+ >>>r.get_reco_fit([1, 2, 3, 4, 5])['JENERGY_ENERGY']
+ array([1141.87137899, 4708.16378575, 499.7243005 , 103.54680875,
+ 208.6103912 , 1336.52338666, 998.87632267, 1206.54345674,
+ 16.28973662])
+
+**Note**: When the reconstruction stages of interest are not found in all your data file, an error is raised.
+
reading mc hits data
""""""""""""""""""""
diff --git a/examples/pictures/reco.png b/examples/pictures/reco.png
new file mode 100644
index 0000000000000000000000000000000000000000..ce5fa3c6ebef1d541d93ef5a8cf563e024c9d9b7
GIT binary patch
literal 15443
zcmbVz1yELP)a?OOkP-nwN>Wfllt#Ki>5%S{?gl|Ak#1>FLONd%Nd=^&MM*)VrKF^}
zoA3MoTl3$0XYQF{9O31h=RD8eYp=ETCQMmT1`nGY8$l2}Iax_n1VQtGuj^Ro@SPi+
zF-GvuRYzHE7X-m?zxYB+WWlFE5E?{IQcT^`V13$4k4T0BbBpOGRvk_D%U9$Bo&irT
zvlA#6e<f)wDkU}e`u!UY@k-8<N)r_uDI2y9?>D6+HXW)RIu1js>dHFA(XUXmkyvQx
z3|{9q$wC!)F15?XEggqu^{ts{Gp-FxD{w@M7OJC*W32?mOa>5$v*R)1i_3c4cOnq~
z5EvFnFia&+85FZ&G)EKU+L<t^EZ121axw@$7sBg=7|Dt3NPiQz;e38K^D1r7>Wln#
zD>F7*dnQ{N!q>*|D;M9X<RRiBxJt=-$%1i81B+H2CzY^M9m_)eJ|~)sV5b<?=;bej
z|2&~l68(9=WKfJ#Te}$63jROqFff@*Nb0`)*1z%DcSJ;v&vkk5mXLSEEuqu~ahuz>
zZ~yr5<JPTPIXO8~Zdhbx$*=j`-ieo$IQ&F8{`@55vtPkh;ODn9ka6sByubVR$J@xN
z#>Sa;6Qs-Y{B$yncf|ei<6gA_#YhU?s9G8`T&^a8sG^Sk>~|%P?C9;5O3W#)7fI}`
zuZNQ}^?xl?o2+y2Ki*%Dh`>#eCm|vlY4C8cXZpO6;NSAlRliVOanAu|D{q;jj36?c
zoBlPww@0TwS>Um0Dr{jbFQOxlZu{2!9^NP7ziqz#d7tDrqleG*+Xsl}GuL@16*2_(
zrJ%6|-SXTwBSX)bB9Wk8i%&`q;v4)tcC)5BZzGOYd%xYESe#EK_m4aOvE8R<iCsA^
zuQMR{`n-|ka>qiJum554j@RKZM>NU)pX#Qm^b+RzS9v05sJ(+ylaxk7%!2ZQrp1ZD
zID6`N^c5Vo^G9wZ^4yA)$idH0Y7<japZzs$b#;csNA<<U5B(N?qTJ9tvJ#G1g}>z#
zPG8C*@t^T}=-N?8Roxik!y^PA(9O!~*E;nsmi{JlZ{~PoS%j?Ez{BpT*Dre&tNItV
zv~C@(*!9C}(?K(&e;t0AsJZPl^4x0v<Kj7E5tUT8@mR`&G0%Mac`TE(C->=tl7VTa
zJdETFN0%z^x-1GrUh7;&J9utw8%?oxk?HIi4NEdhVvj@g)RYs)(g9%$BR&c3Os>KS
z^Iq57*;JN_#^hz3l6Dgp85x_(*F32clOB4N=tycbh2_TDNh7cL-c<^Dc5XlGcs`6P
zxXa_!Ovu>9rA~ylGh^*wYCQ3|!y>(I@`sBD<3Y_!0cRIk&199PR5a1^sD5FpjKAk8
z0v=Cwb3|!aG*8@?tRyIQzUeUGWNAxGaJkf63P_49K{Gnfs=qjiXcyk=*gAZ=3e<Og
zWNdFXoKBtj=sDA{{hj!LkosIlaF>T)3i>Jaw5XkS%zgfbR5+S-M&*RwB2DDiN*h`)
z&iuuNjfEgCmX^glcFW8;!NI?mhaTk(=O|Hcv+m!W&q+GRp;c!p!XrMy)V3SYd8aT}
z^K|INWP5p7pwz;Pg%y@-Xf#X4f4WOTsTOb1txM(IC`VOc2B7JaOGU14uVl7x*bA+U
z=Nqs2OiwmUJh9LJ%$Jh5Bg&kp;fnp3@tXO9vy&-Bo~qaHW@1E#yUIcUL-a>{9OFJ6
z$L_s-I>dbA7a@-m!LXY3d8iu2z1)Hc4b+pg_6H9NKk^FfY(Ey@sm2IsJc(n(PeKMi
zbz7^KTU2hkExBKQ&7a%%h1TxXBbI!>c#>v%`Uh{Qt(aP>K0Kc@XD;E8XG<@gPA@4d
zp%T9}Eavdj_A&zDUh7_IY#7G4Pf~xFXU5OoLX7A!W4u&mPP8`UkFq{~T2B?U8e{#k
z!apb-<(4Mkw!h(NZ(o*`MZJV)ifH-GTB}%f>@`RUC{5p*Zs14c|2Uuo7G@0D=QX+#
zArk{hHw}7O&;xdLl~l`wi(**%6oM@nUtD=Ipmaw*2%Xd?Jp6oqBq`N<`mIyFUcJdo
zqKUSYSMc0P--&6;U$ks1iC$g<o9&PJJN40zndf_m{6*e7)nOsh-aOO3yPqQ)q&liz
z4eXySOb!+1MDlyQ#adz7KBs)_#y0usyA}?Lzf+~Sq=c1;sUSN$U}S^d$zI?%TeUD<
z_}$)X;XbtNaNWSBog>zX41W7*YVn-bYK68vwGU{bVorxsr7OG}Bw=5InkM-_H-8VD
z?l1bY@$*lzcwL@{omJLuEN{ajiljxx7Pj4#I-#b8W6j6C3+wuJB7O}-$eqL<Do!Dh
z!!ze{r6|GJq=ts^_1B}13)P!>-??qX6Smy$eAPSeH^Du%-OMB~Wu~1*D(XVLUKgmG
zVXCCNOigRYz*e**!0uT^KYQ##>cBcd_VjKo=|e4ST+>?Yeg~bmq3KCii5REqT{`>v
z>?UjN<Vo+jxa?F_aoDpUNJ2saE(K3XP0iQu-#u1_^QK0zn&gWLKUx(IwUIwaY9AuF
zGr_upddf%`X2wBQcOysh5jXo7iCvNgwRL%aS&^%4cV-J8TXP2ap}Xq&bJFJbDrxLR
zI0@E|QPzKmgA8swpB=x8b(OfEBaI_HbYf((ooqU6GM8#pbjU(%UBd3`++C~{!ZeN|
zoYtVnrD`F0<za7y_Frp85>{m&(rSzx+_`tJY(<U08^u@ume$ppR5T=uFTpyh;_wZR
z%~<1L3>Enk+TEh9j^`7(*<x+NB|NUyVrp~hce3@+7Mzur*U>b%If`l=Ur4)&p*TY)
zK2Nj|@4do*hQ}SK*=aEQ{{9_nj2Ap{myg)3NU&sau9A{J7$AAQKFa6Ebvvp?Pn$*!
zU34^hvzeHEh3TDi$V9fv9bzO~{SFrcz9mjbX^6iRYI|Eog|sEuz*sGW-AR;&06Q2F
zJu6-7?oM42QyUTcmcaKn4;%5cWgHZnYs9tXde<G>Q=-D0=y-f>f{6AK@`=RYMb9p6
z{inYCM{CX~uG{@_Bda1!`RtiY{X5G!+HWr*%6AeTUA}W8ZTjQPOiR@@@sQb}w=z?V
zrRHX4@(-0Sp#vOVB78{mPB5$ROLCE=oX(s=dK#yJYx7IZC>N^sLroRAd&vcP50oUZ
z(DYJGtV~SkN^-qlwCQP$O3f_`;MM6T;MVIIsHl8g2_T^B3>KxKk-^rs5V&Krw~=_}
z?=L1`*0tPnNvQMZzI4}e_ky!?z~HaTMYm*7rmSHSL$vhTflGW{W5hRZ*p1&ZR99Cp
zPk8kT+c_X+@Hr_CdCyIv(BQDZh({_aDo+n&(&{=wL;Q33U1jGSqe$%wQaQi0=SWD9
zjQMw?iw#h)WEDtzdcJ=@JW`+}M{h^_sjUr-hL%>0m8A3Y=Qv@w2mAmtz-6i5>4=Ge
z0cmZ0j?pclq%@#iZLzShAaopyba!_<G#}<m(*4=rH)!(aM%q7rJ~}=gsb_a~aS?D^
zbv-$7IwIC^uJ@h(;0Dx%>Y?k_-@i(a9-#+pZEcNLTaE}zgy4pShhqf5eL|+&uG34$
z;ESA}xF{<tFZE|MH8u(#Ex-NrN&HAJ*W*?0aK%fdTMJ9&`g>j<M+&N!qy+tXj$aJ>
z;O~0UuCEQXTYU{NX2hrNZ@tUnu%N!KA*HUa{^ZHGYa%DDR?*?%y+yEbQsh!}D!#sr
z6BH6;WMnnxTR+?6OeG~H-QD+pH2F{q`^uW&AkWTE|JJ)KUBSdG$j{#xt8AN^qD1oL
zQ%vmbo9nAS|5(cM|Iy;_Z)j+^tz7!<-8*Xk6Q{mZ?xle&!)lAj!oqxytK#C~RaI4s
zUs6WuT`c$i_VxCrJuv<3&~jP}*bq+2KBq41ZFd%mLoMuG^Q^CDqulM(dGr<2<;&Ff
zP7)>y3HiCyucZ_<S<8<Y=KO?^k6@H@vKI+SNg+E;Hn+BBCJf_?R8&-C)2ZR}J#Z8i
zMJyBR>+7xkbCjQ2TTe84?Wnw!X488g!&PS3@@&2<7M{$5j7?0$sS$E=R(ExES&igJ
z-h4pE#I$tq=T~?4fJK<*??1mf|NQy$;K763Rq5|#`t?FZ9EA)7SYxN@4@|FUKON5P
zetDhVeSJboleX#CNO^R*hlu8PMr82o=w30qCbfvyD_D6`paeb%0_oY4+Q^od<dpHt
zF$VGeGiITX%4PNYTj|u3$ex~_C@O((&3-<}9X7VUFJIJkT?a=-Tn4knSe}$<eQ;e#
z=XZfme_H<tE4{<czByV<>d}D%|8ey2SgNHvec=9OlEEJ9JL5gm&c?oy@;P%#DkH(L
zHZs1436JvP<~<bpGKCs^Ztwd)^Y=gZ`MY6mVZq^th8*p#NE{nBc_&`%%a4(fk?CoD
z#(2Nez3&K3;-eeb!GVDo$jSbsL*_r1|Jfna%K4v8YH9&@`5fg-n3$UOT1^m;u+0t*
z1Nvx=2V2G;n0W1{X3G|1U%&P|J2`}Ki)lb7yW{)<WA$FX0Yivqd6BcT^LvG(z$h%s
zt}&Lo00B>)FvfS+WZg^<%F_GbT9}+ns$E-M{dQeE$MI9pRj3~Hg<_|>BTDm8xCkFV
ze{AqwwfqPwf&BrIbNKvSuV1YnKZ-FB#$H|Z^Yx93i}UsMUB)w&k?Fjr|3Q|4u#eXq
z$#(kn89Vq*c6PQ{;FWX@VPWB2_b)+LiHge9H@jl#A;IA>E^clMWfE6bTh7mj!X)Wp
zyQZen^_=3BZfTx|SVoF?m8@L*uI*6#ttzZ71V54@SX{;6pG)WBO3gA<&?hfcw}#8y
zTuPUclJcO;MmqY|`tRSp^`*b(=H7TycTP-vfLN3W#fNj&t1{~@({Bh8H!e`jP*qdI
zKu4FPQz<ok0YMce`3ArS@0Q@}>Yh)ZUTk=wx-TOwr~8v{q{3@$$HmKzV`3d09V`C1
z?sHyTMA*Sxu4C^Dt>$iSYuxCgsts11eC`HG2TXTrYU<6KH-5+aatwsETt-q-m-$8F
zy3NvG;MLaF!fJ#(w`eA12MWwMQdXG)AYEbST<6-aBCkP6Wy06y3?QuQ>5*F>Z|xx;
zA0HgPjFePfUS48Cf|-Sdb@K>oa<4fDnHNjKX=Ge{yvX5v{MeXQrB=yN8eO1FH1U3`
z%6fLw>UlV)?Aq9$=}spnr?AM#>Gw{iQc_Z(_*Bcfvy!s1So-_h+ec@|n<<F|l$4Dw
zONzA}ZEb<~22waJ`_uWqJbF4-@$<dYuP5aO>%T58A}T5fS3?;|NiO>-J_d$))OQ0@
z0UIj_J-XPZ8&mbysD)kK-0*O4e(f8zcXX8M)yJ_uW{87*KZA%u2q-9yw|+_-xi`x)
zKq9^w5I@SR0i0yy<mF}BY)VQ_PEJfrNegFXeQ!C_N=a%RnE%bQxO0Af{_9b~St5N4
z@si=~9VzC23Yl10(W6I?u(7cnYq8CY#>Sb6v4ho_QI(Y{P5V%L;YGiF`!;_CFE%!I
zC|eBU@#DwN8USC992N#~SY{8I>n%cjHS+M9Tcx^;?hwC<izC$4)kWpd5MbHb**$*r
zNL^{)@bi$dtEL9KDS=_4fXBw#bVKq8pYwc|MPI5DdF+GuM|p6dX8a=O-Tw7n^N}jL
z6chkT`1t#Xk)55smX>Top$1cD)X;7GVnyk#8`j@S$d;ZFh*Rd}aZsdLR{x%EKHcny
zcwBehB+;yBZfQCB#%rfce|kEV7XzW7NE=oqTgp_-pr*e9YFk(uYo#109U<M}p;A<q
zpdZT;kA7m-IfKZ#oY=<^-cWoe#|l!cOEoTRcTmUZL6tL|VqFzsk^{mXsqSZV6ue(W
z`1zOE^W)8TdivR;|1jcbV$+(6wgmYuOrVNvG8?@2#+h9b9atswsqcPH<HI7uILh=e
zURYBs$|NR^xxA0Jojjj+jsz+DDXV-u{<$~dSJveIU%D5VL<QtWy1hB81HQHK^WWu^
z{Pz#T0Swr5YY|a-dEA2Y3;~Y-f2JGQ9U*uWfFBHma-yP8<S`MWo`kzQi+%J|$n0Q_
zgnom&-_|Ueofy`?cSJkza0cYi5I{#fJv|5ewY9Ys(s*vKTYGHp42aC2Y`*s-vWyH5
zzj^ysj)YyiYJT=(ASt`SB?MfA8#iu{kmR>;5@aW%!QCt@Ea+limM39|eIM4r*EmL3
z|K?1y{KdeQEp%sTAYNLh_v5#3wdTD^@Qw18v!k=qV?g93C{v_tI+DEoGV6**hld#t
z-BzJ?YaG82rOD2aqPu?D{$W*7*lXJi+zaYw2gc<FO(mtJ)6>(PkrbGSy}kYU@yvPB
zJ^dneX22I+Ow3TYW3>o6t$g~!CrZ}*Pa$xGY{%~R2N{o=cx(^nJvKKt_dDAC>5hId
zbAG0&rS<;V>CyZ5?_p&tvqa{yEwj{;Y@<6m_?KoB71Q~g{%%aC-<1S_va;OwR&ZmI
zGgqvN9~v4mW9mlb*4H1x9ihAmz#~|>1J%i9Jvt3<?H5<>H6~$vtc!;QZftCXg@x5z
zMkeYU&e!W^;>`D_B;A-Oc<pf!rA$F9fP~}Cmd%wFez#TiqxCdcPRi@o>ORp|$(RH?
zhp-LVFh5O9q^{t+m5CPi-n}ERM32CF2hVnsuuB*2L-yQ$oP;2y3kdRt^xcB!(=oZx
zVi8#IM2Z=9&WCv5D5#nm21^4-E$d<v+FTMweD<cYiw|};`yD?~SJyLbRw_|vwz5o-
z&wq${d%F7y9=Uuf*Ij1j)sw>=K*@>ZaS7Stc$Z(fn<&cc4;ip6bZqYU2(T`~upr?n
zlnrUY%li%S#>U2ri;FIcU%YmgRf?aKLfD`h2PLHY5gY!&nU}1C{)zHlg7td+PrBN4
z@VmAg%YM@*%T!<EC-EozD9nJO9K-g?jls}mxTapjH%l3Sm2)IDl2VcROjv)1AD>NI
zYE>iFhs67+@_NZafI;TTiSN<4ifZ{Q;$rWQE|fpot%lxQu9H)Gsz|wCR8&$Dp{Ft3
z<dZ!8qPbNEe|f>l&F!|PMs5g`fd*~R?zR6k3pud!^9Mu}KIil*@qkxoM$Dytfcj;q
zu>`5efz6i6yTqem{o`a@PxfIzSYV8bqqDPZ-Ui<C6E<Bw&1WbrOjnQ75D;5Bb@rp#
z&!^Y@HwD>zu@U!L_r5gVQS&JT9wgn--96OZ9XnYh;i(Y}&g2i5^z?K}`}oYklaybn
z&w6PE1=F9!(+CM=e6hJr5<V_WFB7GQvepn6|McQA4jmnx?A1^pr|?Tq4xzM0QVArj
zw>U5KKsm;v<l9>8OG9P>))yA`fZ^~Krpt2i^Lw2h?E#&E1p}z3iugaXwO!koZn#GE
zP%?k7r>~ES>**li#Q|e;Dv#X+$Rq#5g=8Re%!!Yvsk6?{j@bnSd^%p;fO}+f{@Ps`
zxhP_$pCh0c6l<3MsJ4t&XI=+o0o44hbR^O`IjPH>n4g=g1wfrQ3^h2$u1X~Qq2nxa
zh+8Ae073%i>QypsYhrz$)vtxKv$K_^ozM2iEpvk-)|Es~zoKk35iD$Mavs~U`R;hJ
zX{-xyg9~?>21yt6irZwN@})Pjwk|r_R!HZ2aOciTI|d-;A^229`{r|#DgNiD9*5ft
zg@uLA^^FkrKY#voRpmO3OG_&rao|S+RQx>Br_os;tg8ovd%L<kmohxY_92%b<zmgr
zwdi8+2n#=(sIgj;e<&#}eI)}U=Gh4?9i0iR*2~KahlCaD+O_ED=+GV=d?<pTcOcHM
zrxZE>#sWM|0#iRf`ztb7mLvh(`0wV-47}jiuj)`KH&%<>H>M7MMT^+~XsnluXN0vu
zwwPHW*}1v%pg|T-j(kf>N`^;9dVx?=@)MGfKq@V+tjImN&BVlOH?a%q1|GWh5?(Bl
z0YA_E5-jFgZP^tnRTNhQvb4a-&E412gN|&AgV2fl8wcRI)hORu%}Gm33j|muit+(B
zw~+st;Hk}50epP?gXuq-t?~{w{EcTvzkxG0i06#eIb=<ow`#-7_#Cg-HMp+K&Ce5)
zkR-*$iMQC0+`Rd&b^@20L)*m2D7(1$cQ6jA>|61i3(*2c4zHu5L%v>8TM;yIo#v9@
zWu~vg_tDdnqwSL#6Lvf3GDe%_v+XUncXq0)s{;cADR^w3T3RkGEp?`tVj0hUeK0yM
zSiBc!&99}S<GwPy03iUn@%K!Nh=-fT_i{sr=?~o6jHFbq7A-*OQd3g$a&rg!CuEhB
zT0Ayy9p)i4Y;2|3*+#Ri&w<v0o~=o<Hol5SudK#R{@l8forv1U=lDsf4x3@~3B<B$
zVL}5rJeJ%r<<j|)dN1vrJ5Pt+NI}Y|3CaA5d6Y*&LSmMp7f^xur9F%Y6)^;S1{AaX
z^aq8Ehuj<-!vh1+HB3xQ9-BX=ffo%7u=De$G<};TP8t+M_r1(?uoa2S^cEzo_Da6n
zh!CZ@u&Q-X)+o<EbC82dS;9iRyvh4A*@wdTb4CL{4|^D}eyg#CLMWoZ5N+}+_9hdz
zFVxc3PE1OQBdKcbZ@+c(=6mbWBCzsDNQ7Z^a&mIO$$FF69xEwPb6Vh9DJd#u0x}~9
zhlgzXA9kTKxspFoG``wL->T{OB6<~zz?VHsr-+>>$K5hLBLnI#c&~378yh>qNr8f{
zja3E;U)DHJr5&8q!Ozdn2mQF_r_I8`;$E~5-o*A}l>KySBRxjqceOQ&hY+1mk$yD$
z3E2IlC%}Rm>cZ>60&{b70in=DrKP2z`lWK)Fc5YFZ@7$(=3{(WgI$U{I)F^+<5l4z
zdeNtzhw6pR@<EFkRv<_lHk*(C&iK1bL{dVD!g#lk-P7Cq8(sr4x}&25s_5b2VQ!SZ
zk<lzf00$>0wUBoWsErH1qR>0*>B+$#AQmi>Z&KdA{RK)`vrO;mg|Pu%#e3%k?82wQ
zP?|t4RA%r8-{y)4sQU1NNN0F!+p4wsB2Pfc19fkKY(0m64h}rsrTs3w<|7Ay013Mo
z$@^n@v{+M7G0K_^Q||U?1M|m^A30JY?EY+RWuNlej<uGd=gXrX&da*Gy4u>>M&A<3
zSIjt~Ls6>ltxxiDak&DkpKkD&0Rf<C=r`Kv<$_@0<KGhY^#m!@^6WI>(E@B1pOWp=
z`%0d_3R=+JY4YA%{n-{8A0H1I<I&`d_wQ2OI(z+k=in6_A|}P$WH)l@0OIuQu0mzD
zr_T-#D%gp(*5tLd$05KQ+}9sMX6ourg3wHFxbZ6lMVaY$tCQM>pQi0oeVG#fG3Ig1
zTcu5AYPEJ9d@Ndt>`~oY@1A*UpeXsWHZT!(_HxZVl{#+}=R2L(eBNssqt2eUMG};n
zCdky;8jGOs_HU!>|3FzWX>$c$^I_3WfBd&6rpOJpR7YZ_u-_wcE1c`MZs8+WuU<8&
z<m2H<+6{a4icX5EKSQ8%B`YQdK`=1!8LkNl30d}~CMGBUU6Vm)U}9=HTudVjYoDC-
z*;^gG&_2I6H<gw&y88OozZQ-F?^{vhJ8%y<dAL4-+S%WaiilukWo4j(qSdyMqs96}
zK$laT{et;|UIwPSmT)Smr$?Exew&3QT+xh{Fl_GEuhIUe@v*UD#~{4Yxuhi}%`7c9
zryKslAu!ZtSH+3JQCVBBLXqtH@?0%n{sAxVOrsa)%3B3C%;pTThpx*2ir>l&{h<zl
zD^pnvQM#gtPstafAFQ0w+FlrP_RV>AlI*#tk*Vnj@N2i<-+*9MCf*p!_-di36l9Hk
zUbL!pTDX9WuFJ<WDv)zqyB=<vZ(s)*kjtm*%&}aVHO#W+-&)5tH#L=TP(%v=zhS)X
zV<0AH6G0ItpRN?%vi`+da&lGgHtayo%$qzy<=+u^6Ca`82*kVw-i7E+$=jfPx6Z50
zI{5O*OaHUTrz<M?BRfQgmXeb4IEi%}h+U}}0RaJxYN44jj<rWfS8s2F$EHE@J^iPT
ze|?Tbec>?M0z~ho65^cxO6Bn1`Imd<FOXIo0}o0)dmp+i%3D}0>N+%mBJ0%bea?z=
z3B4BKDO`w&iJ4trSD+YwfyT(ln2sU9z|h{=xdeE7di*yFbbACD=LZO5qt*a;YEkpJ
zI(?DATXp8}SFh|LNqi5tp`*IHJoI*zlK>rY1b(+enU<FJOR)2U(S3+wf%t@kuaGey
zzui~A!g(mj%LExb{Go$?RP4UQcXoERCQUHx!LHT&v7_Vi05pz|BL52G7o8huuwe#S
z;}OY3Hx)r@jFMkaMpsOt|2?MvkNxHg#dhUMFB$?zKl}UlWR65A?~HG}V%bR_?@Wac
zpGiB`uz{*7AsN{?ycvWO{w+a1Vd2c0mg9{ETYq;0d=%JQ)!-l~{sMF!k9^r-vQ6dT
z;NUQExH?)4TQdc~i?p}5J08u{yI>;X5)$3EgOotW*L^RX4_R5+t=-*RKjG6|wTDfO
z6R{Vh?Mdk>IJKQoR0RrY?&oKI!yD`CY$V~g|51xcY`VCc6hrJ+85Oy4OOo0E>z4Y`
zNTK?tre2&tObmglC?f*{0}G1{q!cJl3;!8_7C}zV9l%Qw|Fe@KjbadIEBRe&&)*;d
z2tsjqF$?f=#{WE{)~4L%dl_UMLq~s>NLFHEM@r2aya^Gb{A(~r;K=2|Lqjo7bN-c}
z=Ep8E4|pFu7zc;t$&;i-t$B(Rc;eZ?ujn*h2hfALdiSp3;Jj4|i--{PJ8}j7p{F;M
zZ20WZSNA3Lzi%L6V)CY;!3VfXM#glBb~Q0}2XI`p*7o+Gwb#&e2$OVuQTg@jR}j+h
z?%;uaXlNKI$yK5qyV_mTdB(~iyrx6QTWvXTOXQiL9fML1mr}H7ot%MWAE(fE<QLFR
zKb{@01I-0fFE}{(f{q+cUP6Bk#RqIU;OyXZeNy}oa`4m}58(M;R(8@o8-Be8_fklD
zv+=KAYox+SfY?|HIJ*A+W(O+0<zzcNsvx_szaJiX0(D~Sy;tFZ^y1Op?;3qp5{y>p
z_A_VgMudijsxyxY#=QhUxh+(>psc6&CG!#N0Wm4*NN7<Kt2T$lm)D@=q40Q&RT$ek
z6+;#b&piFn==CLjRTLb8YKy+8TS7Vp2K}{`78b+QMTc)x!RWd*6uCBD9m=&G83EQ$
zMFmY{Qc_ZGusb~zG$f$()$@Z{+Z)`inQgzV*f60|x?fREk^ihnTHZ4?b<*_rhtFoy
zUu+zliztS4%iwnjnd$?^{}mEhme*VUTOt`Jtn7j5@?eh5x00}zFSFOzC+p6^r0uR+
zEM~w1?*IDrYuO~&iQdj=YK-%)|6oAURjoh1&~@2#o{`ZAqkQ`!@IV*8*kR&>{0AE1
zyy)xd<dk!|;zzC1YYZMapcfiqY|e&>5R9s@VY*uat#{d3RHy=@518a=zCzkk)8Y8U
zL`YB&1W-bw&?SV2hi8{}u2@mc#!OwaRA+D6V<tn$C#g0qC1r?u0wQsdj2qmnW$7*}
zoULei;P_So^_g_vpOnB{g^i63o6K;R@!Zx0-obW(^jA^2ez8iI-k1a>z6+01>&X-8
zV@-ASieiKY+&-W-^Xrq4lAG{#es&6U4&p)cxIcNHPhe*;&B4OVjFQhWyP{&NxV#x1
z5fKqVVtN-iQsnQ!M8Ttd-iZoh3GnQIdW(sD{2WQKzP@f1du^Ne6_<7EKGUS7KM5Ii
zznb^{8o$@}{Nc`$h}U*dO-z9@cJPZD7Gmt$j7Cj`I=P-pL;co&BIW-}Uey?wug**@
z^2`@Db$NNYJD%zL#(sU}$M&D*F(^u~Ylx|+E@%oMF(ptLz}R^mntS&Y>0$w#z*~L9
z5a;p3D|EN6vQoR&RvY}AqraQ}KR*Qn*^4p|jV(~lRS5(KW3V5*ALwQ8x`s)LZ{Pu~
z)~Abgyn=^-w6(OcIR;r99}0$pD?MObIxc!Em`FXf@w{hUuX`7^A*^H>;(##?^%c!g
zq=2@Bij9eh`S<sCr9$;G?dlL=^FMPPU?HE*g|mZ040xI4fb$1z2(4Z;GRBgH00Ito
zou-B7s|DT~y*UQ!xJ7mMop!aFiWZd1-LH)5{;S*n$3)TpyBMTf-_&9Y1z!o8C1LKr
zM+ydj__G@}vyy~cRfYnDzwDS}Z$3arAQq+Fho`6gfPwp;?vbv%txakJ^R7!GOp?QM
z%V^uiI=|NU6I8LW$nXvml$s2Y$Hvt4@D5kTL8#~;EK^t?e;svRo~q}Z>?(T|9Zh}T
z2qAMnvU5|1mnbqeLQx(-T7f?hWA*ge*`cMobD_tKAFqJG@s|`%@XCRR<rwy9>h^x@
z>0yTNg6%`?yY(@Y`kAa{ymENCFM4cBid>zD<E-eS*MzEWFPtBaSsGpJGQ=fa?~flp
zFt6Q;0J0GkrKqF?vKKRB^5!U;{6AZOiHDAk4)smeg2t3chKAsEYU;$q#FZ;o!01a>
ztn3FyAfI&e=FR2RRqvmlE>D!bg6DmV(pqP{siE8Z+@9n$DRhv}AOM=6`7`u3RF=Vd
zylMfan)r1k+}`l8GYpHOe8o|&+Pz|d0Or6Qm>3!1{?6agkk(jw89hBcz?Clj97S*M
zBe+Kwdq-%w?<bZ&6E@{k0`9hu$mb_J*F30-1V9w)qzd`)Wc=_2tTgS6y4?#bA<J(+
zWkLtv)eD@yd1bwpXMRv@!1!otY8w6Kd@vucgc!Aj;K7kz3Wydt^K1=3D=RB2s>Uzs
zsN=Lo*^a~PRyuKz>43tsqm_cT&}*cmq)>34f%m^#_v3@xZwQR;s`7I7`b{-e)wz`w
zZ3BaQ#Mt!o^xf?lzOTZ=_4V|YB|ynNbZUx^S95meu+y}Ge4PIht9C5r^Zfa9a9>4E
zw%VYeQc+R@!ipKHsfXTovTjAJ`QP;qe>SLypcWD4=1x>VLp*;)@<Me2+=znjffov7
zJLb`AYa4#3APl!}gRu+7GOLh~A0V}>K+!=D69F8p+058V-XXMwKh)QwwF28tW;g8D
zFU<(}=i1)h1{ekcUk4sDkg{2KarP_aMKvbC!@$=V^XZ1$uCa6b9&FuUc`_Ip1Qio}
zfV8wpDE+0SrO{LZ(i5<m1X!2)`}<3EP@t7HR8@zeb(`rRd)e^t(i`Z3nf;0+=YA9E
z3qF{vrH_x##s_jB*c5`Ej-YnV&W;^69`1g+@8n5l#)b*<<#F=8y!?C;Qc{SA#igZ&
zh6YW)E6J*bCGXyWzz5Ic7@|R^)>dwU`@6HFt1Bf`T`;`)*IAplue#+)qbLc9h``4=
z2l+ZYJPf&JVQC2?AHkDXo9@d~8eGxt&2gA%h8aGvxHo?^9Rmcn9IX}uWNEJV@16dg
zsRPnuIgrT$QuMu}F^D`_9~6obbRi(72A)IGpG2j1EZjv4v&+lx;iJ&qYucZ%fq>L)
z_Vuc>pMIY)i1S}y6jRfkEM-lBl$e+;AXk-aQ0l3`SDSpXDk4uQlux=qv}(=tKo1w=
z|A)6iAP9PGGtkk!0X_iWV9Y{^v=*urw43~athd4$7x(sl=yj8rxW>Gf1lii(CkgK$
z3}c|7F#<0lz5HtwId>5-95Bw`vccAv^_xZk?-~&z4(|Z*z3@*XE7l%F`ArTpS>V7{
z(P3f6ymTbtn^W}?0@czHH^nklVVo~kmd_<~2C^f=<43G4LvwR;`r$+Mi$=KZNL!mJ
z1gDji)tb)e^+7=sS$3m2bqQ2bpx(VK4)Y%16i`4i;^H{$CTfoVZjPY927pZx@;$gO
zKoKs>Fi~qq;U_C4Ee&$?A*jTA_cnUjTXNT7(cq?dxVeqqAH}fw96^@LC<bE(!1dl2
z7ZZkejE|2){{t5f56I}f44pCP?%4w&E(Wp#YLK=UDlxnmO|R=6;Pi{LvoRpuGcpu-
zAvQM~K76<UF*DGZWNS-H$r%}b;8;M)D5dd4f-0myYrZaVgbvD~2t}!wDF_;@9vWj1
z3zHKQM=*=wl?xR1XIp67r%%AK#?LEHPfipw9(J60C}3SehlVkPkCdb&z)f)fn%lwF
zPw+-yE?^ECKN?v%naXz<xUWiT^fGUO5eO{-N<q(x#KapJD~dD(gOM*K>7I*<7IOnS
z52IX<!EYcJ^xU#3Z<6-lU}5<UX#u7Vs=68n8{2xab{*O`-OR$d%0YpFZ3%bPE>wN5
z;vcZQUefwK4^uORqX7`HeN&bo4$K<Z33Q`i24t+e6_Py^kD?Z&O}@hon1BT%F4La!
zq3g?pw~rJQU`w1vo81;+_DZ!-eceoFl$ni9$Z1Xz=n=f`vV-^f1UC*&IMmECpc250
zLF}6<`Exy1RIIJ3$sGxVW>)&Y`}i4lX{Dwa9P2;3yS~6oyV{@x0@lKQ#{G(E8oKW^
z*A9-%&S{pdsVN3BFFf?CGg?he&CkzINKkNXeLb3z|I>Qu92ir-4dmRcv7*4V{|^%o
zMjF5_1T*&bIvkWMeX7kSS}U-l8`x5S1h7<areRb7$$tNS?DOaRiV8>gHq?*m%1UD!
z&Uet+CWAY{>BB${oLI-u(DfP$H23zHTOZpC(KtG^L0-$1UPi`11^{7}z55uERa9i`
zEXv9<1Y0fdadM^*QEOjc-}JQ5g9ov-6|%IL97~H(?f`B-cx*!bfbxW@ukZ$!`&l6|
z0YR(2z6OYlco;E=i{k+SZQORn`)w|8B3JqvR@^K|&b^IkT?Rsw!_1G-(Zsdf#6%4h
zmDr}HQ=q5N-hcJ#)tJ@f=qhA0ATtCwlo04qOOT_;{$2Q?Q_ZJu&M?+t2+##+37#Oh
zys|j_W@cvdi;EXDY-)%Y>k=d*9K0(%_!~NvCczUrcYxS|iJBYt5ri@r7b)L#bCj_z
zy@2Tq=pk>o+=L%mmv4V1#177t3@R?RfoLEiA~N(pbwgl$V8od|xA~Ei(}Q}t;&<<0
zn;;y$p|1%71>ji$yB<WB9Tep$Ce7P~VX!K=L|`OiaBvU};|g>^MANPt$>6SB7wr7m
zq1pR!a&i)ORZ#E*au1A0cp17=NMr%mWygh{8{3)vz?%5a`TqjUym>P>fwHNQDpVJK
zc5(o-7|<GwXHw$f<z;7Om6wwnuQ0|K?*dZ@`WR2`^n%4dCMmt!n{3M^=tyMIgs}-L
zyDNFB-S&DC3xNp<#D(fW*0rnRVe!zQ0JkMxmH_}otIj?{o%wVj*|4g*`mo2RXnSwZ
zZnQ|l*0yBvx-wN|QJUBu6?mq2D9Yx;MGBb6gz>B0)zRP3m4LPd5->k+wqJeUp=l4^
zsIE+}UR)P3nrjaO>uLArTt_$<{>L4phCy(JgrK*=0YhM|*b#%M@%Q&n)dO4gcXiKI
zLKrK_$?LJ1hXMy_=i=d!#bcM+(sB;!TGj_TCJj!ioR3Q{0%>(Le98qa0IXSl0fG4A
zPB^MPD7ceA5;VYv&~JFcN}}nNFAX|$0w5j$77W?dm$auZP?Qx(7##9I%W2Klw)!j4
zii)8azX3+>y>Q<NS6LrQCefa?<seI_CJg2Or@i=pv10xw*QW3a0PUb&uQ9x_n3xzG
zVAg$O^vN=kZ_f9{U-OMZqfnc>aAEMB%3H8}ajAvVYoCH|8i5CvLPb)b>~9@><i1%K
z0qn{;6UGfAHca|JKs#8hr!{{gR>^`91F$Y!8-7HIEmQ~Y_W24?!Idsz%YQ%d6&;?#
z!zU8;1l}<mn&ov_B$D7hb<OhQvzII=FYx1nC#PI*9xN{AjsntnY<3HNENe8k^$D8l
zS%0T!&vfuXIYT|X@Iid`OC;%1_+2a#Ev<uKnP}gf;hA0cjlsT(scaMIyz<Tfwo6kW
z(12f_5<|Vbhv{o^hFCe}{|_(t$Ik+@kCatw5O~vtX+^<v;o%?`CLWU%KzF4UYqRem
z2jaWa$Gi5~iZ&Et9RKEhQM%aO%^49Mo?7QgEDel*y|J7(Z%nr5yP>fL>cju!PbVg|
zZ!7P0=IexnafylC&}#&6tp25}rFHYh4L5MZARqRAe+Rf!_q!#&YCNaj)@sd9Dd5fu
z<td!K`RQm;Vrl6*ST9}bKEAU2m;ex}ARm0ktIVMx6crT(!r(8==fGq`vG6t8Jvb=&
zavT2KygcY@T+B61-9Tx8*_4>Lu(8pv|2*d5Yd-G$@6eUbSkkIAxjc?u@Fm{@IHet^
z4}j$NECj0H^IgHsPLcqy-r3oyFGZ01p#H0ZH(Sm#8=_=fU3VeyZ%|O!Sic&{q$xF;
zdm60P`!mb`%nQaakawjEG_kT*uHBNi=C7bg;j$Wr^c*L7WXFaXCK(5F#?W`O`L#`p
zqQt6-_nn%acKFeF2=hK0*a*O4)n(O{938+LAeht5ehrHciLfgFgZ#ocBWOqXIJ66Z
zc!SDe<>ZW?45Pgmb4t;>$CP+$S%cQ=tu^R%y4ZKp+!QgfZU0L8OXwe=xOVmXczb_t
zZ?AA@D+B=oFa)E7Q0NH5=7El(!2n1ta2~*g-jR{`H944l`kR7P82G;?FWVK2paml3
zrr<EsfC^vqqPndj?EkGRfG-{r(VuE!^2w(nKDrh?K<{4phy0B{8<)hpkiqD{Nv%nB
z{hjp4zG{cn1TaqY)&5?z<;|*k*z@0iC`vI5&Mr3lf;neM_)CG(PMA6@rkRtMlZm>)
z!?D8QTghDyZqaHB1*)}mH271YSCS-Cqi%k=vr8AZAR&1J<z@TyS>0lvt(~p&&UV+@
z9k77P$`oX$5K&`e+7$MMJdW1(>sYv3y(U<rsi??Eb>`b6MSOf?0(^quH<Q33^QzRA
zsS&)v@uQe|&_3eAHwAzNrVLKqdlp1=A=IGRSDTdt`h6LKQ_wNjVm$)`0?7t4HQ%|n
zHqOb#l^`FP8WTf95<XdDH3D|yHQ|F<0N6x|E6Thavb;g)&-3!kjg50A&chNPEr9G~
zD=i6H#mBUZiY0%r+CCKWeAIbRbp0Clg9lJpY2zu+PId;&`tO?nL<3(1=>`%eIwFGM
zV#pL=_&YG5^|dvy{reP8^W4p1?6HC`%4|x>+RDf$@bhw%VSGz9r6hZ{Q2q%ho&p|w
zgDd~3<qw`ch0!ADH*|k7vz-Ipsiro<sE~?#^L_wv9UY7vYOzMWdSyMB{Sh7`M?cd=
zqhD1)0PE0sDG8&NcZlI%Jn}JZt(Wgl4H!+8(sk0gR^8Yz@!-bF%7efwA#X7-Fj6)K
zy3DG;(Q#Y-8a8r|01Lu@MSj4SK5(^KhmvZPW54}*>a}8n@6NMte};n93L4uN8pbev
z%(#QARr55qprD|k;TV!0O1+x<4P)q+Z2DT=BvTg|-F-h8N6Qqo8gun#y2b>}z%byt
zFW31Xs9JE020A)|?er>)#DL|nx13yPP>?@G&qSNwy?fa_S!8Ucd|L3cRmE=T3o=8s
zM4;ZRtw{Rxct>&c=T~Pa)&D_#K+V$p4Y-Uh=ImSr*7kH)u3S7!#EP}FeF6ss4T14G
z7;?%gn3b^Rcht2JU%!@T6AFDNkjP6-;9LObfo&Tf_xe_%4HKW~9{JEA1l}i`v287L
zW#x0(w8)Y7-?v%!`}-PvPUZ_MO@EY4gonH24u*0+d~ruYxU(>6;!8emRT2hFdMfg&
z8W@xvhHgV601yE>sbHpn!<%g>3KLP_e++a|`|Bwyi!v{12te2=W-6wTk(?RY-$PE4
z(*!tO&sI1s>oabqUW=Ky_#=n=j=segxBl~)!eKB{K7~VS=h^}%_);J`;O1ccGTyrv
zl^zU|OPMY<0yH7?<_9u`AO^FiyJ)DinoLsO2DX({T=vh+@59!eU7oGrSRL;qSN}sr
zn(d^4!rl2*zne^lvujwJBmpCzr3*1>Y0&W-ut%fy>nvOlg$eNB5P#R($2y_zb}f3S
z7memQ%=GkaKwY8Aa&s?ybq?s6DqN9D`Z?37BGrFGd+HVz*&UbzWY;V`_*0NIbL&9_
zyFJ#1n5s-@nA7d3fBq1mXtU_%O|S#GDk>{e1f2xyG;ns;7Nu@l%wZKZi*%gz?$e^}
z&IM}CK<(e1TORtzFWM!sZ|I(*^Vxanb>735ShRSnTQLtlwC9Y=7oG;-ic9xH$3aLw
zn~{lw{O031$B+oBo`+;ILAVhUra`QYU9_gJ|2iR}tvf2FNH%XY^vMqi*%0;P=e$7?
z*!)~i-^a8Cb%fEfk=f3@qZ?{tyiTSDDZ<Okn<D?@-wWKL<lY5m{=*ae`+K~&+trV|
z#vH|i+FwEL!*wAf|G8ut=%|~FC6+pr99v}xW!9~tR0|oLOifHoV5mm<jx2SKvaIm!
zck0Y?I0AvhA^5{KwgxePf0yR{4V`IT!+laxVxv!n25E(~QT<dBij<?-#!4>B{P<xM
z!ynt*bKksq@r)})IAa<gep$^{zQXbzgMC+}k%@`8u$e8)isBROf0gEbs;<k$Ip1zP
zlt7EYVAk)11eC?~E-(a)6*CV%pH2;t4qBx^<%m+q<;S$=6r5bLBEC=kg?2!-L5z*e
za}LW>(IuP16$mwqhr<v(+3xnXa~$8A#*G|WL#`6>WsRX2oHgb+`d!!0U%o`7uf-sA
z#p+od-&NT#lqgBkxEYkYlv7hvdlsCZi_#?bb~Aemb?|q7G!wfL(0&OI$C-F=kXB6M
z-28i$QaN{<5+ym88QE)_J4QctU`RQ)plaLX6>Iff_!CoHgo12dG`kaBJHLRZCfGF_
zr)Y90^CuI2Xh6N4d|WVYv3L$mF){Y#Fk0^y7@kpqxLA?bCO&0$#b_B<hp@V2s@hCV
zQii<Ru^Il&NcdpaHH2o}HDyjvQkL53MIa+igM5XyMERBF>aoxX&LsE?+eRaaP$^x$
zfj2l<D*~ycYv_2SBvysR@6p~qqEwc{Q|u94NF<(&bQkSv9hA?x_PqYn8!EoTZ8|56
z&q;jQ_guSLY&R;{u*d{?o{W{`TpQ0^Yty;^p`Dg%DN1K8(ky6~e2F=1Z@Vhw{%||3
zy1&OhWx*<W5H8sgWoX{SHeJL+C#-cD(c?_4QQZlqm-|Af(6Aj!{LA87&M!hfwX0kU
jy-b#voXSyhjxp}wS8cNQ`5F9`Fd`?VC|M?M^y0q(aTiF3
literal 0
HcmV?d00001
diff --git a/km3io/offline.py b/km3io/offline.py
index 8902bc9..335342a 100644
--- a/km3io/offline.py
+++ b/km3io/offline.py
@@ -1,4 +1,5 @@
import uproot
+import numpy as np
import warnings
# 110 MB based on the size of the largest basket found so far in km3net
@@ -179,8 +180,7 @@ class OfflineKeys:
'JSTART_NPE_MIP', 'JSTART_NPE_MIP_TOTAL',
'JSTART_LENGTH_METRES', 'JVETO_NPE', 'JVETO_NUMBER_OF_HITS',
'JENERGY_MUON_RANGE_METRES', 'JENERGY_NOISE_LIKELIHOOD',
- 'JENERGY_NDF', 'JENERGY_NUMBER_OF_HITS', 'JCOPY_Z_M'
- ]
+ 'JENERGY_NDF', 'JENERGY_NUMBER_OF_HITS', 'JCOPY_Z_M']
return self._fit_keys
@property
@@ -320,6 +320,7 @@ class OfflineReader:
self._mc_hits = None
self._mc_tracks = None
self._keys = None
+ self._best_reco = None
self._header = None
def __getitem__(self, item):
@@ -431,6 +432,157 @@ class OfflineReader:
[self._data[key] for key in self.keys.mc_tracks_keys])
return self._mc_tracks
+ @property
+ def best_reco(self):
+ """returns the best reconstructed track fit data. The best fit is defined
+ as the track fit with the maximum reconstruction stages. When "nan" is
+ returned, it means that the reconstruction parameter of interest is not
+ found. for example, in the case of muon simulations: if [1, 2] are the
+ reconstruction stages, then only the fit parameters corresponding to the
+ stages [1, 2] are found in the Offline files, the remaining fit parameters
+ corresponding to the stages 3, 4, 5 are all filled with nan.
+
+ Returns
+ -------
+ numpy recarray
+ a recarray of the best track fit data (reconstruction data).
+ """
+ if self._best_reco is None:
+ keys = ", ".join(self.keys.fit_keys[:-1])
+ empty_fit_info = np.array([match for match in
+ self._find_empty(self.tracks.fitinf)])
+ fit_info = [i for i,j in zip(self.tracks.fitinf,
+ empty_fit_info[:,1]) if j is not None]
+ stages = self._get_max_reco_stages(self.tracks.rec_stages)
+ fit_data = np.array([i[j] for i,j in zip(fit_info, stages[:,2])])
+ rows_size = len(max(fit_data, key=len))
+ equal_size_data = np.vstack([np.hstack([i, np.zeros(rows_size-len(i))
+ + np.nan]) for i in fit_data])
+ self._best_reco = np.core.records.fromarrays(equal_size_data.transpose(),
+ names=keys)
+ return self._best_reco
+
+ def _get_max_reco_stages(self, reco_stages):
+ """find the longest reconstructed track based on the maximum size of
+ reconstructed stages.
+
+ Parameters
+ ----------
+ reco_stages : chunked array
+ chunked array of all the reconstruction stages of all tracks.
+ In km3io, it is accessed with
+ km3io.OfflineReader(my_file).tracks.rec_stages .
+
+ Returns
+ -------
+ numpy array
+ array with 3 columns: *list of the maximum reco_stages
+ *lentgh of the maximum reco_stages
+ *position of the maximum reco_stages
+ """
+ empty_reco_stages = np.array([match for match in
+ self._find_empty(reco_stages)])
+ max_reco_stages = np.array([[max(i, key=len), len(max(i, key=len)),
+ i.index(max(i, key=len))] for i,j in
+ zip(reco_stages, empty_reco_stages[:,1])
+ if j is not None])
+ return max_reco_stages
+
+ def get_reco_fit(self, stages):
+ """construct a numpy recarray of the fit information (reconstruction
+ data) of the tracks reconstructed following the reconstruction stages
+ of interest.
+
+ Parameters
+ ----------
+ stages : list
+ list of reconstruction stages of interest. for example
+ [1, 2, 3, 4, 5].
+
+ Returns
+ -------
+ numpy recarray
+ a recarray of the fit information (reconstruction data) of
+ the tracks of interest.
+
+ Raises
+ ------
+ ValueError
+ ValueError raised when the reconstruction stages of interest
+ are not found in the file.
+ """
+ keys = ", ".join(self.keys.fit_keys[:-1])
+ fit_info = self.tracks.fitinf
+ rec_stages = np.array([match for match in
+ self._find_rec_stages(stages)])
+ if np.all(rec_stages[:,1]==None):
+ raise ValueError("The stages {} are not found in your file."
+ .format(str(stages)))
+ else:
+ fit_data = np.array([i[k] for i,j,k in zip(fit_info,
+ rec_stages[:,0], rec_stages[:,1])
+ if k is not None])
+ rec_array = np.core.records.fromarrays(fit_data.transpose(),
+ names=keys)
+ return rec_array
+
+ def _find_rec_stages(self, stages):
+ """find the index of reconstruction stages of interest in a
+ list of multiple reconstruction stages.
+
+ Parameters
+ ----------
+ stages : list
+ list of reconstruction stages of interest. for example
+ [1, 2, 3, 4, 5].
+
+ Yieldsma
+ ------
+ generator
+ the track id and the index of the reconstruction stages of
+ interest if found. If the reconstruction stages of interest
+ are not found, None is returned as the stages index.
+ """
+ for trk_index, rec_stages in enumerate(self.tracks.rec_stages):
+ try:
+ stages_index = rec_stages.index(stages)
+ except ValueError:
+ stages_index = None
+ yield trk_index, stages_index
+ continue
+
+ yield trk_index, stages_index
+
+ def _find_empty(self, array):
+ """finds empty lists/arrays in an awkward array
+
+ Parameters
+ ----------
+ array : awkward array
+ Awkward array of data of interest. For example:
+ km3io.OfflineReader(my_file).tracks.fitinf .
+
+ Yields
+ ------
+ generator
+ the empty list id and the index of the empty list. When
+ data structure (list) is simply empty, None is written in the
+ corresponding index. However, when data structure (list) is not
+ empty and does not contain an empty list, then False is written in the
+ corresponding index.
+ """
+ for i, rs in enumerate(array):
+ try:
+ if len(rs)==0:
+ j = None
+ if len(rs)!=0:
+ j = rs.index([])
+ except ValueError:
+ j = False # rs not empty but [] not found
+ yield i, j
+ continue
+ yield i, j
+
class OfflineEvents:
"""wrapper for offline events"""
diff --git a/notebooks/OfflineReader_tutorial.ipynb b/notebooks/OfflineReader_tutorial.ipynb
index c1df4f2..2469cd0 100644
--- a/notebooks/OfflineReader_tutorial.ipynb
+++ b/notebooks/OfflineReader_tutorial.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
@@ -32,7 +32,7 @@
{
"data": {
"text/plain": [
- "<km3io.aanet.OfflineReader at 0x7f87bd4ce310>"
+ "<km3io.offline.OfflineReader at 0x7fb60c2dc590>"
]
},
"execution_count": 2,
@@ -245,7 +245,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [1 2 3 ... 8 9 10] at 0x7f87906e86d0>"
+ "<ChunkedArray [1 2 3 ... 8 9 10] at 0x7fb5da350b50>"
]
},
"execution_count": 6,
@@ -265,7 +265,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [44 44 44 ... 44 44 44] at 0x7f87906e87d0>"
+ "<ChunkedArray [44 44 44 ... 44 44 44] at 0x7fb5da3506d0>"
]
},
"execution_count": 7,
@@ -285,7 +285,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [182 183 202 ... 185 185 204] at 0x7f87906e8c10>"
+ "<ChunkedArray [182 183 202 ... 185 185 204] at 0x7fb5da350950>"
]
},
"execution_count": 8,
@@ -305,7 +305,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [176 125 318 ... 84 255 105] at 0x7f87906f9210>"
+ "<ChunkedArray [176 125 318 ... 84 255 105] at 0x7fb5da350350>"
]
},
"execution_count": 9,
@@ -332,7 +332,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [5.170483995038151 4.8283137373023015 5.762051382780177 ... 4.430816798843313 5.541263545158426 4.653960350157523] at 0x7f879052d4d0>"
+ "<ChunkedArray [5.170483995038151 4.8283137373023015 5.762051382780177 ... 4.430816798843313 5.541263545158426 4.653960350157523] at 0x7fb60fcdd2d0>"
]
},
"execution_count": 10,
@@ -449,7 +449,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [[806451572 806451572 806451572 ... 809544061 809544061 809544061] [806451572 806451572 806451572 ... 809524432 809526097 809544061] [806451572 806451572 806451572 ... 809544061 809544061 809544061] ... [806451572 806455814 806465101 ... 809526097 809544058 809544061] [806455814 806455814 806455814 ... 809544061 809544061 809544061] [806455814 806455814 806455814 ... 809544058 809544058 809544061]] at 0x7f87906f9450>"
+ "<ChunkedArray [[806451572 806451572 806451572 ... 809544061 809544061 809544061] [806451572 806451572 806451572 ... 809524432 809526097 809544061] [806451572 806451572 806451572 ... 809544061 809544061 809544061] ... [806451572 806455814 806465101 ... 809526097 809544058 809544061] [806455814 806455814 806455814 ... 809544061 809544061 809544061] [806455814 806455814 806455814 ... 809544058 809544058 809544061]] at 0x7fb5da350fd0>"
]
},
"execution_count": 14,
@@ -469,7 +469,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [[24 30 22 ... 38 26 23] [29 26 22 ... 26 28 24] [27 19 13 ... 27 24 16] ... [22 22 9 ... 27 32 27] [30 32 17 ... 30 24 29] [27 41 36 ... 29 24 28]] at 0x7f87906f9810>"
+ "<ChunkedArray [[24 30 22 ... 38 26 23] [29 26 22 ... 26 28 24] [27 19 13 ... 27 24 16] ... [22 22 9 ... 27 32 27] [30 32 17 ... 30 24 29] [27 41 36 ... 29 24 28]] at 0x7fb5da3a3310>"
]
},
"execution_count": 15,
@@ -645,7 +645,7 @@
{
"data": {
"text/plain": [
- "<ChunkedArray [[294.6407542676734 294.6407542676734 294.6407542676734 ... 67.81221253265059 67.7756405143316 67.77250505700384] [96.75133289411137 96.75133289411137 96.75133289411137 ... 39.21916536442286 39.184645826013806 38.870325146341884] [560.2775306614813 560.2775306614813 560.2775306614813 ... 118.88577278801066 118.72271313687405 117.80785995187605] ... [71.03251451148226 71.03251451148226 71.03251451148226 ... 16.714140573909347 16.444395245214945 16.34639241716669] [326.440133294878 326.440133294878 326.440133294878 ... 87.79818671079849 87.75488082571873 87.74839444768625] [159.77779654216795 159.77779654216795 159.77779654216795 ... 33.8669134999348 33.821631538334984 33.77240735670646]] at 0x7f8790683d50>"
+ "<ChunkedArray [[294.6407542676734 294.6407542676734 294.6407542676734 ... 67.81221253265059 67.7756405143316 67.77250505700384] [96.75133289411137 96.75133289411137 96.75133289411137 ... 39.21916536442286 39.184645826013806 38.870325146341884] [560.2775306614813 560.2775306614813 560.2775306614813 ... 118.88577278801066 118.72271313687405 117.80785995187605] ... [71.03251451148226 71.03251451148226 71.03251451148226 ... 16.714140573909347 16.444395245214945 16.34639241716669] [326.440133294878 326.440133294878 326.440133294878 ... 87.79818671079849 87.75488082571873 87.74839444768625] [159.77779654216795 159.77779654216795 159.77779654216795 ... 33.8669134999348 33.821631538334984 33.77240735670646]] at 0x7fb5da302b50>"
]
},
"execution_count": 21,
@@ -788,6 +788,63 @@
"r[0].tracks[0].lik"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "('JGANDALF_BETA0_RAD',\n",
+ " 'JGANDALF_BETA1_RAD',\n",
+ " 'JGANDALF_CHI2',\n",
+ " 'JGANDALF_NUMBER_OF_HITS',\n",
+ " 'JENERGY_ENERGY',\n",
+ " 'JENERGY_CHI2',\n",
+ " 'JGANDALF_LAMBDA',\n",
+ " 'JGANDALF_NUMBER_OF_ITERATIONS',\n",
+ " 'JSTART_NPE_MIP',\n",
+ " 'JSTART_NPE_MIP_TOTAL',\n",
+ " 'JSTART_LENGTH_METRES',\n",
+ " 'JVETO_NPE',\n",
+ " 'JVETO_NUMBER_OF_HITS',\n",
+ " 'JENERGY_MUON_RANGE_METRES',\n",
+ " 'JENERGY_NOISE_LIKELIHOOD',\n",
+ " 'JENERGY_NDF',\n",
+ " 'JENERGY_NUMBER_OF_HITS')"
+ ]
+ },
+ "execution_count": 26,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "r.tracks.reco.dtype.names"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([99.10458562, 99.10458562, 99.10458562, 37.85515249, 99.10458562,\n",
+ " 7.16916787, 99.10458562, 99.10458562, 49.13672986, 20.35137468])"
+ ]
+ },
+ "execution_count": 27,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "r.tracks.reco[\"JENERGY_ENERGY\"]"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -797,7 +854,7 @@
},
{
"cell_type": "code",
- "execution_count": 26,
+ "execution_count": 28,
"metadata": {},
"outputs": [
{
@@ -806,7 +863,7 @@
"<OfflineHits: 10 parsed elements>"
]
},
- "execution_count": 26,
+ "execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
@@ -824,7 +881,7 @@
},
{
"cell_type": "code",
- "execution_count": 27,
+ "execution_count": 30,
"metadata": {},
"outputs": [
{
@@ -833,7 +890,7 @@
"<OfflineTracks: 10 parsed elements>"
]
},
- "execution_count": 27,
+ "execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
@@ -842,13 +899,6 @@
"r.mc_tracks"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- },
{
"cell_type": "code",
"execution_count": null,
diff --git a/tests/test_offline.py b/tests/test_offline.py
index 185008b..fa708e2 100644
--- a/tests/test_offline.py
+++ b/tests/test_offline.py
@@ -1,4 +1,5 @@
import unittest
+import numpy as np
from pathlib import Path
from km3io.offline import Reader, OfflineEvents, OfflineHits, OfflineTracks
@@ -115,6 +116,7 @@ class TestReader(unittest.TestCase):
class TestOfflineReader(unittest.TestCase):
def setUp(self):
self.r = OfflineReader(OFFLINE_FILE)
+ self.nu = OfflineReader(OFFLINE_NUMUCC)
self.Nevents = 10
self.selected_data = OfflineReader(OFFLINE_FILE,
data=self.r._data[0])._data
@@ -132,6 +134,56 @@ class TestOfflineReader(unittest.TestCase):
# check that there are 10 events
self.assertEqual(Nevents, self.Nevents)
+ def test_find_empty(self):
+ fitinf = self.nu.tracks.fitinf
+ rec_stages = self.nu.tracks.rec_stages
+
+ empty_fitinf = np.array([match for match in
+ self.nu._find_empty(fitinf)])
+ empty_stages = np.array([match for match in
+ self.nu._find_empty(rec_stages)])
+
+ self.assertListEqual(empty_fitinf[:5, 1].tolist(), [23, 14,
+ 14, 4, None])
+ self.assertListEqual(empty_stages[:5, 1].tolist(), [False, False,
+ False, False,
+ None])
+
+ def test_find_rec_stages(self):
+ stages = np.array([match for match in
+ self.nu._find_rec_stages([1, 2, 3, 4, 5])])
+
+ self.assertListEqual(stages[:5, 1].tolist(), [0, 0, 0, 0, None])
+
+ def test_get_reco_fit(self):
+ JGANDALF_BETA0_RAD = [0.0020367251782607574,
+ 0.003306725805622178,
+ 0.0057877124222254885,
+ 0.015581698352185896]
+ reco_fit = self.nu.get_reco_fit([1, 2, 3, 4, 5])['JGANDALF_BETA0_RAD']
+
+ self.assertListEqual(JGANDALF_BETA0_RAD, reco_fit[:4].tolist())
+ with self.assertRaises(ValueError):
+ self.nu.get_reco_fit([1000, 4512, 5625])
+
+ def test_get_max_reco_stages(self):
+ rec_stages = self.nu.tracks.rec_stages
+ max_reco = self.nu._get_max_reco_stages(rec_stages)
+
+ self.assertEqual(len(max_reco.tolist()), 9)
+ self.assertListEqual(max_reco[0].tolist(), [[1, 2, 3, 4, 5],
+ 5, 0])
+
+ def test_best_reco(self):
+ JGANDALF_BETA1_RAD = [0.0014177681261476852,
+ 0.002094094517471032,
+ 0.003923368624980349,
+ 0.009491461076780453]
+ best = self.nu.best_reco
+
+ self.assertEqual(best.size, 9)
+ self.assertEqual(best['JGANDALF_BETA1_RAD'][:4].tolist(), JGANDALF_BETA1_RAD)
+
def test_reading_header(self):
# head is the supported format
head = OfflineReader(OFFLINE_NUMUCC).header
--
GitLab