From 1425683fc864ab22ae96b6074a19174e0f7904a8 Mon Sep 17 00:00:00 2001 From: YiHui Liu Date: Tue, 5 Jul 2022 01:35:54 +0800 Subject: [PATCH] add: GaussFit, read_data;change: getADC --- .gitignore | 2 + Q3D.exe | Bin 78848 -> 112640 bytes include/GaussFit.h | 71 +++++++++++++ include/GaussNewton.h | 8 +- include/LevenbergMarquardt.h | 39 +++---- include/getADC.h | 192 ++--------------------------------- include/utils.h | 20 ---- main.cpp | 138 ++++++++++++++++++++++++- 8 files changed, 244 insertions(+), 226 deletions(-) create mode 100644 include/GaussFit.h delete mode 100644 include/utils.h diff --git a/.gitignore b/.gitignore index 2c45ad5..eaac86d 100644 --- a/.gitignore +++ b/.gitignore @@ -8,5 +8,7 @@ build/ # config .vscode +.vs +*.json *.code-workspace diff --git a/Q3D.exe b/Q3D.exe index 84a4999ff237f503c6fd8d1b6ddae41f94bad0c6..29f843bc6e0a27843cc8d981083af29ed2a794e2 100644 GIT binary patch delta 53995 zcmb@v4SW>U)%d@=*#?2$WD!u^`g=QYr~x4Vlmo zWjhQ;w^r<9D}8w?Pg^T3)nctSFPMb)fNw$2#<#i#(*Bdw{wvL*9d$MJmOlA|p%)HT`SKKHel z>iKZXbzEoNez;{K*N%G*w@l`$`wq9Brk;QXTH?o^Zr#te7?E& z4Jz<$DQNwA3fSu#<~t#9(mB3)exGle3$2;!kDup1Ki~_*C;6|wXvmX(Uy!%HZomBZ z`8xeR^Mq6ZTzmb#`SQ5Z?;HByxu=|6$4w#6P2cnTN+m<;U!C7qoblkfXFL?hy5gd_ zyb19}Dh{SnxYvB2Z}zCW=iX9vi_iB=7#JCSS8@HOKYmT%>`R(>4Dk31t^@v+Qst-x zZZ_3G8tTg?!P@w$z&Xcvd%1oee}*R-t;%?2{hNXF{X>82_sxrCp2Vef>mLJW4D#vD ziOF?--@fBae>^`}9r(lR@kfH^X$KbgeR{)(wI}NKw#19M@lC-|+S3DJD{|w5GH(vN zzdL?gX3gN;pj5YaB);7e|1k5+UFMJH3<_(rza?cIl(M!bMh}VKFla;| zZ}Iy320iZ2cxc)3_c*b0=Te8NVdhm8B>Dmc?-}0vF!G^sf znS11Rf4t!Ms`JMD%J2J~Kk?wTly3KJNb*l_NLFXlqtkhnsJ+%7f9`lAzVYFZW!QbX z{km>{WZ0d0(Mx*s;h=6GZ`#*}%%Wzaxi4shJGJ|s28a6dQV{<%J^YGxe;D*d?G9t- zKHYxNC~7vtO^Z*M=zzZcXq5_;^ekRJX^)sRLl zl%eJL^!lcdUf(LP>JKINKAYwybAKTi(rrpBYS!tA9`0Optxk*dj8;9|CHV`j3~hv8 zcMPw9E()-(%hJQ`Ri|m^l=`%gPs{P^_07Del5o3Ja9viOPHWB$#piOEm-eAd=u|LHNaAZzMZXOvAV zn_l+SvKe~)2l-Ju;gU$V{igl4S=4RX_hgwxeMWOKsQ+2f1ge2;4@3&NcXgK8sP*>2r+p0Q8FKSh&9)QA}V4h3R-POHPE39-WEuc#~_4>Da^<{@$ zj?(K}E->w~K7+2NIO;OP?nq=rAiWv8+-QD37|l3p#%i;4ZjFqidaO2N*c*ieWjB}I zbaL5^WwUO+>4w&1Zp}AjK)vu}JT&Mu>ke8Nt18Hw(NBa=bQ?5&q(A2&+|t7@Y4?8$ z)lK``w0AS(;-$q`8uoTw#zZSzXNG$hkC^BLFA@5MI~f*Tx-;IdU#*&dmG0;lKoae~ z=P9J;B>h=oA^z3rR((5qV7n1Jk{@O3tsFq>W_x*9tAE9=d1pIi@!p_pS7Ms8}-;% zbh}NMw^W!{7v`-12yogv|<{o*7i zdrV?3&$=`HSSslDB;?|5xI@XsI~Jd$+q-ppOX4Eh(bzL8ncIn)M4G*t82-8Pb9pwL zAy-hoWEdpHrwLafghR$w70gZ%PH}u|V(l-5YQy-AA3eS{ux@+2cJOI|NA6hv$lxOX z89%%kwvPR}pbq32_N%(HQr#vVx+%W*gmbN%M}Vng|CLUiuEdrrv|yNT!%2NnyEH2@ zI%(PmO#3Y}+*Vbn+gdVr@*E+-!P>zhc<#=EP8o*6qzU3;em%B0%diKL)S&r=jFIZW z(eSld+VW>qd#pt;*z$rEJg730xyi4)1tM3rCx%g=9y^Y~(KE=5U4w3v0;OnYc|jwY z^w>)8T3(PZKq*NCU!gqIc@M>g-K0ON$7&GUMJ+O*cJ`TeM>J!T6>Zp5d!BB)8MRQb z6B&tUY6%$OolEJITYS$O2Ex}d3-S{SCI}v1gHeLo%AHhTIyZh+_N>@2Y2{(qw8nIT zuK@I!bA~jU;aydyfrN>6aOZ(>RO6o*3;MaUsLf41XIf*>aLSTqxU=f?RMIjnRGyV( z#P}icz_eM>@XpVr3rV`r=h9_UrvI2SslMEwH!;S~fpM1C-Jnt5jDhfywsIUo5wl?F zNUvEz8sAcNs@tOG1GIgo9#g%@2{!9VG|?+HFQa(O*v|{<=x&<9cFz>EqPzh5zZM;f zqJRj7W{SpWXsW%Ll5)@ghUNEhbJhCuqj5>Qv>^h|iTD?yULkmo@3OuLpZ5CYCQ$**I zx%O|vGA0QcAA0-)}*Zxb)naP0uzW;ZaNc>{^wNhC!{s4 z@zXC*NKw5{sM+^9YW97R>b;V#_jBn!{sJc{swaf%2{)}?zaOGo`+|D?;S3GN6Vr0n zk_-PMve)nT-}5SBs8lqssNEZHXx&;m9dk!_eqAtM%3AipYNWPq-<_p9_b`}|fdSK= zi2N-_jG>Wd>tny2HMly{9II7bh6ssYH~3UEaj*73M~1M`Juu8p8S2_ouO{j`gdUS9 ztGr;oX=g|&jCWO(8te^CLJy$DCO`8{tgi*8_LOGEeqYcfrj*@m#{N{$C~ZO#%uN8LU<=Tkl1@N9PqRp_H z#DGf65z`)zuscBvdhIDj7(GR)exk=7mfoV?RYy$ZTB!&(-Ed@ZitOR1LS_nC+Khwy z_LBGjOf{yU?2pvLsO-|F9Wug)wfmZRM5VpR1#V7W99(2ga)PBY9dNQY8i&wh6OKgf z-Y8?CIT0|L|C(tW`U`?nPd1s(koba=hgjL0;HXgW6K&d#16t)r+O&4McD(!4E2B2Y z$k?u+l^>CySFca_tJ_q+Iz*#{Ut5ZTX)&~?TBTA9T%$-NfM3tp>}GSbd>|dRRo4+= z1^LmU7tPq6brwmd6;zn92vJ z31A6;J^*V0B(zXuHP3mJ7VSdBM0KE9KsU_OhO}_MfxkkYJGhoa8$PVg)MhmGoQH)2 zk(@}kREX9=N9c}fq)?(Hk$s7fy0UFvnRap0qT?xfd|tU$*rX@hjpQ~xd{kR`{Ly3* zkjn*Hq01gxq026&o#G&uqxA5?!m4wzGV+tTO*bI2EVM5+#`Z*9Xc!xxFQW0@d_A_X zux)&~kMt!feQ`4PK}nxV{(Cohu}WTSI+Ng`g(aplzPPO2FX~Kprsj(s@?Ww~aP}2L|M4}cw5d8%Z`f0f zMx6zM>Ngg)CVralkKdR*X6VoURf6kD@MQKGc`wfN`?O~|X-hIU#=W{(Gt(dcOZG#7 zJ+1MDC*Bu$sb&3<6MveS@xWAne9b9&=WfhVBP=tUE`#&PZKe~w%@oVzRq4||$HsS` za<-LuI;NOdKc>Cch)oU2eC?}(yhA71Z$#Q53>(g|VBN z{0`BapABMIHEhyW{g4|B9+6c>?Dj6hzAF@s%|TJk!_X;5=()11(u~b5D7HjQEfEX5 z&&cx=Z8WU}Jx9Tj)s>;*VgdbsQ>5m$RJ~6}H+5-=$@^?nJG4MG_ z2h|8{(K8fep7)wcz(CvJpwU$k(u9MV>zH-#iI)#j+{ZO#)4z1_7gN}C!L`G4@i%`_CKQrtPg_aP{9Ze~k z%b<(e?+MLiY^g!FR}H!^jXAB`|~ zw7MfqMr}SLK4%ExqCF+!gj6eURa9KA2Be-$0%R$S zGbMFLr)F-`Gw=b3UP!i@8HZhQE#la1h9`k(iy3}Q$#iYyeTS3G&RrobMuUiLoof?U z)BbsZC~P#esJ613#8IYebDr9X4TE}I?gVvtu$r{N!0#T`=q}-9B)gJ*0D-9n2Y<`MVXE2MMzQ{Pr@+T5*n4cq_O46_(12Ni$T1 zrrc+o!=l95<1)5h0wZbm+OtVs4%&#{yzwTI+ z!CadgYl~^`)C#*z?b=<;b>UYr%$vS6?N+TYVUW&rrUYCA{>#&Cqe8E|jDJ_&77DX= zG0(!mze8I&O{S5Iw>%@fWocRewCsQB!4|Q7v};@RB!j$@yfxePSO)AKJB6S20Hr43 zdM2#iKWxA0cgJ#7A$~yP`I(uknf|swPQg*gjHih&q$H zqK>2%8`#K7RCcPG=eMpu>%7|nr{sy{n*Y@qQCp`kpP2TPl84ScH9qOQ<1J*`7^TCb z8BN;qw&~@D%Koh8Dp(R`(GlI56lA3!!(t-2C0cX@su)GD zMl%i@EQm2o&XIKttO1#y&e$rp1SJPGu?Dn|EFCrHdvHLm{nn{OzSs>Oai zV3jYOrp-XIhGzBwF3-WR-pkFboW6MT)4A4fY1-0@b$dcFUTn;nrf3nf{UJ0v6u|RI zYVyvMwK6FSTI};zzaxI~Jo+bgVxefJ;+t>cp{%%UgUeSlL#8wIMcVFm_-dwI7+yoq zl)B5bnP)^<9$D4KYs(O?*(!%ymQ5Hkkmw~k{-Sgo69etM#VqP#2@@yQ5tH7#J;#uc-`>q%oVF?XMEIoW5%I9>C^T25l)9$uuomyuC1Y! zE?2uOoeb7aLbT}4j3d#a*Y)Q8(7MHN#ud+7|HOH3_^rI@$N&sd9Ac%&6Y*xV2&lZ+ z>Z$(i3_SxYdkU*#eE}T}>5TNo-zIJOOSFQqS0Y3GYxxZ2Vy*BcF^{!Eg#XOr@l8~o zz|CMt!?=4=uw3{yfLlWHYEH=F{_6~n&w{lJ#dim-ryo)3g*6&P)e%El%xszIU(rWK zvgH^zmR5)hcV=lXgQ)yCerJ{slKT!9j=QTecpT8+Af^WSJarrR7ppfZ*)Izop?Q{8 zxz!ABtIgA!-^nEDkkh4IL(Y`jGsv-7DX&~y!<8MvXwGiB0Ah@yEi7mYf2m>jBcOj6{fFT~y zlv{w9xaXOasXUcB_cO&24?HCBQk7DGAyC*P>_^hgs*^vUwKao5gz}5r*oyg#+-i1f z4@z#SF~~%rzDb^Fu*@Xru9rcG`J4<3q3zf6qxRQIn6i=QbO4oSa&v|mZq-(8B#B@{ z&6h>|nD*wVQvtVa@%j%|n4AxNSPEl39px?~%FOE13SXxGV1b`$(yvP0Lq%nUEH_9R zxsZTgHwB_l2AqqSF3F-^^6X(i4L*f}hmCO8(h*7k2KSgOA**`(`;JU%vZU`wBI zI0(8it-zHY7ajQ@4Gz=pN=-2FV`|IA^wvUf(ab%<{c3 zspt)1W2VNHU3d>$gd{-^&|!pIYTo3{?r3;(ZH@t#if%Y1oj#=)L@#tK{f=TyaUBI~ zrZ7iR1{%l>ei)@yFglXtRe1BfYdR z){<_{F|CnQsxg9(RAX4St!XmCO-s*0oa4rz0r+#3M)`bWQfAoJt=JJE3jG}7b_Xm-yP|PHVp-h;7O6sea1*^ zsp%Fm5)HT15=oOHUn;p+#wp=~R{uV{L4VG=fkyOe^_sbcy-9b*g@|E*0%&^HR+{%A z$eU5xE9NzlrJ|Rg(mrh(z8r+dG@(Qmo#3fc-f-Va$bga&I|YcGpd5eoMl@7S5INCG z(Jd%!QY+%BGDIwps0Xp2XzS`^k_>sv{#_ zkVOWa2FCv|?DBJmexV+ufoBX*Pit%Zq)`_Po54D}KKeyH7D4j#iUP1_Dg*cuSmd%t zi*W&8H)`71WmzOEdP~-8O5qYsukY>!Sv{vFbMGqj$D2mYE&J9_kvjUaq(pKplEvg? z?gtmBb@FJwK6+asd7SutZ>-SmYRu@RvrW5tzE)UWtUI@sGEZ#$Fn;u+)8g-4aB-mP z*7&J~cb~a&kXNV7x9gj!^R3>Vlask0pHK8v;a%tVtQlA#)6Zn?{PR_X&I<4$buJm? zk6$rjRNzog{I(HS1}4po|7^rvXNNPrf@RgEzMX>K@9oKT3mEQ?Uwq+*=Xd{@vW;a& zQ=K8A>&Aj(%f2({kB=)V3jBRv{98q5-?UDubj0fu`w6h3(3y-L2cFV^J`D&KH5jf<>{YfJ+R!B9wLi$dj& zg-|{mteI;%-~23M|HZUlg@5AK$cfoTcw5cyyxj4fBZo}hF7vSP<{AQNA|<-5XK~Gk zVdz_wdzSRXjTK>TBWo`?Rd^_w5ETkf$`m0VKLNf1-rG|{iLY*mUwqN&)-1!het@8o zbs_C37V)lSF*LZza3&KI^;8GXETZzvVwQS6$b0Me-ZK3OZi0kC-C2}x$TfVp=69@% zEoH6#H4=3icJL^PcFJN6@AetXK7*zaI8@R5_vgACS}EMp?gziNOFi4Ap6ybv5b8O% zE>LTw1S!F;litTZtS-7Du;%9YZ;HNt=41HZ86c(uPtg!)15ZJlt@-}-MI!_Lz@ML{ z<3H5S!=uos+tEU)@bYB5>!N}{)4Q<#6@k-kj9)S8_Ol1WbpK+gU2(P=YE}u>>F0*p z(?r-vxoW6&(8qr>)Mmcpg$Ie?mW4uC%=vF27gJZnp%|uU6DGG&n{u*hz8!U@pyBFB z0&@i|k#uelfnk?-inw4sA6WTLneOhSjPlNqK;@-?W!}3pB;Y(390F|aWcpdYQ$fs5 zqPALyU)WXixWsLoNj3eoxD|*!9o|{f@>2I4tKxa+Y0JAB+{gQI) z-?gwDljMW)gLC<@-L%q@T|8`MMKbo#+@T zD44#)j7@Jc?H|>V&y3xOgi+0rWoO-16^odp{oRHv7))P-SYm>Waew1#Dix(T2Ke)9 z$j@rKzjo4fPW-Xa1t%&~B^8JnZyJ5(+_lNv>rP`(ey&xQ(kh%+|IjSuM4DAF+jJT_ zl+-4)3lDgrWErEITp5m^2Q@{A5M>)FCn&Q;uN);q#&X8gzNBf@ZiJmLdZb2XxR9k%Z6m7 zDc^h9nS-GG*||iOUw(?^1LrAfm+B6yhj3$KqIm=bqLjgFWP{wy=LWun$Y2$Z=@OH$ z=gGu>5Fs{T4Ou^wL^Auel`k`fjc{KT9D;R}%w01C_0y27)x@RW`Ydi&2v*u-WlDTW zsf=u0EuBvj?IixIAJOtClLg`1RO3+}O?t)_SU$R1&xq2n2+Wk-p~td`O6MJRz?Z%@ zUq?FfUgB$xNnoW4Fkgp1R0_Nt$JPLSgNTlt^4C!Tt*C(A+RBAWBz7%5|I3<`$`E{A zPqs4ktdYj8N#_3e6d2yHBXCxOa5zO;uyDu`J{NYLEyl?LvMWZX%1E7zIMCzB0E7gg zsJ+$1#)V1`E07ZnC0_d|nFJ3C&ta)!-$hcO8pDx%H42KO%7HvFWc*Q_)BHmyZK^sm zo?UXXKQDex$*|&of^fp2#3q3S@)Nh!z=8utKA@KYsSePR%J`inr&+}!wN<>gLtLHD zm-zGfDhbgEI@%|qnAni8dWuv_iQFI=lz1ROi#}GtTW)3nFKA=UAJ)UdRI=u&msQ_X z;mO={Plk(ZN2jJ_?vRuHC%7g>GWYnC{qal2j zST%zpDmFw&%m$qO15Dm&%H-XERPJDnE%B>g_s!_H#t`^L1|{7zLli+N0||eMsEeFv0=?fyjpzp%V7N`A zbu#u!k>ZVlDH5SiND*15Ft-3^*TEJ3PQiu~)f2b(B$NBTE1eXrBS^jkwk`DQ7V`?> zEr`6`v{A}L^hWCvE6-C}kKHOF>)`7p!s@ypRTM`Rx{_o9o2}SiVt+*XGa+H%(6+?) zcB)|zts{E*BpCzI70KLxa@S}&kEtpR`-XfYWG^gcAI8F7!+wVq*c)5de`D-#{eio? z;zPp|0-t=i{?>4*KQLiu{E4f+6*xuw%t z8=|u;-@qiraAXNuw-lz)bvyIN5vox0%skybY}(8z7hxJQ`u-?rwZE2?)+$UeyR`d7 zshRd*30lgOdHvUM`n(ZeG-{-^@~>oYXa4!0oB2O06Lhvd%+>9iK_&s>i@$lf8X=k@&P@^nQ^*lFU(SMM|O1R_IHlyd7pN_xK%&j^UtLDvWbCBcXgwN zU3|xT{}Km}iZ1~i%+YQ3wb&;~b@82}PMyG}gHvQld~%cO@j$nTNuG3%Zp zA8*FG*jsvQpy@nzSaaobwO5K0l<5TssolAg@1z!CA__;cdh%jN{aFSsdx+n@4w z@fx@(XJ@!6Ul=zhFzNaDdE=*@QH^(kwx}}i38Yir74WZraQxD&jNkY%WXw~n)qlsP zShlen#FX7EvBCrXm5S(>wbx%46iNuivJbkH7|+bWpwvgBAsQrHmVjW4+31=myG~?a zo3gMBQrUjF3_EiplzHF?N^b zstN6>%_2h)rl!R2&BcIprXyCP3j~6g6JyPkAl7K5$lyXu=|)OM63>=r9W%O)=fu!I z$jUql#*VDjGyZ}2;MNx#vsw?osV)BzZKiRC{W4Cv6;Y;hC__xG=v?nh&HTi$+w`LLdom;#ac7_CmqJW$!atm= ztO<#RE}85nmktN>rAWPXSigM|9*6UR7*&}o6BkzB&hC|D%^umiBn7YmBQc4bvXVWY z0AHdGmadiw7lGedbb*bVpp)vKddbLO4nUCMY4kgw#|U!FVfRSeQnrsB#|~4BDe z3HBxO%?>iX6^kjjS&s!byN)=BbR01$dd7BgLyd4-4a!!FI+0|QLF?R3{RHEq<0+<9 zLb%ypgkjVfb%yRHzzjb@H_!x{sTY#St1OA=hKr|nfbM;u7|PukB4Hj0^Uc_;g(NH{ zz*t*V#gtWT>WUJeXIsdO=ey~nDPb_72Dvrt%Y|LwZm1xHFIMa!d{CiDx>nl>Us&!SFePq+6;5|7ctgZx>0eTM*$XZ^gVh~gPz)IBz@aFnNN|hcT$ki- zW{?QN77#$J%2To@C?^OycaU_uO7~KK8s&;?Q+svF+r(kymV$I^3#!bi(xjB~wvZG$ zunZTRXRtgI4YzA6=VHSt46!h{l%0fLGZkUo(%rhi!7y6^-jNbeUAv@=lcNxTy`k>HKlY-l8dZRj4Xz+2OqGl}!B&C^4JX1L znnO|hkXrF4fE1@*o^j|+FE-b5vKvJ5ro;q;Q*eR3ICEStp|7FrUT!Dy%wh^}7`roJ zU?Ozoomiz4vrL$_z9nBJNkw=(tB^qe_EP4<%}}6SMMFs0n-QK;ijl!;j8)mmcG#Nt z8F#Bio~zNsq+3kYhXoyY4@D@WR2UvIC>R;D2nuBo$n<2@F4BNJnQ@8n4$l%3guG`z zfnr1Pcb6)0gnjFCt#vW>wa|1dqotMzqiAtIPg1h8282jwPZoTiX4;k5N=>d1Q(h5~ zR$`;AYT*XUUD;NOTURvzfdV(Md@mu=P896_!(hcKYX3tQ>}MoFL^}j9oQZix(rU)i zSBpgh!oF#eghbCI9zzx`4{%#16aZhL!8{Ou;Zh?Ql2=)J2K$DZg#e_GIA+$Zp3`B( zbCKR%I7a@Wh}weuBR!G`917lBYu-oxi_EYi9SZf*;Z9ysZceC&fQl%q7SR&MWg$3| z`4d>oXDZ%HKf#o^r*y}WmkN(0)D?134T8x^^?+vNh05Y-0Ti!fFbk|o-h|Xsj+4(d zvGZveNvcsxJMs6?FcD^xp}_z)xQalB_A+gngx9r|X^*tYMu_7h%dw~1C`GXb#;H$1 z-sFZZcI)sf@)(=GP8LqAlr}+=HMlQ4&2nB!v!IDYrrAWMMwe7p;-9)5G_80C%ZA7wOW9SCT!q)xo5O_026o_13Q z%BBma3%jD;tYOmuf$KQ>#fKJMiQR^s=IgnH%nEjN(fKLZSu816=aKH{ zqRdF_{GJn4mg6kX@-vy=8TK0zWA8jk8DgL8+y757i6UzZ2d}xk6tib0uh!c{H zs>i%uir*xxi9srH!%gb?^MVepI-3| zc3I+LVvBVz0xcx*fq87HNgsfq&ZyWO8Qo-v74J7=Whn7778e)eJCj|4Gg$F3n02-o zMQxLvkyn`EzQxzGe7HEOg)HpdRxKw(3t{NyFbR)M2#ih))0k@uO)ZXydWMbL|T$TMj%1$4zVX0jKZ3asByQ1z48`}mW zL$B}Pgn<_qYtNo4$vG4vjLI0d4{My(Zb|%!Yu6jdiSJ+vKtBAOT1YepMB~X@rmYyz>%NC z0=rC6-)mZ?9wuiGX5rc}bRbSaS2FQ%+pzchNhKk&?i_=@Y+ z1cKMbFPwa9;OYJG2PY31-SDy$>SVU>SACMUkW0+~fi?@Yy(|P8=rP*C$^7yDxL#4{ zcXQ`@xvekFea{zje|3MnYw`twSN2Kef&YFfet*TVz|{BSzp5yhbrFlM(m&1z80Oxd zkj%|fcW#qZoNvC`s%o%X?lZgzl0H9~d%(RVeL~&2>G7c;?v799zVz4lk&1$gp7h=o zrslW%CKaW9oigKfGZ>*{?sxwhziP@b>n3st#*?|z29zsFQ~w6E>3DGpUnA?_|H0Tsl=SV-n>W`NLdyG}j*$$gqMKbiYm z@3w*yS-jg4X%pNqRsT=gbe6Q~iDs`&<4baBQ7M0llcZH8(yF*Wdq})t%B7b^p^Y9> zVMB?Fc1ok_**7?_&56c8`_!p9cs-N3k$2bcoAN&wSa$SUv=Ui^MuJ* z;F$2JTooS|yDXeSwz|3q!Ezgyh2ptBgaeu=91M-kgt2Qg&iVkK52i^Me?V`Npx?wO zYT-FUVb@BL-z!!`VV0c#GFux_RZ#J)v_l&q%X0SXO!Cp;{F(*%$kvX;`qj$Xm`ZSH zA&!Vyld_aUUA&zqk(L(oz zw8D>g(W@1*jYA5M?YjKUopf+diQr zvAkurGV>D*992f4lwr{d5{jHy0G{cV>S7^XR2j6H?OGw-mj|s7j7V0gNtko?&0r5) z>0(UeCmwnUwd&-)1iQHa`w2K?0MAZAHXaV$eUYr(a~Oe~Co@%!D8kDB!htfwg#DuU z9UAra@vh8iAFpJ5TLxm?9T>ZxqXIQhofDsY9wYDj4|_bB?VNZu<0(-!0Pq>iocLM_ zV0p*cjMANnoLsVb(bb9S9|++(#L1U~-NNZF-Q9Xj@1wdx82&4&wuW6**)jZQWv z#S#@JXExcp6|?7SGw|%<_Emg+HvPAEh`~v~q8w6XhMO1Nq1$YrV2`=|s+tVQN=>M0 zQ&tjle^xKdhe6$DC7yp=wA^_3sIm{Id#4R;E|(Gr%k zgyTvTPvVf?Sb_c;->0ywjxR{>11YG(5m|?MF1|XkbWAK3=N7h^<1LaO3op+-1U$HF z@x`)`hhCHeDMH4dv9#nLr;9XI4VUv!*uUFHpe6fvF{0SdkHy!?T9fc}apIJrGB4aH zOBjmnixbP{sTBn=?=k4}WsXrI9us^)_*!*mamcjBPRNqO6HclcEw!8@{>wAGs@N*r zq&?oWyh&>~gZdH_upqqZ&2>`aZ)C|Vev&ARJv3)Tki7ZZocj`G>OdQet|59 zEMbavl0^7O@VO=KDfid}dkfrDn?J8n)}r%e-YX$J9cn>0bzp*?J}_Z+;^7~7OQI;* zd&Ji)_$-M!FVs<+;k+PfD23R3mIY?B_8>dQ|Dnq788LcOGPlekqwHnJU$s3k_p68r znwj5&@y$-3S<(Z}`H-*spRG_OF2^PKr0RIxNS0pao=p_k%v79m zBAr~BsD`55)D{ItAR17iRJj0132Ev8=#r*U(p1VIxt2ht0vG?bEcFvT2eFn&do0b0 zP}_}-Ob=Ik&d0Tk(*!eK)lpb$cH$z8eyy-qA&zlUOv%MP7fSVX42ttCj#I~xEMo0} z@7IJY^+d~k9uI;_4|QPFl$a7Z{RiWgLsMShs4<~}L~JJBV#ihwHfk$cR`joTRFAE( zMBYjH@RPbgD7Y-mTf%1@@;c3DGLF3s4yltlcGX_sP7ZZjWn?xK?1IbKS}POK*Owm8 zj6_QB`HM7vcH*+{drXG;b~>#jq}?6JOjTrDr%o=BU51;gc%NwH7ilZplP5B*RAd@r z%PcA@G3OYCBA%iU#oe4A53d?xO{E}4mdtZ1WJtsE6~cIxi7U|Zy1lSeFIrflHy6$uqCMj=#?>w&5vUDK zmbM%X zo~3N_mM=NBaAKdFB3R#q3g+0tA%alZ+1HVX+%J}atC!!Se04BHTKyuQW0m>t$yn}u z)>)hQ45}^rpqrVmHe$gBZN!qb=(y6=oREdH<{o4b%LGNoi8e~fvr>IdFXtq{M(2tc z`@t_=CQwrB*~zXHEA7dYm77xGBt4aJ37+m1c98XC)>pa|6E< z7pzWrN}~{}qBP{9Hal-mTy>kc`ywl-%XQ9vAofAd(h3k7na_|wIycmTLpo=vX04TO zqE=Z2f><;(+C)Vi0FYhJi${{K%$9-KL~mxNU6QOp>>vQgBjG!2bz26cRm>8#-i8KEr?%Z(H1qj9QFu{-++3Gs=60$AEyOdlc6Y>D^a zT1j{V-jzGEJfw!=m#jM5`X>u3uH_y4 z!~{)N04HdY@%|&}8Pr6XsGhLEQ-WAuAyB2_qPkKd(|w2_7MNNgYTs3Yr7>;%&Ma0- z4yOI^5-Bpgvvw^0!&kr(^|H9-=D~H{hQ?SfU|-z-z$w;FgD~{>;Yz04pZ$WyqY&L# z&mWRRXow+l=NdV>1u?HRjHPt3s<55JUXxxoCgz5OD}?c;uT_ABrl%&Mt}re^Cel5@ zC#qRIeuTrtqU#Q1vaWl|*M@>ZERl&wC+df;cI-S-O~ zqpWZ_(Y1-%?aO#1S}y<3L7jUp#1@Hl?ixmDZM6`5=q!l9M8L==O=ZKNpY2 z#r(=7%*lQ63Fd1}m$Qm$XH!GSfZsckxhJI)tdRtd4@eMF2|oOpFjwyR6Gdu}LCO`| zOe@TMi1soh8=-l-m?oBQF1KdvNwI0%5v0r>S__Xk-9WjjB?fu+!7TC|a<6J(fnZwE z+al)~oonRcWQN}t?Vy65d0Of|tyAw@^NEJ+DPQU?;lk*=~M#(V0c3!EeAJ$0+6GWF!O{o_A*mfY-QZlRkD#h1(-nz}Dd z=7!w+Rsya{3Pc|haaGAoaiL5X@9mCgh2wFizmwX6!GU$Hi6`gNoWv9A_uND>x&AK? zp5s562B@GD=5fqtowNRSPCOL9;rO>DaAReC$6M>G;t{`fN}XV17P7rwe4zJ|>wf;? z{QZi*Kk(Pe-wysFZ~1(WbKS$=`~3Zbzfbwg1bz~Kr|>tFzv29i;_nLnuH$b8f3x}f z27h<+cMpHQ)BUFh_=}f3ywW<>82NRU@cjyOoMe=fv4RkxS-GE*#JS?Gu9}2WN3&<8SBOGmcx% z)p?t2tOVn-=~znBZN(1{HBpCm3UR5NO@*dLkLRAJirbRdb4fDkfu{i;;(<5zgOOml znI8Dbe()^72iHj8KCbc~^uuQZe;>Gjt$y%az+E1Aem{6V;7uNQN0>SN4M! z0REK+KFGRv66Mx;(2g{8IneKT(BGw@ zD}YvcQ0uWYcqQPwJn(&KXd}@19&|w(`hB2Z_n==(LstWxoL)QSU zSc9AsOy{M+KL^wl&J^^{xUL7URh-vhnWgZ?rN z-Jn}e?nn>({WSPbBtGAR)~2CNK=VE5{4^AYOYRvS^!hYZw8Rh(dQ}>_73fJG^n4d; zp~O1@XL;aII`PXu4?pg5;^9%LK6e5Az=OVZ^`I}Lp`}1W9&{~GOG(d2IpDz__@Q)S9ca*l z)})~((4#+g8FOA5S^>1rgI1)W(}3>xpjW1$l|c8pP%BZ82G0V#-AjB@8af+jlLzf9 zO7(dz(BFB`x6{!1K!5H*x2B+Q(%>q<-|-TEFAZG+^bQYN zm4?;;o##PsNaettaP)C^Ijaj z5OY|j+s^KX4<}Zf*aoG!53S7r8*3FP} zEPIf|s_NQ>a#rg#$%hh337=wLi;LZbdRW?)~Np*XA>GYTa5Y- z{nd|)1E!dnxH1Ih9Mr>j$z>~%E>=#7wlXTlg{=%ooi+c*7WsxE{~s&p^Hc1yjbd~R zG_K-*et3*FLL37eLzk=D+Y)9t9{52baJ>=#;0Jkuyock19x2Rv;cH@oZTMRJ@<%QR zTytLh=0|b^U1Q^m9x2LhrlS~((e}iYvGFG#IWPMm5Aa%I*x2~qN6rh(sf+tpSKV+L z(nb=pUG}lKo98%=I`wYF@r#-sPe|Nm&(X#^AJQRQb(VzOi9sk}BhX!gNjIaISM}c=rH4{afQt{Aiv3PF34ePJHo?f9k(O)%LfQ@sh_s3M_mvUbJRu;N4B}U#_VRM1K)K z_b0W(mVR$QZ3~jQcgi!hHF>oy{if75_{e0b?VtR{_pH7$Frg!U#ZTAy=c~HTiN%lm z*-v?pp53%U)%8wO{JEb~SDUKq`Jbq|K5RI)t~ca)KwTf+Ds^SuK1u4jb<+CSFWNLK zb->oL4-zzqnLNDsjgO-CTlmXhd(*k+6Vv(j5uWgSiG%0Y>P#$?K;Qc2&(t|n;`z?P zkbh5n?S+Jh$kERE@Es=_FflSYCMKc{c~{)p#WE=_VKK0d%9iFvMq8}T3}GU&nK`Yt zYWN9wSuZ4*P5BJenR0x{VX4pe$B^L5US5z#vD&!JoTaMAvT?wBtU-lcl7gd6g-0Z! zTj(iKGsTia!n%1@G2LA0jeAOcxus&cT@cfCUf}68@Qo4pLCfb+sm4^@9x&w0mG4}~M^Ftm!<$>DA1>t5^Rx z>gVxgd87g3Dzu!hb<0Dimh<&q9EsNhxlqfwJEY}&hjsXzaV*K_Oe|$HR~^R=X)Eh- zW~8XxGYCx6+Nd0zdLuGdqWEVdp;Z%qo1%8BBEvpt7`4a1bS&{%?EfU>{T4;n@qWYJ zk3$I`J1e8_5#>r?u5&Lzrp=R`tB&f2-kj_doQ$g%0<{`N1k4OP{i|Nc#3FvE{6M$d z-rOg@7s%l;yAR~5J%RY)SmVVvIN7oOrYR>PSpyd)>^;()^s`kHZw$l8&xl#F?o zrJ7m&PP3KajH&OB_XJbdE7oZ4=Jsp2+CsyGKs)8}n?&#uKplVgPX!kJ-G#8xQ}z>Qxth|Uio-ri;zTX zz{)0a^&d;W?xAAHAvgk&r(_@g;I?DAK`+AF7~(9rmaO;(&!1xvWK1<3s%S;y_Wzi) z{iJM4lhP&PATx}<(i}TY#uN!B0qLjU2hbB0CT9W4E{u+-6Kscw;_xhHODV!3B$PxP zu>8zf)`;=V>=@yEHwgbZ@gQ^*X;2~JVSlgORx0t3&k-@|aX!L9&R82L9O98bpJoZM zb_q>B0yCedB)H$NMu9uZ<;=L07zlEJAtx)lfQU*BVFyiFy~; zwMPO*50V<0t%BkYy%Ei5Vk=Px-VY8`JLZ(9F}Coph#IhsyaA8uz-oXAI>4D{p^1I{ zKwk=$rN=KcPS7PY-=SdH!p+8O@1SaeQjOh!RMdin4^Aud|{c z2R(~PYTA5`t#6g%w)Kp_kDJU5+z>B+@rgjrmGKW>92y9JEq+4V$iN$q#`8N*58VA|d|F#^p!m`F zeQj3_KI0a(B6;s@6-Wr&F zSA5jATZjM8^=i6skdL3)>cABiy#>V6d@v<*p~h9ubh=9VCTebtneh>y_gRY zT;0qqXD3e~LW>uV$%`jggynv%+}}*1ePei;cwZw>r>ZhZf%j|_iKAJ9)3U&-9H0-~RFUOzV! z#((izv32=LfD^^?OGmozrz3}kOK)H1mfjj}Z=aGEevF;MF`|;iN@~5jvBfRUSrlVV&PiZZ`u6TpmO_t9`9YvuxI#OeKYe& zll25n8^BJHu$cu8XB03y)8&VyNOdgGw1W9*3B%dFSr*^IF}2D2IJW*`;nXF3OGR6y zB59nZNj9}*aen;3Imxw%%u#LM8Kf=$8x3M-h_>oIxgjQzl|BNVSY?$jgYgX(mhVG6 zt}}}GFr|tNSel2?0B8fv>?27?|E3(lz?o$nKSuHnb`5ux9c(+UdWdYjY)-ISrc2Ig z54X#VOyYzGZkUO#HSN@XnQ!CLk^M4P4xB`>>?;P-^aRuXj`)K7v#mF4IJ=v@F)YNg zHRe1rahi-w6IsCCApO9>!U}I2qV0&D++oyzfJO5GBF1yS5FjfX)L4{PUw8yz652nT zHQ6<)5O+n8r91k}$IK?N19k{_A>L6WvT4PlYXI$Hk$p8D>(%{7O*qjNAYg^T*E!oG zjj)JX4a!gK+#z+lrDrUyV>>CKf~v)O{olA%t){`HrmgBWGu}}|e>O`MA??MBR4GAy z@PPU*Ww=jU*+a05h&nY>CKjud8s)tDN_G9c`UH}QK9LF{gcooX(YHihA5-5>Lfi>2 z%NCO_W%)8;dR=hGI+(b_V$(RpMRuJa7r`L=#bwdkWp!OPj>`y|O-hDGx9R{Rq)UCV zC3ZW*V_}mT08o_kP!#j4c$A1XV*Cmfz&(1dKU{JU}J4Ws$J z3?rjgcFitVdL|ZGpaK&Z$$aJGTAV5>^k$lP0vSW>oG#kL`Gg&G5QMm{tPCRkywkVT znl*cJWoq?CY}qaFwV)ChV6YYsbGO9G?5Eh4_`KG?_-&wK~N-z z{bc-;P3T5=ueNe9=BbP)ox`{kV%$VV{omnlg+#>*Qkxz{;5=|1Jf#$NsZ99c;tF+z zYq^TnU!tzRRvY$(dspx*+}kLp%7~+MwOnzOaw)sn-6Mq@*Hu_qNKU}fV@vYIO_mBS zTjH{@R6uZ92f<|<;`d#Dfu-A-)9NT|w499K+6C;W`cSy>qIyFUu6`WItc28zaE7JNaBX06$;vm zHS3sAGa?rVaFVFAN^{QmLIa*u|5W5NfmUt(kmW)m`=nI_592#b`eM|N3Nf4J?j;fZM;gh_Zvxada5zSkO zIFW(2CDG_4llAvnBN5MW;V9=`SaxYCf4Ym!NT)=RN#uouiE?h{fZiXawt>Bmb|!ZI z4@~Gp;6Wm8dHx8YPlYTha+DNtW?hKvJT@D4OGk+Y;^WZxHLZdL@)5wWRB@(+C$XPc z)iWUId$p>}J1vS~mbgu-!Y*E*it+-5)G?eaZOPeBoNLBmRj!U?3oJMEPei3?$z$S! ztUfm+*9#CjOIjRa;&N=vH%ClG@_XKO1AuxHl?ASkzB8G7^)=W%W%fV$h75Eltc|KB z$#DnBJ@6r~hh*xexE^p)(QtR7^Jy)IOB@xMdrBw01FcUtgP zr-baDx>j&{9Rh<-lSYcMCrYp{iYXslN@zAwao8YZx$z-91&!AcO6Eq-r-SD`G=>M- z5VSwChDji}h^%9-EpyGt3TAPFWWE@RBg*|0px%X!8qR}JRn?e2Gz})Y zHWkBH2WvJ;_f8d_NsUa!GY{N=1=%0tR+qT?2}LbVq;WtR5}?%M?={El;61; zVjwB7>USifZ;PaetdUR2-ARt6YstAz<&=YMyxi+V%&5#)56nE<&Ad)!CezXl?U5#g zp06#k5qfj1VH83C24ll?);Ef`g#ERR`dEZf6%!baw5`ta>BxL0^C~fwl10bBXuiT=oUtxT5wPyG6${dL3rMqdrqaL1kOtc z07rJw$cn^j4GTvdAA>(>ip0o{&tmIvvsQl(@&Z{lc4}em#Za7+!-y6Qk${iDLRL*3 zwJ?_bOE~{jZ6eL~sx>*q%&EY60==aFUc*p?Pp<$`3;}f}q6cqS1I&y)EjA{*8M1t; z#h)peXH%Lxnh-llHkqp348+Y!B^>*AExW|nh0K&)DY>2#$<)MC#(@Z~ zKpLGRP3>sLhkB77G~w{HjutLVwxHs-<*93hj;{bTAL7W$2Q&)x-fT&OH^ zK6eIqbRJVWo_MkngLrc#2M>}C|>4$ zr>Wl@^~*XG5QMd?{tevr_VyGR$+?qm6=E}^PfR+WG!ns;7msw=ft6+JcXjn-1a{=b z|Galxz|Mi^S*q zIs>c@44CJRCO6bY_%hL{CzHLXTUNCe@GPrZv$^Va=3@(ZqwJ4L4$O*-J{Yn8R?bnN za(b3?R$Df2slm!8h-@tnf0PC1h;uEv#>wY4@jSprNg^rk`+SO@lj6T51w6c3CbzMt z*K#2!avc{9Y?ZG_y-=r;G78-H64vIu=LM%O*+=1+%O;EkJcAPnW?{c@EJ5)+Y>vN* z#7$33*#UEzbj98nX?XXZvNBwL$xR!!fwQth-Oc+mv}anHx&uyjD6pyg!sfENreb!9 zm)k8Ip!X^%%i8=G^;|aKSPaQ+{b$+fUM213GT-gR0_R5jXkY%c~&UKB%*@*wg!12U6wXVmkDoqTJELG zVd=FbP(7BsUXMRxM;kp>dD9&pYPF-4Y*j8sy(E{t9z{)ZI<@xQA5ssWUI=1Dt5=R{ z6DKTke$IIJp4i9B$rNr|B%f3mgdZlNo;wzwNLvSPNeT#4!Yym#ME5M2qBI-(|d$bE`96Oyh$Y=}J1Q>*j;d=0jH-VvNS?$YD}qEa_hfr=WGiw0SYn>IR6-=(<6}E) z{ilNS|DwkwZjXhsUyvipQcBBC5hN+u`7Z^@^JwQU2$H=m|5A=*rsW7nC4XLySOWyf zuuH-J06}u1dj6sy5nBL3lIaPOQM5;$;+DJj5?;ux1Bfd}B&jP$n&e%5OQy(?gQrRX zN{$d9m=Yvzv6-tS;W2`Qbwwpn6JsH(D@gjlr6)+T5I_i$FP9^dj>BY_ySjp;{iCw70Lh%q}rph!knpolzKp^AApM9T)2E{*wo>DCwI z)u~eaHjwfEPoOa9_h85;Y$?IhURWjnN}gE#;;!;bn)ij)-7jSqtFj#>M)msV*O4)# z#am2mybE&QEEl6r0kC@6mCG+XDa)Q6DSiodGh#EAB&-OT8;v=y8FRLBz|&3^=>|hA zvA_szAX(xd|rD3J2eieyq9n~^w^#y^`6=O?z& z3Qoje*eSDzPYG2q_sd~vGiO|uUSF0K>dBUllOU1!6K=oOWk-aPY=~rXeeho82S+Yr0ZAWL%&3cFCtd*PUc_Cc7#T7 zwHjgV4N{n^Jfsv_-ftqf#ed#J_(JIx`)j3mO%*SZQd0a{X=J}C-=?@2*5iazZ>CVT z0mRcmM>rvII*35uoEZD!Ylx9J3|Qi59hr&9q-{3?NG5G{5lg<#bVD(-HfC=3miABD zj+vGvhc|(~(9lPyC9$Q zzBwm7rf0dw+tz=;G4C8^+^fEl^oc|kR(7x8i%Q`)tImbaa&&{8>qyyrM*OJjduLxP zv!ueF*VS>h?*3>w#DZN%#pY}E?3+T~`cR`k5))>^-@#Qk7!>QJx$uf1>a zZ>q}nU;E_IG)-SLEom)~LW=>Rgi@G-JW?o3^&(a*OvMN33yVS zDYOyZxQh0HQDz<5Nd9XUzxyTra+g#1r|5Sh8*pcUOlB~B1!&Ik=low+_?YX)|1;<} zM;;!)m1txc?mCYQY&XlZ^Z;H^I~c0&&kW41%N;8D04#`unSaUSb~=WL)@MR^ zMhf2=68ZxeTZ0fHMf_K$44IiaN%EL*1~V712&Sw|&I($1j)FAGH<=7Tr$(!L^l5Kn zmoD&c*ni}BH<`L?I1G%!)838*k;l8a+S`-aBoJ|0;D<{30t@M<%SjC!Hbo;(w>w#*yZr)QJt zFWykwC)u_h_$}CEL^c=Yg}Og+SWZ`SGKG1ezE50%Gho!yN_m^Tx6O17v`aS`(-7O^ z@RChXK3MRPGnUOZ(@QKzDe*5L1(lj{5Vs}QLV{P2Lm$lYxc54aA>Qt$8{yS< zd?`mXs*FvWnMJJJo_v8C#(6d3OyxM+lN}HZW-e9ZG=RY_we3T^){7BX$Jk0Afio<2!C(` zO};MjbLb&ow*y6*c>a0^t&D)DUxosc1n zi4)s-B_cr7rpU9b)WPZ-;UF@6PLoH8+MC89F;p_YMJUpSAhNIFbFR$Glgx4%DR-uU z7eY7ui#MlhR?wI#1Wt>KM>sdwY!Vf3rF*I#3@YVA)4ir8hz{MP-0@EbU;DrBe%r^i6y3motcn@k9eN}^L1&x?4Gqdz~_x%u75 zZipgZp(%L0ZY&}HXuiTahc9OTr#;ZZq2JYc;$05C5q`-k>I0KpQsiYsEx%-raT|FU zOnI9l(he-YX7XX1C6a$?a1d=Z6Sadt%RKzT$jB|wnD)VQA~GMzw;eo{4ULhi#>MhY zj?lTvx|02P+1K7k^R=NSbMAmjQ7W^Ye-azvV zY(Trm8tKAN;pby3`mRCm{IaNrH4k=+Q{N2@l++hZufX$HZ@y1S{r|gk(l}W_ZGbhyfxg+3%9!A*u>(PR@+z@F;cR+7*;jN_k=y5AJmesNc zVqCN**D!yT+9L}Q;WU1Y#N-0beeWRd|4b90fGMsBh*`#9F9v$cWP&TY@uZdPiRfY! zSrqzN#GWzdp>u=lTgq#beIU?fCO>p8bX_D6xCPmt_5re!boXD@CzSstv*>DBPw2?% z92cB(AtbcftW*Zl??S-?5bsGL%xQDMQyf_Er z#Zg(I_fAh&F3k!VzsOho%R&>s$epx2lhqeud$KSMpEdy3)>A%}>BTC|mz|O?%M3Mt zF*E&Z{O1ji;qvnD`cUSRm7!O^xL`E*TC5l^+PP)A#X%gOOad}|pBp&zuP?^U-G+dv z_xG@509}mnafz|3k}vyn!V*3>9`Ur#ISX`Z1pecTmo%)}Eh9tAza00o4^`hI!+d9R zYi~S=H{AZcc|h(9xXtC(2pt>$CDirhHL0&(i3Tr-&KpcEKWuvH%238vC0S3SX-hZu zin2|nM?s!SC~Gr4up(6b)zm39;18P|2*fQ|wkbIiWa`VG+fz5*rwRyySInR{jrx9R z=)hN#3g2CdDZbZo*z__81N)fvgDU;@pwnNnvDbOn^f*`^(wj_=Ee%}|%}oqK;Mp|l z=A|J|)DiH5OMRf{;9$XR8Jyc4(tS|6tUXNzyM$x>LGp z!fPb07QiWKt91DlJc9+xrCTN46>vjKzAkQxa^09{BgI5eI%J5v`gcU zzBtRw3N49=8ZgoSItfP0Bb`AhTzaEkILIOyAjF22|TW*<@y2{JWUfCFo?nE*@XN!-Am% z+hhq)w`l`TD4}jk*xw!x;AGwy=UbT9Wx&h)iFk%OhSfkn`yU$e4`_`W6SQr1Vc5tF zbV7yC12P$PU_`e0BhqA!!JZw})&Ez*G}ijmjG8*LwZlfx*D= z0KY?knfd_8z*g)N1}ftjY+=4E&Q~#C9Op}z&yVxRn77CIv&w#$4NoGuK`Dh$qwNCfdKIqv1{syEL4pVMW8=%HtL${=caNy{O?wHT-K0w`h2= z7VtU^H*5GF4a-u%xoxjnpl_d7@i%)^d{pDRHGahlEDsRP&;mRm1w@}}_<1eBXf41s z8g^*@u>@~;QB7c}hRZd4NQ*b9@pCm?tYHgdq)%UJiht5@EaK6cVvHZVPfh5gCjUmm zjr&!3Os{GfU!jGwX!u7hTnczjpA^j_7J*ZX;9d>y)bLvx&eE_&!xJ_9gBH&T4gW#I zZ)mtv!*^=9Zl93p6TSx^(=OC-mWB^$IsI7UV;Q_tOYrnC`I_g{cs6VJ2@U^A!{r)| z#k1vk*d6&-rM0zR!(|$dDdrC2JL-goL&k7_qsmj$Fsw!J?OIK)@h56DJ~Z?DF&ViU zcT~eaUMoUg>K`ByhJYOy-o|L^)W)pgG5 zHEWt`t5xOB|4Fa8x%yUTWBo18HK8d#Bo?gmH?MBI-sx{@a;|A=yk4d>$w`Dis%);W z$1N4n3v+%LW4w^5MYE)W_G*|DeK-GL4dZ0M?bwUB z{#}sk)(5J^s+`^}YASsiZq@SAvoZXT5INSa9zlKg2^*|ThJy>;HT>Zy$)EpO3a&qWvMEN~KifbG5-R_}?R7Gs4gMxV1Te!Ne+FzZI z2U!r=MnSyuEp)qUS7G!wfVVKeU@87#XhA9K_=d??*a8{y5z|{H7n;`d6>-kR40@jI z5*gvvVKNqW4wEf_tm}Lki>Uj2nM`9pWPy0D7t2I?&Qn}06?@Np4%R>S0nL^fGjtHS0JF zM^ru_1^EzQwo#iN3i&J4K*u1f$VQ+RE)OQD-GWYYB$CfZpCYn9gZ|4hwObXj%=jOI zaN3T^rm@TcnXda4v#k6%^xH<*VYR?BSV%mufeO_$iKm0DF+RUpj$~QMvW}Z&Y0?g0 z27$?>ux6d)sSl5=VK5c?v`G2|c{z^eb&?Aq$3YN)R*da~s1?C5ZGu6yrC7 zGJC&(qm)y&)kxN2Gg*twKK3()1T`eX$p9!R#~@N>VXDeD!f!stFa)lc&Q#XPSF-U@ zN^_CoX1}tS-(~Dqt_Tiy7c4<)_y!VLR?4chQC4vhWfhHxruC=v*?dXfJd;AELss(K zPqgd-A}_+{p*(tCpV01DLVG!(=auYiE=Zy*NS#-N=UD{)8;;2Qck~1;u6`XaGAxu) zoJtu>ZBc8##h2hUSLox028f! z4-r+&l5P>ZJ<{!{VY)-orG^`rDK$LFEpmkE;1BEyWwvA*oy~$KN%A^3$Ttn04`Ip@ zScf*4(#>m49B#6i-Jhg8<1&Rxlyc)%%YBOVI@|3j^@g-8>-(k=RrZps`HqsD`MK`V z52ekwPEHu_a5*MA=*1B>i_K$=0@2hb(B%Pa$-+@|Xe==u$i$^LW`d+t_* z>~A)*Hk*BUBaAenE{RN)R_gfx3)6>qWjWOHY(}!-Uz|)01}VD>@*vC69au-tpBIit zH~(uC!V%)7z_5tZ$0X61hLJQPA`S1a9t=Nb5qWZ0?uI~C)m{)vimo29tj7G_2Vu-nctsQ1uoEsDx$qBP zV_OZB-GT%gOc)6%WPRT3OP`!hlUpw$SIZ>IZmu7;L&|_LZPkd1l(mZ8pH=ACZEqO; zz2ZD+yTdT{R>kSF<>x!lH5<}wSvDJm`xAwIHW!1oV#QB}Gam3WfX(XLhK8Hf!>+sJ zPR>IiJBrb0E`)-T9>ks|?6!*BK(-rwq<9{U5yg}iP3gDwS$ovkK^?v=(8ttCt_>hs zIk2q|`x~G6*HowA@geM9DaW5oIdwM5sZ64rV(`sb{n<;>X->;>`}Bri*eCns*C@B>f)P+Gf6Vp(6)gxvx`XpVVY88Ew zrg9l`xo}=NRU|7+|5dun6d5LC5sMl!ja8iQcB6BhplfAQB%6-~;5GMBdQB-a7MLst z4Aypo?XHyPEE;)t&P=mqzH_H#yz@@WqOqQ|6njcmiX$ag#buW1irJU{>i9dAxv~CP zk!`n-z39BUluN~{C|AkbaNY;>1I(7R(E-^m80W_!P(+aTndWckUj}o^d27&Fczb_h zof)K`w$mQc22nN$<3yT`=|C>Hw7qWXD2^l1lxpG5yC`|2>78l?_Oe(sON; zi@uSIzLAT*ky|v>`W!AmeX)SYIp7@%jZ&0VJlySwTuh_bl$DoZq>M8~WI1s$dB@^z zWE{~61j$3B#jRO1u4N=$(BPo7qqE3zXr}c4C3McP4&q~dJ+HV3u=vz; z0|k`OZ|dPSQSD%3k-#|Oe+uwo@9KF33D}Kf_ve$P?jk8`g}j30NSJl>ygr;r5R~u+ zM!@JK8r_hN8FC_KM;GOYd`j(4_Stw9o{0J92qJTR1m`n=*)qa|nNBE6cxJGaCeo)y zyLl|q3NJ_!d82s7cRY#!d2k8=Ml0(-uJ9}u-jF0l%Efgte7sk{QF2HWn5aNh(#UA0 z&(0ljS{k-+`B<_Gkg@TUed2P;JX%5NhpwRHx}{{qR<3ilpf3THh=L2M=fxx=O%hXS zv~N^J-lb`DX-fuO++e3dJnS{*#H0W@j!s5og=DRpg2>#|F;md1fE`+7dS1-7(rnoM z3wup8N{phcsJ)*zk9o`Ow0ysrD!LNr(ybOMZB3+-7Av_MlE`^9n;iZejKB*p0&}S^ zNzlVK99m@ju)pX8;irqFO-3V9948Lm=6UOt{?1Esdz}S2a#6f^rj?Fxy9Gtv= z3ol6%BMl!Y?r@(?fHALE zXs?-TgR_D&{Z?Xs64?nI6Juhr# zvQ+}`=*rMR@dRe~OHebHVlLz`yh5BYO8ys+9n<1re?6}clb*{9Z~^6t-2cAQ<5{I_ z2WFMa;CKOM8}DEb80Rzc-pt+UpDALxoLjcB;_u*m0dS=7JS8{DopO`h>3)dA1dE)Q zmP!@wiL^}ahl6h_74X?$Txej84aWI$rc0C?R9Cx+ey1pn#^;qs4D=60xyIUQ``(a} zmG&QsgFYTnWIvac;|REvELuDQ``>il9@Hgb4JjRP0>ylK6UXmqjzHD1J=Nsn@Hth#e*RHmo!lj;BL_i-y7g5 zyT$up%vcqlF)M7%NKZB6Fr==kSbAj+&p^fX6P(m1p(~FZ9(m;O$a=Ev$9b;Mn2WnAvgU9r_L$dh1rO@=78da=WM9boSY@5~T6IdNW% zG0unomLbN=X>vCr`c$A*_kuVK(B{^J_}n_tE@tGo=c5+`jzNWU7w5yb*~Qq=oIB6c z$lWDyegV)H_s8Yp9)8&_3g^=AVKf8O!JFeakLO`tX(i`uqW_!|pE%=uczmW+0pce*gg)U=}V4#Hg2JBq1KhoLpM7EW( zv36!-?aaopBpb()?4r#6a4<{c{+we|kfjN=*e}%y`VO%%t^xRUe4-J(XKXyH`|1Ui zt`b<4B@Um;5;Lc9+0+uMSjaJsagND4bwrYn{rT6qO~vKAMfi0^WLg_wXYBy&-^~{y z0k(khFLrO?0NwdpDH?LpK=EY^|MYt9VS4NzyAjgsrD1pt|7rQcj`$gv43V*bnyFZs zuw?}{oh<7moTMS^)rFk_SvO}!r?(ukKIpN< zI`-@0=OKc!HU9UJivNA2;{W4EDqQNI_XI`DO0%B}+W;n`fRkuQ;w>y0CP)+f%x5K3sSoH_e+j zRrOsaV++q2C8jR8SDCl6YE6CP_5OyH_07#q&7O7sI*;e7(khOcqy2Z*Xj<=I*|ch9 zb9LkO^}{vxtsNyUG$hWUwe@S)Hr-qw4mrgg`7>#y_GH6KUv+c!+IoL|^UB7i)s3r~ zfST+5>zf#duMu0B;^I#xs5%IHzN?@Wga@Mk>Vao_`&i_kA;7q~J*}%Y!aq37o|? z%sB#w%nU65*yHe{x>shxXa$qC7x)~Y14dKI!V_UIa@-DUS10f|-0^k(7WZ?&YjUu~ zfZPW>CYNY8XeaR4(YPwYw!;m)5-Xp38{!0xS=N99Lnbpeo%I(G* zm?#Z2{uP0iU9Lvt1>Os5r;Ida^+Z2l>*oQ$9?*B= z8i@Yfga;y!>3-nj&DctzCi;Of*eMzS;{&niHgqB+TmjtpOFZO{YWD%R-%fN7=nmjt z-hsril@0>mc&AGHfUiX%j$!D=Pa1gk1Tw>2MEPhY#xMT{LliXQ6M!|K*^V5D8I~D4 z??=x>6)~Oz*a}*IUCh_W2LUZ1ScHH-1Ka~T3Vh3h_`fW)eGu3k#y9a$Ta4vq5|^@! z*8+;6&l}Ods&s&O7damw9~=O#1^6CB>46^u><5jjSNfwy^P%oDz&6OEz=;nN&4yi$ ztpUY=3M9n>p8$l3 zB5nb~MDz$z!V4(z6EM31+xFqv63}VDQ(l2SXuU;!fXbR^JD;2j5u<iJ@Bo0#R*bt2qU}L5u04cy1&s#+=>tr{H!KCRJF;IH1p z_(lhe0>AbLj3y-31N_`Wpc!8W=mV{{ z1?cSrdK&`U7Q%o%SlAdp3+M&S_@97dpc&Wv3C-~&&;KCy0A?eB9$cwhV0sKcmE8^(`ejAU#e9ZF)aP5~E9ms?a z`1`;&@IWu6f1+c3hoJ?!1333bmDbyPYaxr-g7tP{wjY0MfUOi}U&b#dATyvD&%{Qd z4>aQ=uwwUv))}uUz-I0HVK=tGgYmNf4`{vZm~F~i3ZR3`G42M;2FkRkaxRpk=-h8tJwzjA6q6T2Y<;VbNH}KeNFiOxI{J>m4U_Ua-_}{J-bOQ9^YXgFw z0b|>R90H#P@UX%2wshn=*ft<%dPkM~_rqeg@Wg6RW|#zVn`$Ax92jDON- zszq(p!RUanjPC_-!rOp9(&%HrdG%_X`M`l^HKrT5e3hWx2+Vj5fJ?ym0gc`Y{DDRv z1wI9kwE{e^LC{fzT@3sPfaTr5NvqXbv;lhn94F%uH=q!^a6M}W@hiY-1a1W`z7a=o z&<(&F0US{$@W%jd4#ti(D(wV50uz)Sh7UY&O`ifRpz+W(nE`2_@gM^20+gd4@a6Y0 zjpqCRJDRIt8^HX}&jCCGV4AP*+04W=W40tQEwM(k#fa@ZCCCtCw&Ac2*{yYN_U z(-y(OqYE1$VdZ?%2I` aH$MSq#~M`88SL!p?CXqlzEFyb>Hh*%Ain1S delta 21632 zcmeHvdt6l2`uAF!85|LqK^a0s9dt0fhJuRs!JxCL!O&1KZ>XqLW~75sse>UUZl} z@{dZ%pkossm6AaBL_aDyK?kmVH0&Oc-z5zYba|Lg9~1c*w=F5AI+*d~WiwXvK)92Y z{d#zDklD@#upUx$e|8?6ES>w!d8Q8JgQRk<-oz&jPanmZi32L+td(dUXP$1lsEA8A zD@MZBhqJ_ONxk}&Bas3|OK;9xR3P}+IBS=-B_(W^EMrHtgKI;eicLWYojOq0h*d?4 zJc}6n_&V@VkfnlNVCod9-$+s_Av>t6a>iVYHH;|Hi-`uJq09og6!jsge|I!rluA`! zMGlP?xlwh0lLqm`k(@27GDm|no_JkK2xqb{Dv-+A*%O~L)nO`?uHB|;rh(S47-w?j zw;Q5lHqQ`^;>hF?y;vBtM*5q zVroqI1Z!RZVC4nU^`5`Ztj-M|EXm<@4}?F(!wz0%YKmogm;$=a;tX2>y2p|nHWBo= z<@&hr%cSY&j~Uz4P#Te=oDA4Mzr<8$chA6yLpl2+_s;_+D;@O#JaT2AG!i{Z0_30c zF;m~@o})T#Ws*x#pR`>gUF3BKZ9{oN#2CihRWCZrF=#o;$Fgs$NcyGG>gAXuQ~wi; zsf(&()g1}L5^wsMTHjD=mMcFmjgXajPKR7^B&E@R_&8JlnvkA&spG#0tb)MsJ_+$N zAz*#AMOKan>@g%v3D&{Xg?+Y5$2Y4#_IW^Re6a4m#NFnwzz;CPfcPX&1U$xSGoOI` zd23Aeakv~%4jVo0yM7!)@$+lsqB+WG zSvfyViEWvxm}176?+0dmI3F|2t&EE>pTs6H{%eW5z{8?7vmG8-k$<*CT@pHwcpKO{nvxqQcB; zpsS&gDjbbuXN5Kb?!at9g^84&(8$#Gw5ehJ;EtoymAce13%T@?Royn$EiJUFe;+$t zZOL#+ZA|SsuA0c+adpx;rp_7vTwURWBB{=lIZxU?EnVv9p^llhs!#qK&^`1^+2?k; zzTTE`cUcR7uAQS_fvE8`%8rub!GnIJ|_+e{ca!h=RkiB>AxmmFZz_# zwdG4(vR0|z&wR}K!6!mvZVEMZ*{kZ#Sx(8ZUwv!VSm_h9dU4hrVdMYA)DgMyVJ9)g z%4Zi1EM7|`f?dWm57?JtTB6Ij4pP!_fT=CHiKf@~F;@3=!EKV7H*c!9FBVs~GQ*Oi z+-Q?4k6AZGf{~-d%gQYlx$=~C)8%jyxWSfF(ZQ{iqp=jq${3(Fw=%(9wLVX-npq%M z%`djjKI5*cEXUxc1#%`dTPq(R2e;bdtg}CttFludTMpUHR$F!+X2VXaEw|7bd<+#2I+(Ye*&w(8s(aNBim zEx1@EQ7pG%nKkx=H72(KD6Q+z)B@CjB&y79J#Pk(%GHPF_x5shRQ8E_Wnb=Q|3kg0 z<2$GU?y8x!Ze>@!NY1Z^6Hx=?s@!H-5mmVrmj&%^$x#|*WoET(efe-s*kM_jzZunx zJa3(o4*s7{Q=y6g-!OTa&-52$SWZ~8+ZXu}B~{0QUaD(BpQSqk_P5@L-OgvCnDe}~ zU<_7|F=r>4(frvBKU%%nG)y+Pud*545G!Xc^ABCD7qOX~vT{8dv(8x}S7kZNWQF7_ zh1V;&7P-n-O|`wfAl_R79Xfxe&Hu?p-B(z^K8&hAORJDvvB3hXFBYrLWHL8;StwVP z!rT|+rpj`tLLpzkjJs$Z#!i@9_%@ct zE++PF_Is|m!p5?)gKP{pB^wKKRpqvzk*D2?@Gf^%ZYxY~K?R1%{cc}k%zrRDMkHkS z4$?_>cLnW!lP}WZR;Lvv@M^WBuunK_{?JLLZY%8N?SZk@?VJS*`evfRx#j-taZsUT z)4=WnjR%EC!n_z)Wqyi#LZjzqyc<2eW$_NiYuGD4Dqp-L)2!Y456s$}38y`94Y#u- zV2?csZIz#Iu+nR+=_9-&;3*|-X1Qvem)vHav&>DNobN2bv^QC|ZbJ`pee$nZOKRn? zlko4{&2m^SHKkluBAex^$X^5&IP{-k_b(QKR~(Z5n*f{w%Dokb!mgHpEltfZWY#b7 zSdvXS6MnXCU5zg0Oz2oY=sz_nSitm)9B6{ny3Zh09kBO%500IK2?39U>ky?yDA#6< zDRt-`ldr$(Q|Zluk?_$RGp0|ja>DJyWb*A0OrLm@_+R)j5WpHM#`W!YY2vt1#F(iW z!tG_XYVlw`R;^n+z#9uNM|nTT`f>{lo8iA>K3u8mRL=091`Dm6&#(g6(w#@MoGDyx z%5u`G1-HwUX*rgaSr*ynlfm;p*URXZtmI)KzX3xn+c5nCr6#eAgX31h-O5cCcZ|{t zMs^mueY>cryK0`w?})q64f|)ftK801cU8VK1!i@TOE^$lhFi&&{Y6En0&{{21?&%C zZ-z{s$lQYrpxlgredP(JKD%VF_3)Em`QJo*f}Xlhm%YNJ#pl$m#lMrDKVSDzv5QM< z|D|4Y`(4t}bL!^X)1;(N)j!|9CcgJvUDv*V{XI-pvTeYAq>0ttbjK9V*VZX3@P@qk ztorbs{iHw-_2oO~#xGtJDsmA!F}=vQN0>Tbdvr=NLvrP+CDPDO)DKqWC;Vwkr)2NLA)V2O zSzYGp^_H-Q8klStS*O3C$V zrY|<_&o!jfXC8~?Os)!(F_}?nZX=8V?nvS`W!lgjpqcO^b9BAB*>|nG23I`<*xuy+sW{fiIS59`Su0eihj=}nR= zzJbGa57_q{P(Odj*>}BxrBG|>(FJmB0sGad3g*e62*ZxKttp5CcPT~eJj>eoL} zgE@AD7XLHkgyNI{9j;rT#UU%7`qLiLk1?`>cG&CcKCOQK+HnKG43>Lm4a!|;QDK8refUMg5UFi%tdssdzwWj-C3CO48v_BX5Y1S%<1x8i{=6gVhC{2> zMa^5JZ8hq*&38$T8g;?3rPkXvVLioZESe4kC$^}kj=81D9jfDU8mvr3|LR>mhY!oRB0a-dQ+W>wS0%)N^vu-ee5g~w9X8@xOz#al< zH3Q%jfcK>Ta%He5`8v|eDDsaGP@n^n0lX6e=IOxCFf(=Jl@L&<1MLQ|GXxarz}E)g z^@fmQ9r>GqEDHfEbl@Wc$PED{I&i`O(nEkp2O15aPY5W}f!7T{s`Li?zD`G8F_5#H zf@$Ft5U@!Ho-%-gAz-r({N4bb3;|npz-IuJ0KB>nRqMzm16dI&@`w)HX8^Z^ zfSo$9$^gcPfG2c-W*(U&DFoE$K%oIx1mGpt-m4?C4dnZZVBh!Yz%&E+I0U??1K9@f z=MYe<0}~BkZwNS`149i!F#xZxvR+3L3}jWP$eTLQ-2mo5z{ds<9Rj>3bmT7v^6&eDeQ(i$0|xL#2>4J3_83582so_++YMlU2>3z=stjOT z2x#SAF&E1Xq%4G-(Tl7#fQ2EzuLCO$;D!+JjSeg|fZ-vaT?ghFKx_#334j+X!7T=I zX=AYO9eR-)4dAN~a9Ia34d7S^pld&<$}oW15Ma`Q>kL2(0Tuyx?Sl+tLkNk~i}W>s zr6Itk1F;6+4gs+`(8B=ILO{F@NCuD)0upuLXYc*G5AhJ<(2;fn@RtV%#i;{d8^FmB zkfH;BGk{k^z%U*7$N+YQfK(kgaX<5lNxCV7r0d8VMv*&0K!y(NH-H%-V4@B@Y5?Ox zfJ+B98-ODOWa+?K0=!rvx`hx~N0u2yz6-66ZXK9u03U^bJRQg~fHy-xz7C`rz@895 z?^t-s03Go9w}g;+I?~%Ha(4(Q)B&3T6oh~x9q48NQ$j$o4*a%Jx9*@2utEp^Z2(bT zAg|F?xqr({w2O8SI3G$5vuE4rJ?`DTqBzN9Jq8FaxlgI0Oag#p^_{qqv(2 zEA6zSxTvd}5LEZybyT7-&A$zXlHn-gCM9m=qVo6L$}=10U{jW(bjwa_ka3j41zh#% zFsdkAR;vBzi6V9-*4f9-%hn;ybQN|grwf?JsTVp(H+!c7jq;M?ITKE-+n=Qb-hcw; z%~F1ZzD&C5eJ@bm4)x2d*BxCGmKkWu3>?ke_#1V2-50@=dO_}zi`kU{Pg2^kfW7Cv zj3KkZO1rziz>8@73dHf+S>UenI17Wtis>*`c9So5P3CK;g5@el>EK!H_kQ0UyHMj4 zS&pldYEu0AJ&XxCe~HkVUB1kEm-$FJV{|RF`8C?F`lsiM-5CZe=WXQiv`?}2&!X{( zb}JjHE41BvFEh}x%sS%ey^?2^|E?l}(a7^ttudx{?EkWy#nY-LB87u(PS|_PtScLT zwEhd@4FEPbVN`KtkmHNIL20zU+z1?3Q(f_#z|lbC#!F(D9*}XOL0_|X1`{<_+j?Z? z5qEyGHD-Ow`79jDn-Rxkx5_@$kBgyfzcB(k-8lJ_VaHq=g?_n}^)<3my0dK<4a#$6 z`d%8hFz@+qh2O&AsH(PZa5{|K&Rpd;cU5VP@Kl-3Rx(|hLHO5d?AOtaGSs49Z~j1@ z5wpJJN3YduT#a3C!JfJXqbSA)id!gAh{S0!6#98D!m*N(mob z2%%qSN4(sW?FirjGevK$XnK9CkI||VcTo)GL z!i-TXsRdQJ#dH`|W;*mc(nhyZ8t>t=1909ZI=Uwa>7pKOgOggYg_=N1nleR5QotU2YM6(7Jl|I=%86rkKQ-a zT}Up5Ph1V=POm)Xv32&p-oKFUynW<{GwGJ5-oIhJxa85DvIRp-dm8kJw9&XLxeYRI zPx4Ts`TRO@y&C7uZ4u6f&Nlnk4H6@oTOY8$UrI@=0w|WfIcGvsX#}oVs>uO8eKSS> zbCsqXxZF+mbNpw^c~=f-vR;k z(7UnTDd_6iw~<48cKV;G5hH<*3zs~)la!;=7KG8~@LlA*bwP>Ent~yp;{tPz`On16 z#NC~fvJJq26&&*aWg#3prGp@o6N@pHqW9K@FSAZR+Gas{|0}>nIjc?tApZ9zkwbL0r%V|r)O+nx93RT;ObGw6`scfxiuBr2Un9aQvH&Ykq!K~t@ z=f0;Hdk{fIcmd%cLL*uns6pj`q%Ah#+e_aiWYU~ zKi5lpKU4Srvs7CDCpGC`_YJ5nMoq3w0Z$KjM_IsLPI+1J(Mm@j(`>&}r?;g>yfq$d ze?4Ng_^ZkM_v+ZMYk9V)>YC+hH~(&)C90aWOl|i6AoXfd-#EKo+Vr71{#>c_ZLRvs zxeWt;UDBm0T+V+>c~sS4RJDGpsA}jKswx?=dg)&i_%^kvt(Ln)RSzy!wQqJKL%lk> zNK{pET)pW#RQ0~7YW|C&syi26S=D^Xcd6>zNmP~P@sU*3j*)e}+mBigzTJZlkCv+M z(!T4>pLNvt!vN5SS9%9lG~oMaYkQ9JZrfnmC)sawX$KPcbJ7oewfTK`qEsBCc{nw8pxK73-tJK5r$HtAIM!HJ1H5GG|lr%iJ3TcpHiwCxVQ zQd;-GK`Du+;ktN_Hg*7CDz$FXo&-(qo*}w&!4r-4(_1oNzc5zsO)(ZMve5%^TFO8^ zcJgq1#>Zu7)THHD5#H=YU#mwViLceuJ7`sliX>n;?Hqhj@su75p|j}#{XcZqb`9hi z(sN6-PX_WA9mit-XZL009i#hm$7p+;e2)|!t9O6rXl+?CPmz`_)w_SGNP8ifUn?Ch z(mqV)E2I+>4o)4!)iCL~UD_vudA_vqac%ezbo%GY_LQc$?^lx)pD9K7MxYZCXozNpG! zXaw6G2w-NJchkAU`*FS)Tj;~!-G(Lw?8gTp zgI{IfcPMr+pAgJfLrjksBl8@Xs3>6nJ@I|^gL)pnjS1K*h=ntPof=@1t0p~$I^-(z z;~;JQaN&jXOxm$b-*g_ZJC~Dim02bch&$7c`S+rIn;P(m&llAfR`q9Tugu`sW?0Z{ zW9(u?t0#yyMP0sy${zEC2kbe-p}`D|4V0V#9XyA`?2 zdJ|-1CAL{H|1=%>Dw7S>jsqkNA_tQ1ggWEO-#nf8=zUMfn01*Hh`df#QKBxP4FP*H zRj)*SM;$;90rwycomc+mFqBk2ok5EFW+{hYm3K(k7geX1am+=b<+`P7h0b$HLg(Yi zz$VHE0edZRd{NuFlz2BuTa?f5|HBj=1rKFYNmI)`zq_?cTb z4Fsd|Sm}*T<`2L^2K}?Ym?XMc3F1!>UGqfL2Y{XZR2NvH3rr+|AF=E=nRlQxMv1jq zT(5GfE77s#xW5@JzZLDJ`CO5~JpHg8^e0GK?Q9`FpM)@n=X`^{p%r-cIp-GUx50R#De9{@b zOE5i4m&5F9Z}(ThDDcsgS@@x!6OQ^GU_Uh!o<~c?$CSTbiSw|3 z@$}Ya&BF%cf!^AS^Y{(kXiQ7F3TOL;7QaPL`~LuX;>2M@M*Edd>3*dbeg`F9X%Z4R zS;teh`StE5M!q2y&qCti{A+;&X1sO{1}I>k2Ob8p+iUPs_U>TxddrRI_4CMy>!l1> z z>coj!-bkPE?*^|?P;Kc#K5%j*idr9OAl4iM$`KU5=dn;P9I~(P@2-U+-y|1ETTlAQ z0c?Z~!oGOzD0E!cM~7ptf>)4%D;h2GgukykOU*MrgKLgLKF>Nf8wO8qpAxXjP3HdD zno`I!;*#LwnH9}EtI8aQWLvi^#T@0~!@SpnRncr|GLJ$cj39F=QkfST9iEJeW=E5`H?r_y0M{T@ z52=={fv2-3oGyiPreNiR(q`l5JsQTEjT=Pe80Y z8nE^JC6>Bovg1ZLko=y=w2N2*E0BCbq~8+pJ`q1A;x|S7n~41)9wOw2i8xNgrv?9y zB0W#U7X&9&q{oYRqKKDKyd=3yWNs93k%%7>af*n)5`}d|BLpX2#M=dborq@%PJ>9F z5^=2HyoxlI_Y>NCi+OtSz2H$!e`IA#HOBO8{zCmkVC3-=i}h0f{38ps+n4a6-at=!alu&nLrcnh<6I#7wxm+g0}Db`({<CDrYnzLs*^R!5u^>pS-NMBKBz8o*xD-;NYKg~BP^ljO974D+^rl0~vG!SE9J8J_E zW066)sS}qefAO1;z5>dh>D2EnCdP-IjB?5!>CB74JkXh^kvh9Km{+z?1)>4f;GqF6 zgl`Js#^8Id%BP3&-Yud5&LBhC;wJ@JSLTJaQ^A8t=+xPX`?k=%^w~~@TWTpE(&yXK zO!-)pb;r3^q9@y@XvM3!Yzyd*H%PVC)%^L9ef5SF8XXOmHiSZ?HPh(8 zRPA*SUnse_*5Kjx9_>AKKI4m>E-tsbaSvty|l;H;Lo0FVUd?VY^pBNi*fw6fAmo754wu7;) z2>6?}p{(^1K>WR1`>!Zr!&jBhu@ZiR-6X6G=VL4q+ZABDC2i0>I3_#3#g`S&(lU?p zPBDq{fVm(R&4$3m9Tx#l)Wy|qrguS@v(N1Ov&LlzkO2J!+mv;p(M<^_~ibu@Ct9aN7v*Igx#9}3a zhlNwGI7d`cc7>o@R8oE=PsnY)f|n=cs;}e;xt&+=@`YRtc%cq^3-X0vEtoCc&=}({ zIfH)|s^^{g)M6}#jDFi?@KzD;3HSzWfDGQ6yHoT5KLGB$ zRWDNu-cBpZloff!->u5{Q&RXSG036dVT{pV(vrU#PBw~o&NR-j)erL0h?fCg7iUZ- z9=*{1&^;tKhIsViC;%Cu|II6-3BHppHWIv(RK!@O)j-wXeONRnCa2osA) z4`&gh!b=GI0_@^{B5ET%x)3aBJmy(%6UKNdR&6lY07Irb5|UHu9WJr(QBo>li9sxx zFj?YM?tYZ-P4#rQh>A{p$F!#PykGBKjD_PT=2?(Ug5K6pUi)o5zouuGI-h_@;nM@O zigKPLI?)`g=zXdvTIxL|hK<;q!3LL&W68ziSyIskmXyApCGj<^S9?sWwIQk`G7t0m zd;EG}GQ-x; zp9Mcj19g}fUrN8|B*XP@Nf%rnM(9f>ir+O;{A)PH;~!$TN+Z}``%C4lpL9D^n&ySO zEGdywqnC!+=9_yii_Np;NO9}pM>*^Dm zz!E;?L&7wK(a#E#%n~!Fn;2bt-HVV9D2V~MnrFz_AT#7DIfEC5!rJ5qdB5wK14cwR zL&{=JEY{;>mZG78R}cPAq7q%EYrac4?fVCLQePUr70K|6QbhM7m?G%o-_z0!7%px2 zW`2#8L~4K(K7&X~&~C2c!@bpiWb6Zkec&aG*Sd zQr4N*ZrsBA*{J2szzW|*B*ioBi!Izah~@>mma`!o^P)QveGtMUZ0L5TX)<}qO6vO{ zg^%h_vIT$_qeu!Z7K=zWTVhF)Gzgp1MujIOlnh{O43^&&z^ECWdF?laPflpKpRu12 z8VtNMug%)Zr}eQdU~C5hOVP`B$!oXxc+wDR#UlXU3ASPvGU2-MT zek#>EDis67B^n!6D&~9ZKcyNVUz(PNk;Ip#spL9I-NkfE`IDKgm!x9To6<93>?29K zHR7?zV`0&8wr639FM5uOy+@iKmk^&K#W8o!NNkj1V#8T%QFj)bZeg*bY+{i9n!B6rj&38E0ZFbG5<^S=^M^chcFby=%F|nxihc5wUzfvbFE{{z8-r+1MiyW zaot&*$IRl2Oe_w6lplv5D95!jO?wze8){%A{`^2n-ai^FPG0d<5S$JPG)JC4JbHGI*bc%4qj)=YzyTV;jNQSLAxV zkr;|#xL&u;h34pXJ|>|Zu9Z2Iv563&7I)^gIF%3RM=jn9EafFcZ;&2KPko`DMs1#@ z@)4758_+9+W<$pnd3w9#N3cAsT5&0zYgK1nJFM~>gg2I5*JY{Mc%9bg5uPO80{Z}n zA$#RI?WRYtevt3BlQi}4XC!*y%2n^-q3}=-S+p;VU6wZXH;o7%#2%2!n0Tp5SS#Yf zF){G#v!XC8USLdQY$!q^)X*q&&5K3O*uQZJGdz`XrseJ6@g|bgig)nDo3~(?poTsA z19mFsu`{3@3+O3?GhOn=?9S_Mg6#+3s~pOq5A6D%uyH|1Kn_$pzJrgmb89%W79pgY znYCw+5~_}i6KRiHEO_)3kL23hYq>pQGd99vsT>Cv6(c5LBmaHfp!a_1C} z`>0m%EM0gm@LGtM-dQ#hH-Rm<<}L2RI|p8B3}^M7JTJQpW;vKCT`K-*BF+=wHFpZK zxJjHSC(W7@lvT?1Y{M!cy0INJYBOrfvDgKo~S+h1RrrQeit7UrpYyYop!E<_v;qN zSKqzL!wx1r&FgrN<;=5m)t$GkT;y4L@Ppml7Dngk?poxtd}7aFCa(!E?{@X#V?s+~S>wx3dE(m0oyz;5}OndJ8@$ zw&FNf2>KgBGtvy-v9V`lsd&M3Jsr_0X%51hIP4Hzh|q#`*OL~VxSmIN1bm_|BkV)E zHXgTOXOVFu(++yMKVt{5Rc4OPlNFt?f~O~yc)T^7NE2OxkdAcc#-ED= zK}Rxv-hxd((HrpM-HCMP3m4ay^d>$%7#B$B8PP%51;iuWxw)kMD)!q5#n@O9Z9`ax zbmtC%i=D&g!_XPL`H7xGSb=mqXu@?EUFdXxmX1X4(aCbqkGyGmrWN$|(Wn4G3FxY^ zP=o^IpnNG^%w0Ht`TIfO7momMqZIdt(cwrMAi{L`K1#~Gs^qa6JF98j}{#h6y zb|e08nj0NJ6_(bQP3jbdU`>|VQ;9Utm(TD=Z zUPPK`IzkK5F3=`~y*Tl;fR@jJS7)!O9iFh%KA}R9_VU>dNhJABK9K8Bc}s> z7oi>cTjqNii^8eT1E9-^Pn-m=!Ja4`<{+9_hFu-f4$vv*0b99p@2hU?!p_ z>;Mmt$-3+YX;*kdIeZPnPITZ#ToHJYNk!to{n*u_K%b9$sJB0&NDL{_B?yT~6Fr464Cz)% z@4~hlL*02{h7IPEjJ<-84?fYe2*pUZgC2Mq>jcsbpi_2Z=#h4Tj@yHojf!2MAMC|3 zZYVzEwIXrb^XMh)R04YWCG-X(!Ttcd|4~nOy@2R?Innjvf-W?cy^iwm1ENFT!Ze0X zqR%&A`XNnp;$fWVkk;=hV1K%}pbLxvEoc$aL|q7vkPv9%8%B=k9)zEeZUOxof$p8?+eAJBEqX;>B+|ddT|ee4(Rm0YMD$Gr z5+XYAJ3Z|LU5Y?5#h}LxKJGZ#S%d!#=iR2jZYvZ>g(QySkU89rEw%81@+q)R~YUm)35lqLEb0yT*J2>T+?sHA`{LLhz- z=t~Hsvles$HusdC_!H*;hsd0OBBC7zKs4#!x{woefsrl*U5h|U%R!Ihd(clPdm3~w zz6%k5BB+KyBl-kr3j$R}baw~*X(SHlp!T2P_#?6Y$0KnBAs#>nXs=&%MX8`G5J(Zx zeMXvS{$(t;kjLj=wiqD;X?y}^4Gdf6_oPg z?P0^U=mUI$sgSevv`6lBu<>4JsEpkdkXdx?kV0= zyQh9n!=BbX{ypt`EPHKx +#include +#include + +double Gaussian(double x, double* p) { + return p[0] * exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); +} + +double* GaussianJacobian(double x, double* p) { + double* resJ = new double[3]; + resJ[0] = -Gaussian(x, p) / p[0]; + resJ[1] = -(x - p[1]) * Gaussian(x, p) / (p[2] * p[2]); + resJ[2] = -(x - p[1]) * (x - p[1]) * Gaussian(x, p) / (p[2] * p[2] * p[2]); + return resJ; +} + +class GaussFit { +public: + GaussFit(); + ~GaussFit(); + +public: + void addData(double x, double y); + double* fit(); + +public: + double* parma = new double[3]; + std::vector data; +}; + +GaussFit::GaussFit() {} +GaussFit::~GaussFit() {} + +void GaussFit::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); } + +double* GaussFit::fit() { + int n = data.size(); + double x, y, sigma; + + Eigen::VectorXd mY(n); + Eigen::MatrixXd mX(n, 3); + Eigen::MatrixXd mB(3, 1); + + for (int i = 0; i < n; i++) { + Eigen::Vector2d& point = data.at(i); + x = point(0), y = point(1); + mY(i, 0) = log(y); + mX(i, 0) = 1, mX(i, 1) = x, mX(i, 2) = x * x; + } + + mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY; + sigma = -1 / mB(2, 0); + parma[1] = mB(1, 0) * sigma / 2; + parma[0] = exp(mB(0, 0) + parma[1] * parma[1] / sigma); + parma[2] = sqrt(sigma); + + LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); + LM.data = data; + parma = LM.solve(); + + return parma; +} + +#endif diff --git a/include/GaussNewton.h b/include/GaussNewton.h index ce17b48..222e7a3 100644 --- a/include/GaussNewton.h +++ b/include/GaussNewton.h @@ -3,9 +3,8 @@ #ifndef GNM_h #define GNM_h -#include "utils.h" - #include +#include class GaussNewton { public: @@ -70,8 +69,7 @@ void GaussNewton::calmH_vG() { } double* GaussNewton::solve() { - Eigen::VectorXd vX; - vX.resize(L); + Eigen::VectorXd vX(L); for (int k = 0; k < maxIter; k++) { calmJ_vF(); @@ -86,4 +84,4 @@ double* GaussNewton::solve() { return parma; } -#endif \ No newline at end of file +#endif diff --git a/include/LevenbergMarquardt.h b/include/LevenbergMarquardt.h index 25bbbb8..660a1dc 100644 --- a/include/LevenbergMarquardt.h +++ b/include/LevenbergMarquardt.h @@ -3,13 +3,10 @@ #ifndef LMM_h #define LMM_h -#include "utils.h" - #include #include #include - class LevenbergMarquardt { public: LevenbergMarquardt(int L_, double* parma_, double (*Func_)(double, double*), @@ -21,11 +18,11 @@ public: double* solve(); public: - double mu = 1., eps = 1e-5; + double mu = 1., eps = 1e-6; double* parma; double (*Func)(double, double*); double* (*Gunc)(double, double*); - int L, type_, maxIter = 20; + int L, type_, maxIter = 100; std::vector data; private: @@ -75,27 +72,27 @@ void LevenbergMarquardt::calmH_vG() { } double LevenbergMarquardt::calMse(double* parma_) { - int n = data.size(); double x, y, mse = 0; - for (int i = 0; i < n; i++) { + for (int i = 0; i < data.size(); i++) { Eigen::Vector2d& point = data.at(i); x = point(0), y = point(1); mse += pow(y - (*Func)(x, parma_), 2); } - return mse / n; + return mse; } double* LevenbergMarquardt::solve() { double v = 2, cost; - double* parmaNew; + double* parmaNew = new double[L]; - Eigen::VectorXd vX; - Eigen::MatrixXd mT, mD = Eigen::MatrixXd::Zero(); - vX.resize(L); - mD.resize(L, L); - for (int i = 0; i < L; i++) mD(i, i) = 1; + Eigen::VectorXd vX(L); + Eigen::MatrixXd mT, mD(L, L); + for (int i = 0; i < L; i++) { + for (int j = 0; j < L; j++) mD(i, j) = 0; + mD(i, i) = 1; + } for (int k = 0; k < maxIter; k++) { calmJ_vF(); @@ -111,11 +108,19 @@ double* LevenbergMarquardt::solve() { if (vX.norm() <= eps) return parma; for (int i = 0; i < L; i++) parmaNew[i] = parma[i] + vX[i]; - cost = calMse(parma) - calMse(parmaNew); - cost /= vX.transpose() * (mu * vX + vG); + cost = (calMse(parma) - calMse(parmaNew)) / (vX.transpose() * (mu * vX + vG)); + + if (cost > 0) { + for (int i = 0; i < L; i++) parma[i] = parmaNew[i]; + mu = mu * std::max(1. / 3, 1 - pow(2 * cost - 1, 3)); + v = 2; + } else { + mu *= v; + v *= 2; + } } return parma; } -#endif \ No newline at end of file +#endif diff --git a/include/getADC.h b/include/getADC.h index 80d778a..2633e83 100644 --- a/include/getADC.h +++ b/include/getADC.h @@ -1,192 +1,18 @@ -#include "LevenbergMarquardt.h" -#include "utils.h" -#include -#include -#include +#include "GaussFit.h" #include -#include -#include #include -using namespace std; -using Eigen::MatrixXd; -using Eigen::VectorXd; - -void getADC(const char *fin) { - // read file - TFile *fRun = new TFile(fin); - TTree *t = (TTree *)fRun->Get("Tree1"); - - // data numbers - Long64_t ntot = t->GetEntriesFast(); - - // five X-4 block - UInt_t bArray[16]; - UInt_t cArray[16]; - UInt_t dArray[16]; - UInt_t eArray[16]; - UInt_t fArray[16]; - - // calibration coefficient - Float_t CalC[5][8]; - Float_t Calk1[5][8]; - Float_t Calk2[5][8]; - - // read adc data - t->SetBranchAddress("adc0ch0", &bArray[0]); - t->SetBranchAddress("adc0ch1", &bArray[1]); - t->SetBranchAddress("adc0ch2", &bArray[2]); - t->SetBranchAddress("adc0ch3", &bArray[3]); - t->SetBranchAddress("adc0ch4", &bArray[4]); - t->SetBranchAddress("adc0ch5", &bArray[5]); - t->SetBranchAddress("adc0ch6", &bArray[6]); - t->SetBranchAddress("adc0ch7", &bArray[7]); - t->SetBranchAddress("adc0ch8", &bArray[8]); - t->SetBranchAddress("adc0ch9", &bArray[9]); - t->SetBranchAddress("adc0ch10", &bArray[10]); - t->SetBranchAddress("adc0ch11", &bArray[11]); - t->SetBranchAddress("adc0ch12", &bArray[12]); - t->SetBranchAddress("adc0ch13", &bArray[13]); - t->SetBranchAddress("adc0ch14", &bArray[14]); - t->SetBranchAddress("adc0ch15", &bArray[15]); - - t->SetBranchAddress("adc0ch16", &cArray[0]); - t->SetBranchAddress("adc0ch17", &cArray[1]); - t->SetBranchAddress("adc0ch18", &cArray[2]); - t->SetBranchAddress("adc0ch19", &cArray[3]); - t->SetBranchAddress("adc0ch20", &cArray[4]); - t->SetBranchAddress("adc0ch21", &cArray[5]); - t->SetBranchAddress("adc0ch22", &cArray[6]); - t->SetBranchAddress("adc0ch23", &cArray[7]); - t->SetBranchAddress("adc0ch24", &cArray[8]); - t->SetBranchAddress("adc0ch25", &cArray[9]); - t->SetBranchAddress("adc0ch26", &cArray[10]); - t->SetBranchAddress("adc0ch27", &cArray[11]); - t->SetBranchAddress("adc0ch28", &cArray[12]); - t->SetBranchAddress("adc0ch29", &cArray[13]); - t->SetBranchAddress("adc0ch30", &cArray[14]); - t->SetBranchAddress("adc0ch31", &cArray[15]); - - t->SetBranchAddress("adc1ch0", &dArray[0]); - t->SetBranchAddress("adc1ch1", &dArray[1]); - t->SetBranchAddress("adc1ch2", &dArray[2]); - t->SetBranchAddress("adc1ch3", &dArray[3]); - t->SetBranchAddress("adc1ch4", &dArray[4]); - t->SetBranchAddress("adc1ch5", &dArray[5]); - t->SetBranchAddress("adc1ch6", &dArray[6]); - t->SetBranchAddress("adc1ch7", &dArray[7]); - t->SetBranchAddress("adc1ch8", &dArray[8]); - t->SetBranchAddress("adc1ch9", &dArray[9]); - t->SetBranchAddress("adc1ch10", &dArray[10]); - t->SetBranchAddress("adc1ch11", &dArray[11]); - t->SetBranchAddress("adc1ch12", &dArray[12]); - t->SetBranchAddress("adc1ch13", &dArray[13]); - t->SetBranchAddress("adc1ch14", &dArray[14]); - t->SetBranchAddress("adc1ch15", &dArray[15]); - - t->SetBranchAddress("adc1ch16", &eArray[0]); - t->SetBranchAddress("adc1ch17", &eArray[1]); - t->SetBranchAddress("adc1ch18", &eArray[2]); - t->SetBranchAddress("adc1ch19", &eArray[3]); - t->SetBranchAddress("adc1ch20", &eArray[4]); - t->SetBranchAddress("adc1ch21", &eArray[5]); - t->SetBranchAddress("adc1ch22", &eArray[6]); - t->SetBranchAddress("adc1ch23", &eArray[7]); - t->SetBranchAddress("adc1ch24", &eArray[8]); - t->SetBranchAddress("adc1ch25", &eArray[9]); - t->SetBranchAddress("adc1ch26", &eArray[10]); - t->SetBranchAddress("adc1ch27", &eArray[11]); - t->SetBranchAddress("adc1ch28", &eArray[12]); - t->SetBranchAddress("adc1ch29", &eArray[13]); - t->SetBranchAddress("adc1ch30", &eArray[14]); - t->SetBranchAddress("adc1ch31", &eArray[15]); - - t->SetBranchAddress("adc2ch0", &fArray[0]); - t->SetBranchAddress("adc2ch1", &fArray[1]); - t->SetBranchAddress("adc2ch2", &fArray[2]); - t->SetBranchAddress("adc2ch3", &fArray[3]); - t->SetBranchAddress("adc2ch4", &fArray[4]); - t->SetBranchAddress("adc2ch5", &fArray[5]); - t->SetBranchAddress("adc2ch6", &fArray[6]); - t->SetBranchAddress("adc2ch7", &fArray[7]); - t->SetBranchAddress("adc2ch8", &fArray[8]); - t->SetBranchAddress("adc2ch9", &fArray[9]); - t->SetBranchAddress("adc2ch10", &fArray[10]); - t->SetBranchAddress("adc2ch11", &fArray[11]); - t->SetBranchAddress("adc2ch12", &fArray[12]); - t->SetBranchAddress("adc2ch13", &fArray[13]); - t->SetBranchAddress("adc2ch14", &fArray[14]); - t->SetBranchAddress("adc2ch15", &fArray[15]); - - TH1F *Left = new TH1F("Left", "Left", 300, 0, 300); - TH1F *Right = new TH1F("Right", "Right", 300, 0, 300); - - for (int i = 0; i < ntot; i++) { - t->GetEntry(i); - Left->Fill(bArray[0]); - Right->Fill(bArray[8]); - } - - // use matrix and least square method to get coefficient - // SigmaClip *clip = new SigmaClip(5); - int cnt = 0; - double vX[300]; - for (int i = 2; i < 300; i++) - if (Left->GetBinContent(i) > 0) { - vX[cnt] = i; - cnt += 1; - } - - VectorXd mY(cnt); - MatrixXd mX(cnt, 3); - MatrixXd mB(3, 1); - double n; +double getADC(TH1F *hist) { + int n, cnt = 0; double *parma = new double[3]; - double sigma, b0, b1, b2; - - cnt = 0; - for (int i = 2; i < 300; i++) { - n = Left->GetBinContent(i); + GaussFit *GF = new GaussFit(); + for (int k = 0; k < 300; k++) { + n = hist->GetBinContent(k); if (n == 0) continue; - mY(cnt, 0) = log(n); - mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i; - cnt += 1; + GF->addData(k, n); } + parma = GF->fit(); - mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY; - b0 = mB(0, 0), b1 = mB(1, 0), b2 = mB(2, 0); - sigma = -1 / b2; - parma[1] = b1 * sigma / 2; - parma[0] = exp(b0 + parma[1] * parma[1] / sigma); - parma[2] = sqrt(sigma); - - cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; - - LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); - for (int i = 2; i < 300; i++) { - n = Left->GetBinContent(i); - if (n == 0) continue; - LM.addData(i, n); - } - parma = LM.solve(); - cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; - - // TCanvas *c1 = new TCanvas("c1", "Fitting with Gaussian function"); - // c1->Divide(1, 2); - - // c1->cd(1); - // c1->SetGrid(); - // Left->Draw(); - - // c1->cd(2); - // c1->SetGrid(); - // Right->Draw(); - - // TF1 *GausFunc = new TF1("GausFunc", "[0]*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", 0, 300); - // GausFunc->SetParLimits(0, 1000, 2400); - // GausFunc->SetParLimits(1, 100, 200); - // GausFunc->SetParLimits(2, 0, 15); - // Left->Fit(GausFunc, "R"); - // Right->Fit(GausFunc, "R"); + return parma[1]; } diff --git a/include/utils.h b/include/utils.h deleted file mode 100644 index 06f3942..0000000 --- a/include/utils.h +++ /dev/null @@ -1,20 +0,0 @@ -#pragma once - -#ifndef utils_h -#define utils_h - -#include - -double Gaussian(double x, double* p) { - return p[0] * exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); -} - -double* GaussianJacobian(double x, double *p) { - double* resJ = new double[3]; - resJ[0] = -Gaussian(x, p) / p[0]; - resJ[1] = -(x - p[1]) * Gaussian(x, p) / (p[2] * p[2]); - resJ[2] = -(x - p[1]) * (x - p[1]) * Gaussian(x, p) / (p[2] * p[2] * p[2]); - return resJ; -} - -#endif diff --git a/main.cpp b/main.cpp index 02eeceb..81210e2 100644 --- a/main.cpp +++ b/main.cpp @@ -1,6 +1,142 @@ #include "getADC.h" +#include +#include +#include + +#include + +using namespace std; +using std::string; +using std::to_string; + +void readData(const char *fin, TH1F Left[5][8], TH1F Right[5][8]) { + // read file + TFile *fRun = new TFile(fin); + TTree *t = (TTree *)fRun->Get("Tree1"); + + // data numbers + Long64_t ntot = t->GetEntriesFast(); + + // five X-4 block + UInt_t dataArray[5][16]; + + // calibration coefficient + Float_t CalC[5][8]; + Float_t Calk1[5][8]; + Float_t Calk2[5][8]; + + // read adc data + t->SetBranchAddress("adc0ch0", &dataArray[0][0]); + t->SetBranchAddress("adc0ch1", &dataArray[0][1]); + t->SetBranchAddress("adc0ch2", &dataArray[0][2]); + t->SetBranchAddress("adc0ch3", &dataArray[0][3]); + t->SetBranchAddress("adc0ch4", &dataArray[0][4]); + t->SetBranchAddress("adc0ch5", &dataArray[0][5]); + t->SetBranchAddress("adc0ch6", &dataArray[0][6]); + t->SetBranchAddress("adc0ch7", &dataArray[0][7]); + t->SetBranchAddress("adc0ch8", &dataArray[0][8]); + t->SetBranchAddress("adc0ch9", &dataArray[0][9]); + t->SetBranchAddress("adc0ch10", &dataArray[0][10]); + t->SetBranchAddress("adc0ch11", &dataArray[0][11]); + t->SetBranchAddress("adc0ch12", &dataArray[0][12]); + t->SetBranchAddress("adc0ch13", &dataArray[0][13]); + t->SetBranchAddress("adc0ch14", &dataArray[0][14]); + t->SetBranchAddress("adc0ch15", &dataArray[0][15]); + + t->SetBranchAddress("adc0ch16", &dataArray[1][0]); + t->SetBranchAddress("adc0ch17", &dataArray[1][1]); + t->SetBranchAddress("adc0ch18", &dataArray[1][2]); + t->SetBranchAddress("adc0ch19", &dataArray[1][3]); + t->SetBranchAddress("adc0ch20", &dataArray[1][4]); + t->SetBranchAddress("adc0ch21", &dataArray[1][5]); + t->SetBranchAddress("adc0ch22", &dataArray[1][6]); + t->SetBranchAddress("adc0ch23", &dataArray[1][7]); + t->SetBranchAddress("adc0ch24", &dataArray[1][8]); + t->SetBranchAddress("adc0ch25", &dataArray[1][9]); + t->SetBranchAddress("adc0ch26", &dataArray[1][10]); + t->SetBranchAddress("adc0ch27", &dataArray[1][11]); + t->SetBranchAddress("adc0ch28", &dataArray[1][12]); + t->SetBranchAddress("adc0ch29", &dataArray[1][13]); + t->SetBranchAddress("adc0ch30", &dataArray[1][14]); + t->SetBranchAddress("adc0ch31", &dataArray[1][15]); + + t->SetBranchAddress("adc1ch0", &dataArray[2][0]); + t->SetBranchAddress("adc1ch1", &dataArray[2][1]); + t->SetBranchAddress("adc1ch2", &dataArray[2][2]); + t->SetBranchAddress("adc1ch3", &dataArray[2][3]); + t->SetBranchAddress("adc1ch4", &dataArray[2][4]); + t->SetBranchAddress("adc1ch5", &dataArray[2][5]); + t->SetBranchAddress("adc1ch6", &dataArray[2][6]); + t->SetBranchAddress("adc1ch7", &dataArray[2][7]); + t->SetBranchAddress("adc1ch8", &dataArray[2][8]); + t->SetBranchAddress("adc1ch9", &dataArray[2][9]); + t->SetBranchAddress("adc1ch10", &dataArray[2][10]); + t->SetBranchAddress("adc1ch11", &dataArray[2][11]); + t->SetBranchAddress("adc1ch12", &dataArray[2][12]); + t->SetBranchAddress("adc1ch13", &dataArray[2][13]); + t->SetBranchAddress("adc1ch14", &dataArray[2][14]); + t->SetBranchAddress("adc1ch15", &dataArray[2][15]); + + t->SetBranchAddress("adc1ch16", &dataArray[3][0]); + t->SetBranchAddress("adc1ch17", &dataArray[3][1]); + t->SetBranchAddress("adc1ch18", &dataArray[3][2]); + t->SetBranchAddress("adc1ch19", &dataArray[3][3]); + t->SetBranchAddress("adc1ch20", &dataArray[3][4]); + t->SetBranchAddress("adc1ch21", &dataArray[3][5]); + t->SetBranchAddress("adc1ch22", &dataArray[3][6]); + t->SetBranchAddress("adc1ch23", &dataArray[3][7]); + t->SetBranchAddress("adc1ch24", &dataArray[3][8]); + t->SetBranchAddress("adc1ch25", &dataArray[3][9]); + t->SetBranchAddress("adc1ch26", &dataArray[3][10]); + t->SetBranchAddress("adc1ch27", &dataArray[3][11]); + t->SetBranchAddress("adc1ch28", &dataArray[3][12]); + t->SetBranchAddress("adc1ch29", &dataArray[3][13]); + t->SetBranchAddress("adc1ch30", &dataArray[3][14]); + t->SetBranchAddress("adc1ch31", &dataArray[3][15]); + + t->SetBranchAddress("adc2ch0", &dataArray[4][0]); + t->SetBranchAddress("adc2ch1", &dataArray[4][1]); + t->SetBranchAddress("adc2ch2", &dataArray[4][2]); + t->SetBranchAddress("adc2ch3", &dataArray[4][3]); + t->SetBranchAddress("adc2ch4", &dataArray[4][4]); + t->SetBranchAddress("adc2ch5", &dataArray[4][5]); + t->SetBranchAddress("adc2ch6", &dataArray[4][6]); + t->SetBranchAddress("adc2ch7", &dataArray[4][7]); + t->SetBranchAddress("adc2ch8", &dataArray[4][8]); + t->SetBranchAddress("adc2ch9", &dataArray[4][9]); + t->SetBranchAddress("adc2ch10", &dataArray[4][10]); + t->SetBranchAddress("adc2ch11", &dataArray[4][11]); + t->SetBranchAddress("adc2ch12", &dataArray[4][12]); + t->SetBranchAddress("adc2ch13", &dataArray[4][13]); + t->SetBranchAddress("adc2ch14", &dataArray[4][14]); + t->SetBranchAddress("adc2ch15", &dataArray[4][15]); + + for (int i = 0; i < ntot; i++) { + t->GetEntry(i); + for (int j = 0; j < 5; j++) + for (int k = 0; k < 8; k++) { + Left[j][k].Fill(dataArray[j][k]); + Right[j][k].Fill(dataArray[j][k + 8]); + } + } +} int main() { - getADC("2016Q3D/root/raw/201609Q3D1002.root"); + TH1F Left[5][8]; + TH1F Right[5][8]; + + string L, R; + for (int i = 0; i < 5; i++) + for (int j = 0; j < 8; j++) { + L = "Left" + to_string(i) + to_string(j); + R = "Right" + to_string(i) + to_string(j); + Left[i][j] = TH1F(L.c_str(), L.c_str(), 300, 0, 300); + Right[i][j] = TH1F(R.c_str(), R.c_str(), 300, 0, 300); + } + + readData("F:/NuclearAstroPhy/Q3D-Calibration/2016Q3D/root/raw/201609Q3D1002.root", Left, Right); + double x = getADC(&Left[0][0]); + cout << x << endl; + return 0; }