From 54c2b1c4eea9660ff43b60ec7aa7b3e1aba5ef88 Mon Sep 17 00:00:00 2001
From: Pavel Batyuk <pavel.batyuk@jinr.ru>
Date: Mon, 5 Sep 2016 14:02:49 +0300
Subject: [PATCH] PID TPC:: recalculated probability coefficients have been
added (from G. Eyyubova, SINP, MSI).
---
input/MPDTPCPidResponseVhelle.root | Bin 0 -> 43748 bytes
mpddata/MpdFillDstTask.cxx | 479 ++++++++++++++---------------
tpc/CMakeLists.txt | 1 +
tpc/MpdParticleIdentification.cxx | 355 ---------------------
tpc/MpdParticleIdentification.h | 9 +-
tpc/MpdTPCpid.cxx | 285 +++++++++++++++++
tpc/MpdTPCpid.h | 56 ++++
tpc/tpcLinkDef.h | 1 +
8 files changed, 587 insertions(+), 599 deletions(-)
create mode 100644 input/MPDTPCPidResponseVhelle.root
create mode 100644 tpc/MpdTPCpid.cxx
create mode 100644 tpc/MpdTPCpid.h
diff --git a/input/MPDTPCPidResponseVhelle.root b/input/MPDTPCPidResponseVhelle.root
new file mode 100644
index 0000000000000000000000000000000000000000..235ccee44f1332175ed7c795d493c12babf6cd68
GIT binary patch
literal 43748
zcmeFYcT|(zwl@k03ZkOYR1i>96cCW!gCZazpdd|36a=IbdP^cND!qw_)F{0py%UrU
zLZpTsB-DfuDIt)M<OcWNXYX^rGxm4x80Vim#>p5@vXZsdoWD8OEGz4Irk96@4+F!`
zN(KgomkbQWOAHLE-t@Ku{bHtHR!47)3=CqO3=C&07#O`IVB3m_b1NyljEoG>Haxxe
zzxDd`huSlf@@o$n)Q<=rabsX$d8+N?YR|9#`2N$!nvb1c8rge$dboSro4&Glb+wo9
zIs*6adl^`cXquGY@S~HX(l0Ce^@lbCL($6r4X|i09oMO_wm-Ne|KK`2n$P_7f%N?o
zFYR9j06m@l@%ii#Yl0|4TiB@pM%)}*!%Gve&~N8UuRL?OV{W|IJuQfqjR^gCL)2Dr
z<MZd&3rsb#gXOVvSDl+fh6dIijUk33`U?}!<q%@-#0@B3N?Bl!ODw#T&~P6Px~Hz+
zedVvT!n<Qq1WXpqjil<1w8WSow_VgpK*!eA!%a1CO0k)w+JXtpxy+yp?BRr)Iv{gc
zx)_!CK6f!$g%82T%O=!aBpP?$<9OYhrELpBU=sqt846|vJ3;sffy^LU#SR9NX_nby
zS72E-QLG>^+SgcB>#>Tb;KwapDsa4R@a?x{27BMi)Km!6TFZbnPSj2;S<MN~enx3U
zM4=T_9@dD=e6LBLeLa)b3>{(Ia)vx2wpOM{c-fV}d94dn4B}IEY&*Jb>Uv{+S5j#F
z<c<JeAc#O5@h_i|HGAK-Gf|Ye2&C*M5X>N(a_iCc0@Pv`03mWp!CJ~}HfmCTmRvFR
zY-ei1X8fHBxsF^=&53NGWevf+4Qx<Tk}l;i6=qX-ajl_qtwU3Afw7yfqa)7$Fb`=3
zTHms)vkB<}6~KH}iYnA-5wYPg?`|!6&+f3JJGd{;tirOM`rCI#7`ranKuYt)3bs;y
z;==j8dS9#xHf|+oA~WF-v>ngrRZ&<QVcT~V;D^rfcW$Qs#)W6XU)U88hZt7HiI_uL
z;b9YJu~Ys;aZl}@_1r~lF`=Z+pN49~{4|NrQU+n8<vkjAzq{mn70hi-q|R_y`0f{X
znEAg0n`hRKPGDzoDxxhlqpy4&<Bm;@*qQp@Z?eR=dUsq6KUmY5KqIQtgv)o7tD%oh
ztS_Be>>DkcC?ZCASD5*3Kg*;UX1bSQN~w)d#%GyLC5D+c;9y9rGZHs3+vMF$hB~{S
z@FPtMPcK9ERDFZssEWE;n#(HYhg-0T_$QNUL1py2q%&z0gP?s_*4_kNw<YUc)gJrR
ztKj`o1mxFs$u!E4Z%SUPVc@csKYaD`m~)vEZpLZe;(pl#be$hD*Tw7d9B=j8WzN@A
zdsPt?)W^o_HZ?LeX4Pj^JLatRs;SoOSz<=`$75+TnAsaDv=Kk5<@j{YPS7TBn$p}z
zB&F<48n%Hq5h^3<W&26AZ?nOOp{Ry}q`(7EPo_V48VOzsD;cq*5~I~UvFQr@<Q`bV
zM#%}v!sYc3K(!YG(<$K1Nkf$POJ%2>ydDUtn+7_Z^EdDnxjcsxm1wa|1MhdZlSet`
z6V<5z=vGy@`i9`Bj(W(qj?EUD+2{!xDG?B|ek|@GbRVUDP$Om-tVarqcR%bmtb-`E
zdXUP^RDWf|^<)&<eod4iGm;%f<dC$5PLwWAMSTYgR9;XjqD0DUHa^Gk8j%NKfz;X=
zc;a+LgUd>5dRywWC6To9u0U11RuX{G1OmLl&cMTAO^)CJfY0K)f=xmMP!FfqwnswR
z5Q_~C;%)Yl6y`IgEd@eGvGO!iXOCu_+X+J&aSXC04D^}iE7!xNpnlJ*p#lw@X?t6O
zYKP^`#$*6$YjOY)V($#6o=$1mi{yBSbAkGj%F_*WaWoRlX1fl7dOwXn+%6Ub<d347
z;lvz7!GLk_O1j|wB!@uev|zxhzCtUewTSYIAh<tO4A;fQx8lZE=MKJ%QPGxga?=U_
zm7P?YH(J3XP{$cDl;yKGC6d2K3+h9<t2^jBqgseKyp6xTGt>zS(Geh<6;YDoLnuRV
z+Q3*F;bR4KX|tkW&y2Qct$?bOuZ3?>5vuD56*O(a9_q7JPrF(RZ$<jIQ4KjB;Najv
z-U`wz9Rr;FO9TlU2JLTpP(om}=|*J1@`LRuNh*1aI-iapuP2_!+rz<$j!0@hO&^Dc
zgD{fthAXxB5UL6AS+I^X?e`k7Z3PKPoCc7#!fp1Gjg82#3X~T)jSk*0j`TAYMbdJn
z?RMLFv~9PKwP8W1ro(Zwf<45%_tG)i3TccQh=n<#cB{syT^3NUxw1?gZM#V`F2t$e
z5J-=YB8MA4$J$nTjiK8e-lIQ?bwU=yJm7Z*H?|t&;DN+UCkKLV1(oREa(i&&U<qp&
z+<9W;hm8o7WEXCg<*=ttdJ~3+Xrx)};LFtxz_3R$<c;PY1Rj%1JlMqH4>p=df3(3u
zda7a6Em?0ks?if1^f;}NW=unmBl+Ju2z*630U`VPnft4>?KBZLeHo6BWnKV<)RUKb
zLTx3hgoE}<1cDYrDKC(n9N=cI;7uLn8H?ZrO+4|tL2IobZP?C|Z}dkPLUB7Z9WW6I
zZC)W7P{L$)5^j6Y7VFG>77YcrC-mbV$dKWOxD!-A$3rr0dxH<qe0z}p_jYvr{$#WA
zmIw;IM_A2E4Y#R9m3Sytae9z@a^P_F!NYzFBz&u%_A?Tu_JXDlFqGj?Spo_i98~x&
zrc+7H_+Tdy>MqR2Ye%!6z+bQ(p<-%x$QLl%Y+ys{MS0Z5L5G5=YgKwc3Jz^U9LTAo
z&RAy2`=74MsWW3{Ju9T~nt_pt-nd{8+E9OSX!s&`=&$o%ZTwmR`y10bw?bc>{b~^Y
zLh`~fMkciQ0%{60-@H;c!Zd++$c}}6VCHES;IV{$;nBnS7vXr6B?vb(h2nA|1~Pqj
zIYEsSC&TCdtRH$jQf<~g1=V6r+J52+b)J7<{un$fSr3Mny>nBgWWlI@6&j-Qf<2a<
z+2-%tAF%FX_ivMLjp|=`B57-^cEGHC_w%}<#Lj?^QvFRs`){vD@5C{BOsuf(Y6-=m
z#!}`<NZptcsM;AFL$q?mT-nS(tyy3U+^qFx`~_T4SJwo79OlKOcP;`gWjob#Zzp&h
zYhvSWmewwH@<R`(!so{?8<fC#i@8sJ4{M*X<0@X=(0!iEoEEo}f7oSEZU%T_)7|Cn
zIEL*x)By`5NmD0*=j`V7*Ygka{Sr}l^SB%oCEgO+!c(u)Aton_H4#8Mwqmr1Pb5dq
zrElkij9iWXHH%(`MEX1nMxphUR3>yS#QKE+q-Y;wacuH>zTY<6DWlKp4U#=)fnPGb
z?(!pAr*$Z6lG9qxu1TF_0?p>DPD)?)jyR0a#|rlkOpX~Y5u=8I6+!)M;`g6d_840L
z%A1Ym%}c{0?z}6c?@prEByFQkYQ-EDnsFh5sa78)WiQWIyrLwRBV|#^^W{{sSb~>-
zrNq_CmXQ>z?n2&!#JCFw_JAh>2hp+ml#i4Tgnln4rf9mxZrS$7Vtv{W704&4LvY-4
z8%h=IYG#s}=B}lwYkz!DYk$ssLziwQ(TIq)ZOw*seo(L~;3$>GkjX9Ti$SGUI>Y7s
z&?%3xXh=UPF)k=EZt{cRGMa}hhk_7f+}qeQ4g0wV&}{~Qmg2KLzL_v+mw9HKFIVN7
z!k+P1#pn~s9@ZxU>u<Yla6!qEJ}Mafd<}lUIdnjlzj|tM2xeOm-$}YTz*jzFG(I@%
z(SV^wUYPNSvTrT+-9Zf-ddFt!(LrVE1s8mrF84<`pQ_eq6<9pHV-hQWIRkxHsXHig
zyAa4-k0F@0ATIWpT=cHi0`NvdK2EDoYl)Mi={jH5KNv2LX~hslCCNzwF&h2pMRuiO
z5g=;BZcg#>TO(MnpiSQIe$)@xTpOqf$XcLBv;W{+(6Uaew03iEyiaE@uYY`d!$x=@
z@j6H#()aA03C#@PV%v=&u|$9gE}rX7`9QW!r+K(ix6LUe*F`05<^BQlrzsgyp^5^g
z&eOxte!jwZd`EE%*H%P^x^z#^%08b`w{m1qVpjihrTSIVsj{)NSdd(4&P_?XO?R7!
zpq-OEsg+}NA>0!?jKLaA1cH$juxm@8hemQ1THjI<NVaRtT0G<=*r{p3L0M5d`dL6Q
z<s(vuA9yQ8N7;5L%O5scqB<SP{82c;3(bYSHqmv*<Swa<f7FMlr7T|n`-C1=a=$rJ
znPr>V>1$(B))wE_X9~79FE1!WB?w}#*+!kwViFTd@Pe%?l<?k)qa}Hf`bOt?*OxNX
zYajFkMP0y(4oQO4)B17f7H*_t<(WmuO_(EQ$7Z~=4T7ze$W~P-L1^k9Jds?yht&du
zs+cio=^>Fyn=JXcw+|fi`5j9F{Nf4>&DLd_z9btAt$!E=2}D8As+>r-%B;l}`Q(`{
zf@NtzMjTqvA649Q8C&Muwnj)yKMjsJ;Gk&qguJd49N8MbKdO5KD}x`FiLUli^{A60
zie#o22N)>Sj2K5<SWOjKK8&a|xY|uIXmno;F)oKK9va63wo3-bp*n#ta}Td>3twJN
z6m&{#&(PrKC3shk<~K{N&+;P>pdHs_sBT;$^hNQuomyt{03xe=8_~P<#=R;zZeU$x
zY4;uffG;N=e?MkNr&3@MQjn!HR?_5_$DgAr!lyYF?f2x)2;C>`H%x3MfUZG#WZ*`?
zjWWnS5)q3qu>PwCR2P-vqbX*cR=mqihrE8;WHtj`-P~uYENOCp__0;Wk8WFanp-qQ
z@B9@obQZg`U<DQJQAF)Vsg<K_=k@QHq?aHq4JNKkk6%13$@8q@OCVw+c@w_MaHI5I
zKtTY&9_?c$WSxcnSR_hmnJM0%)i-Hzt&0N9H@-4XZu`pjY<;TQgjBjN+I_IHn-NDy
zYU(N@9azpHJIShBhqosNoBl9_(FHWjEpgt7BR^dV*BsikT=R!hHY<?#Q<__s>+o`I
zcT;)@kD}0C-g1A{EZ1#D*Wf;issjot(uzk?{OJ;Y97Uv-ZA)g9FUokU`qw=K<nTSa
zH8Ei9b-<=1kUoGuuu?UQULB4fSf;N%>Gs5v439;1kqNz_c#0`Nu&}8F`NTSlo<RB!
z22KxqR6t0ztu`}H+N44!Jt~1ls;H*LI`<Rn7{22Dzffh1b;88~<38Ll1<S6%bgv)l
z4X{851K7-M_q>9F&1f#9-Il8IX3B_<i2x$UZmu)%#)R%gyH@@&w>EU3g^h^GW~(}a
zXBB4kxB}K2e_>fGhH9$ZO|@9_J4%%lgNLneNqMK<FE&6R$(44kk9=iD>N24A20&%Y
zWsj_(J0@G-r|248zNmB=ku$LRv%t3m-vEnTWvp~xY_zbF>?Rd`^bqR;4>ipqp2mT)
zB7Ia-;s|Hgi%D_5JnQ_yO8HhEl~+iy4Py|VW@`f8QfRLOoU+-ge_*c+ts{%NwU$Lm
z{#Y+jFj;R9XtphW2>e(a@Puyf-5KRbT>}8P&hiUnlYccSz=nr-XRGa1X^LIbAYI=&
zrRfdo(Y{9ccA1?7i_+)^W^?l8@gD|+-`I8bN@BGj7LwatTW1Ra^|elQDqno`JD=v&
z8<=c0)|+hg^A!hd(RC_g*|p;pVjAb@S8chzF$R8us?SK<LAMx{#|OZ6w;F9SbQU4k
zM=A{#L(U=XuwKnY{C1T4L2g6wet~LTAOWIxqJ-UeP>=V1zE8y=J;Tp9JtX$Wdx7kk
zZE0ipW-IxbxyM1nb)szx9VQ`EA8yI!FHyV*)2+YU*Hi&n<s+*<Ot&`NU*&llCC8DI
zjq^)(Z_WLDD6-<y*+7jVlXARTPb+nkA_!+!(;J|v{RCTZ=YhzC4n1OraANbEO`&nO
zwUm*Fz9!Y#tMFUetr7OE((p$qgp_ofb!v)A+Y=`skOs7fRkQ|M619T9rJBSLE#(bO
z)^`d@E?|FUWT>n3$asZ#)jt!<OtYIs^!QGKlmPSfF@$r7^+jdzP0K|pU_!6cToMGi
z6~EIP?@MPp9U_)oc{3g@p2fa|ArP!`b8K>X`9_)ecBhBD26z=?B1?G6*CU)Ux3c_O
zf4_6rj>1;x(bw0}S95hoQpt3ns9GSuWIMueXdUexLpY1HHpG~E<SE8?63j}Yk8;@s
z?A0O8L_jh*8m&5V4jsa|T<N}jaGn2ZcT=I=IHd)i1L4X&%>F^wbOVbmlz7QrIlnlu
z&V*gaT5f#sK}vzB`^G!dAcr_d1G58=fYPm1dXDI8R3+U;@PRV94*0y96``Fwiw(fY
z3q9bPcr;YBKHm2O@~-IUu=;8bD6@ElaF9KOpwQD#xw_>2z|O*i#jf*qQ*9gYM+$BZ
zUkGW_hCz`1yax3s;`Mr0%t)@N&!(nBiBCveRSPYG9H%&Mceyk0`VcP&%-sXEGVCHn
z>{w;!Tm{7K6m3fR9fmy06{cs+9i7fzL}ix$_*@(UvAmsO?ID6mxuex}RlD_S^PR>?
zQtahh=1)Q=4GF6W^=?>qNo#q_E^m5{$&$CuL#v)g8pb1*L68WvN_i4J`5g#N=;opi
z&#wD6!WAX&$B+*Td2!jQBbDmK?ijS%P<2sn{DIDoE&}Dzl1|KaggqozSbjYEFga>@
z&&*jRY1%`%BrR}yubG~N7N<nBGI!8z*WzY`y&!b|>e=G&0cD6>y}S$xNm=iWH+Rcd
z-}J6(b1#O4DU>C>Lec$30BEmV$9GUky$Z|_4<!YOm3Iu4tiTqV#MF`8*sJCBiy`N+
zSL1dRA;Q+(=74iKkk0M>tbTOgLD9nr<P#9Tq^vf1R^RR>*`@dbA-m-f-J7<h_ex+@
z!1X+~PHG`kU?XsJE)E0BiCHxV37%gKGgnyN$_h`sYG5LX!HW!Fzs$QQmha?%;}Qk$
zXa!b6t^zR2&N~w(#;t|GZu(x-VqR=+to#<-4BGD>udz|F`(tX8VdNt}dPfvwlZxS~
z$KbyezAqzzhv+6kirn4~kQfrYg;_q#DAkOSl<xOgDmct8YWX2qQP96FP2bA8LTWgB
z^gdG?zYaJwbL#aNnt`kpYK%GqP?8f`onDuh2%J4;7mdf=z=H3Sd~ae?)(bBj(BsAj
z9amA0e#x2S^{UL>0QOFJ^?LTy=Kax?BGGR0rjAci*?!xSqK#3qHB>#f4Ig)?Siyz8
z!7I6cU>qY^$A=rZuzYS{yDNVH((hA2<*_`H=!Ec37kHj(OwB>vW=Tz8Qxsp*iynem
z3+;(LvM1PJs2LD($R!zd0UdRr%br3HScHE1zmjpW8LglqTz0GCMJ7;2ZghS#rvbu%
z9xAx)GM9Gl*MHqHo0C{ID_dDb61+O9Y=+7;KOir%hn$<ScuvR~&xS{+@fAv)W4eWq
zwpN`3Cl3`#WALpvCo}`WABr(56NUyl<{Mi)1YvDJv)GvdML*=Bo4;eIB6ak-2}D_@
zXPfe&#eH#Kw=w;J&X1{`2}28;iRj&XcE(V4w~BdBwwQiu3@AU^1IsjlV&_>+4nYYt
z`BJU|`i&vy2v@ZOQ{xWbc2|Lon||QOp%{}FZ~J*whPW@ibLYL8;MH1l@NTn&3(>WB
z;3?rTU7>%_XmOdoy>MS(3BAHqi1LqLwe5Kt=1K~j{jmDe=P-qym{`-cb6S9)$2055
z0D{ac6M-_&Wjq%;ff9Ris(HPEz7f}|?q0IRD){BTZf9UI3bFA*yZV=GwY}`OTl2L_
zS+-a9ueN>hbH=hVaxLH*GiGPN+o2g93q!|$>)A5zrf-eMEMGxHMdfQHi*?UW-ez*r
ztex_?8@EM$DWG40t~U#;(A6SGv|=RdLh?w~8;WCs9HAu0_41zl;=nEtPbJcIAjrPS
zjJ*7jBPe*_a>ZWGt;aCnA~|u}dntda6Lhh%C@X{see*FCHY{W09|Qsz57mm|LV*4q
z<PY|gxA{80bFr8P-X`E!vq$KM{AQN*izb1iK@n(S+bK(9h=pNV`%U&S{3@(@M|^QV
zzQMcbsU)L&>l^Lr7AiQv5VmRXdHrxiR;5QPDdrHwRfv!@%$#4|wYAQ*5*ulD&GqbK
z3Q2iXpM4rj*j}Q#c0)Gf(G5iD;vQYFn%}mMKQF;WU;yXXi0h#Q$iU}>6BfzV@#*tU
z3~|Uz#cQlSbfqFDw=`Cm4^@6CP9|$8gPJV)0cI8x`u5tek+=c%yBLShOzXSO-djpl
zq3*NYtTCT7l<k{T7WdUZU4R_mp{a2}_ExW2t1@LnpPYSg08X*N+H<Yg#xR>c6WX*i
zqN-_&wYKI{GgH!DE^#Dj^y7`Y5m!{)!=3OK2EMZ8juxNm!fUx*0yf6=`xvVnFn8+;
zWttD(XfrEva@VsK)N@11o;hIlE6`A`oer{!#|>4UvYna03w`ME=Ub0i-&|TG*2mx<
z4aT%@vy7~d&KWF{#dCCMXTsk+cplz<4^Q;3`5o21bh0xj)rwI@<f~`sev7>^JJMcX
zdRyd)P^Ro$W}WYy$Y4|kuULXX;WKl|=wxH)d6vbXA3>>|tsX<s+!6TZwOSkce^I|?
z4Ce){(0#~2KH=vFK$N$f1Jl4!x4AoDn)4HL=<XZcPkBn`*Wxw$;iy_b>-yVvddQ(r
z0}vw{th+iwnNoAGxZQitnV5~zUhAe^4OM6B+*wlhE#PCayS)%h>JB2NS;1HqiDr9^
zmUHJW=rQ8qG0X09Me6LGt5%D|NMd8~pP#Yz;{GMNWeFO}&To|`Zin9Cia_-!yfc=$
z&IINx!~-*cB!}XhE$%|(0B4VXoP+jf6qA!RJa|knbc@pyauE)8))o^<@ItC{RC{#?
z2>?RpD8XIbhGQ|^7DAWr3{}f0Fj2TY)h2+3t+Dm5cLoCRz-ny8>(xFAVQeOVJ(qx1
z15^wDq`M0xrGqq_uFe`$uQ6Sa23zezx55U?;Y^!EGPbCGi>t7rM?Plt@JEyefCYoC
z3?omKw6Vq%X!MVI(WJA;$WoG6!h|U$?-*&OVxBdoTEk0CeyABaOS<+Lg%=zvIOD7N
zJcb2B*ygiS)ioCWDH0mg+5lZoDs-nqUG8gyUR-T;b*G%C`$jy*ES_{ls*Q=TvvhqP
z-WkLPydop~Rr50h3OARW_;yCPg7u9?Y4lQ{S3(5yB|4lt*u3EfjG;#Dt1lj&b7I<I
z!H`pqt)}!{(Rl7iPpq!at0A|Tm{-`xuw7*G4ruA@P3Dn~l9&bwv8d~8XB$#9N;!#S
z;=Qf+8ey+@^p3^2G7+zwT}^I?Z;ThMVm;f-bPJ#phn1FMW|d_>+Nj73$|au_I*Yxw
zX%}%;=%@Hu<O@)iQ+d%FhEUSsUw<)#W-$?hwhk(@OhwN=NEIp{_}eJrK~ZDWc_!v6
z4X~OveQXxfkfd7VQ=!M~od`Vgq79wY&9kuZO7}~VTxSbeS`00W1MYq<cug2ddceX|
zt5M2}w3CbErHefB)Gm>UiTPxxxnxmq$~BGWZ5mebgOZ7wuQhHhP$)+C(=@8-l2pf@
zXUfya0QQGC9f{wGuY9l*y>wRSG3zyWc*WVoYd^&osCiS`Y;@HjOg>&6VNPQ)@=l~f
zJQ*4Qh4+5Erm>*Tyh!a@y^TMn!xMV1{UDDa@!D7R(ES7OYKpbct*Z$vixe=n>&o5F
z^>njpjU5v@%XG^of|+%fD*{`Y#cd*_@p(Du$6bx*znGkUnT%hP+F~m2xj{Y4@RgMj
zX&2pY<h#hSNPaFNz{C*B!m0=tXVB$lvJg_|I7ZjdJkx}php^nyY3K7Ty0)SfMs2-V
zMwSfbez1u@6N0pP4WE`1k_u&-7+~J?F}A6CO{jXQd|KqOkp4o3ermfBi>jbyafV;g
z&4r!5qz2?*U=fo#pzf&>?#7e$1N-yDCM8Le_MOVSekX-UrZ#`qqhC1b%$8X)fTr`4
zJl#6&#MpR5n`}}t?*lP7o28JHJX4lP{9w6kA|q=9^ddu*eoX9ZBC*H8ERBBq)GGco
z!_<v#Z8~vRzr4s(8tp-UA3~6pZyvdSTxj87k|Dv()|I(5rE)6yhBQN{>}$74=1|Z%
zPnxtu@Q5m1SSc%}P2f8h=!HQ~_5difffzr0I^{J%+KSxq>YZgfh3bhH%;A$@BuaMF
z0%9B*SV~iV%>%R*tmu01#Bt<YbZ5lc|9uv|p1+*(Lbl^}!n@!AZJp<Z%p;FuyPvb3
zs@G6{)zIVU!~9N{PEiVX+U}xXsO-iuS#>7OYYRJ4nl&<~-qMel?XE3Q?PL7%i%J{l
z+Fk6Fm(r%Qxqy3QYyqOP>56naI*8`HwXpM%vvEd0D3Zx35+^GxMgKf6;x@$m@#aF1
zYouqIe(q~T$D<{8h4t$v1F~A~D-Iql>(N~Bc~9Di$b-2@Z|E|fYo5ABzV^Izm(56l
z^UCwqN1g8Sed=tv>VzD=qvLKNPY*`C*}P?k%h9*B6b+1-4n0~hGC47=yozEbU8h?#
z@I5J&nGryfW(%I;XQS&gfDw=9n4jg*F%yd2{280k-p(?M=g62ne$;u6$?2+t(|cx;
z@Is&$l%H_#2$KqXZf24ZGwa4Zx_MYH+`LLh8hV|bZl9?JntyQo?fm_?c*2<rzmL1c
zwX?D}pag$qI<6g^$=im?cG}6vGIiXb3qQFRvGsRRla9WR>@^%n=m^~|sH|hY@kmI2
zGDJV6-HUEJw{;9>bU7a+$-mgSe#9@7H8<3?OZ$)ck1)^T4W|VZP|K|5k^{q|ej{`X
z(w#xC=>6?Jz^U}t*vc0(!rd}In4CVHaCAM=l~fho*|hzT7x&%I(r*nMYtlbYH*O~8
z5#}GmTy)trp8uv>+Hnu}=p3=<tudQzGTYzS%_YMW>br)yALOX`2#fqpOxBpU8keFi
z=*jULD_zW(>av({jcG^Kl@~lmeEu;Au5pQWEA2{X8(VH0Ay=3QElD@`eY~B_2JAHv
zNOCMH*y8<a(+w3&e?g<er&?f~4(Z0bh3av5K$my!%EDol%9xw?Q~Hlj=DfT0|0dgi
zF!<95{YJMRRdWwZVU5@AxvaUmsfB@bvHlEw_mpnS37Z!>E<FsD{g5Mlld0p@!l85q
zxvTJ<G2Iu3D7H_J*2kA|HRTB^tQ+uvZo&!Eb#fTGO&8+tq6Ks5Hu%78=4ieD++LyU
zNQthak2wXm|Ni~&9TyB3iR6w(V2rckk^jv~YRhFD`QIjQ&c96l4=S(Qzp~L`y9Rd#
z{8M3PUmcnLHDTBhuaClX>G{>%j@IR9+5oSoT3v{y+M_r+Bqj3nqknWTv&J;NCjY^2
z_)q$Ohx?VE5z+Gxw6bi9aM`}ze=OoD?RLWC^9;8mx29X-p|n%*YzUpI)CIaazA4U}
z{Ld2pV_E->XNRQ)sD!$I=3fTUl}oi(Bt}0u_x~FEFGlg77w^yH3tEgYLxHYVhab&)
z?Lkx4+jS~Ybk*J5l(W0>PrV<B*39Dxhr0;o{abpa8Lx?AW<#cke>;kf<1>(H9)(ZI
zh5XA3{1dp;{83~lQsq?s>EGz1;FqzN-}qlh@z0t68S}qG@ebHdv_3^sm%jH;9{+Cq
ze+l%z&+gCkb|IES_2|3KGd%yYa{pwiZPvcCG(6|8ZSa5C(*H+Bf49JY&U#+S)CP5}
zw8hX<>i-n+=Tb~<f(}wcb$`P2KgIZ;(A57S!2cBJ|E2Rk<Ujwr%KxPg|F2s9S1tcz
z5%&L|!7B6#=RPgg3Mv2R(OhN*1_kh+=U6rV@xU&_OO~gPoV<Nb8T)wIzi_kn(sg(6
zIDg;6*VZ*q$koZ+#r`Gz!`p}cV5`Gf2CxDH{lQjCYgc}xuR6b1<6-p6hhy=}xosvH
z5!pARb7`VeRpDG+&+5eUeu?In_D*>EzvFlDf0vfZt!n|22Mr$wh2FS(vOXe=f$e9C
zdd&ktOr6GX@vEq8mbCh-hK$opx1*Rg^@dZJHvNZYQ!_0j`m9Yh1(iDwH#5nZ!9$I!
zS+Fm<2v?)fAlCG6CtP!G5JG;Rs~fT9nlv3eBM}g17(fDFl8x;CX=St4RSEZDyDhpq
ziV9v=YMBW$^<SPF81NW69$MRnZ`Exjo!%6*loWiF9dyQF=%gEjbeyml@RsCx=z3E7
ze%JY++JWJl&zkHB2Lma3{V5{*pfkbl<ZfYTh&1`Lhr8dFs=#aQD6jSl1=upFUEaR-
zEp{e~IZxBVZ$%wD*=18#Nw3SuGZC+N$S*c#+_fC8F+9rVJB<nCS-mdpdDya0C~g2Z
z_R1+%%MQKT_C-QUDk5w3(tV^S#HiC(2q&=kQ1L;Ul6ZyT=YwEH-e-9)89t?MC-Yao
zA?6LrolDj{erYxVl<_b$;?>W_#m0Kf-rIwYpqsvpi0VzwS5x8OnKD{Gi)4_nbfoE}
zm!*q`$1f_{aCdHr+kQ*+&iUP5#dnJ}|JjolFSQ7d%gA@%M(o+*nWHW?OhsG_(7O;I
zrL}(j(#=8bVJ&3!`b^gj@6UqJJ>H*LN{Tc6cUy$`-%P&UFP)6fpUWOh=m2N!lt84O
zb|s7MT!%cr7{Tx9-pu{o1n~U0Z$#w<H&fgN;UzsEQu`enU$J55n98A)Ogkx}pj+oG
zpyo7@ch)&soz&M^-(Z=v*7430M*WgZ{;nYVmG2I*A}O@lOv^WJxe)<ANhB&a81Zb(
zgCjm`3Z1)qyiiX`Dt_LrQ#D8@z(?`&DK7U2xm|9t1~#<Mp)sB+!<G_MzIb|K82D=0
zC_Bf)S>krmi>fOxu_~#4kM0W_%$EFi0)||c(QUNpn%77gk}Jv{t~##Zh~<=hF%hAS
zeQ$i>TUB>p?5Wx0mrs>C)Jt`0*XAq3IWJgozObI{ddZbh?4gx#f>L;8`!nEt;D+JU
z@fDtjymL}li3sX0N28%E^|b7}n``N1*{?%yY`W##_2F2R!*mUw%ZBNauVt&{X_0xB
z9(UQUk>-UPN@i@P>OPj(D(`yBsra8Mt9PZ@fcJE!!u$GfE<ooZ!naqYKN#kJz3`A-
z{5s_mhwY<|yp`9G{J^5pWvP>;Ub!EC2^4ybpX?CF^$|*-#1y;;GQ^;JjQ!nSNc}m(
z6kyUFwT#QgJRf~T7J@kP!+xYJn_juTW6yGWGnICIK58bS^?39t1{+SWV%C|zKAS9N
z^6xHW-$Xqbb~a`IO%~rXu>N2f)ZOs%(|vo!;Hw<n-dcf#>8i_88P>Hfzr*i+Vyrvy
z?0})KbQ&~xyB9Mg3avcnGS~8;aJ~9-bSL;KCQ?8<{+;e_k+$mFGc^yxU&{0~#PWFl
zb;?kgU;o7GGr@Y2LqD5EfCs@fIL)M#VXmEFMdb8FK873RQOAj%<!2E+$T!zTMg(+&
z=+6v+z2?h`Odk0`uYYxJUc<6iKz|==Np+8G;~Cs*Gl&LRiCdb?e5?5#cCOKM>oke)
zXYjNywyV7b82;sKJ;p)UhHro_2-5x)=hPSOOe6D(G7o{92Z<&o(S%FUA#;9ara6%%
zYfWN`h8^~l?j(x_#?uPkg}3K>!W$~5ctAyZ_X`JlA8S7K=zLWYG@7f$q02$z(*F|F
z0Ua$4DqSlL0&0f^D^^5f#_y{m953b{`x<$EMO8ePPqc^QL*69MQy-t)?gMutU0F~@
z#Lx?bn{tRA?E3Vioh`+w^ID6<wNAY?$Es_&>RUE`w#nRArY*{K{47<<=2mV$IG*r1
zMtZVJs7oDUaFHbvTzf~o(|vl#bK20AoqmLv#U^-n;p3ZnrRp?Lo_4LvSDCcE`swV}
zG83Bo#da%B*;1MPj)~~jH7n7Xap!<+CyNwoNq$%sK)q;*S5939YL(gGh1qsfd5y0e
z;$JaSk*d1!<Wn?Sq9;@)3n+Z7M{$0Ltqx^TB|q!wqr0VRvhkI_!9_~xz7G#MzS5(a
z`(brJU+4Le&*Ax=^NcIS?N*eA7dD)qkW^zS;D~K$Qs|w7tKHQWcg`qzn|e-x?x~P|
z?C*|^jlGdn=&`qzunZNgc*ba|_h_`Zb+P)?=~na@L`=f#g9s1#w!X%Lce-4bZyc+c
zPf0VL|4WSrt#@kg^MaY&;pxzL<aaA6xrGekbA}DiU87OW2-xG45|%>&w)AR<p~Q+#
z-5cJ5N_iKS>e`aLqVjLf(qFmHZ6sN=VUj&phWVcbl*GL{H}Z(5t2|W2`Y(lPiEnyy
zan+;O19B+mTV#jjY*pgDL<goAJD>Bpm)WRfBQSL)Wi!$*CW;|kAir`wkRKB!=+3JU
z2*c;*ZNJIJr4zL*TQlAJBgIYr;&Hz=N1^0reCg^~&;)n>xAr@xO9k5->Z`$tKn{WA
zQ)N*hCZLYm;#xMMHRSb+CO0#SBv?<I-jYz}k7Fk#h~LyR)TBMFk}vN1@HC)8$qnk>
zvhTWLzAY_Bj{2SP4LilFHBwx`?$z>sm47AMJ219Ub=n634s5sExS<>bKMe=0Hru>$
zWj{{2e&xCKyS9YOq=f3n#B1gpf`Ee-NSG2(tEFAl^F3_5Axt&DswJv8+^>CentgKW
zqgyjUFYd+n4_&5@R_da0UuDev`BfG!xZeksNErm|9{Vg~wy2^c|KMQ&JMZGpJ1N`}
z&h5fT>8Dm^tO8;Bwc<XZZ*Fa<obi`pdw%0*)u+qn?ln6lt+`BpM89cPxf+*&Y&p?q
z>v4{4apx|>pyDMl&6az@QOjd@mKLw}8u^{wzgt*%%+Q5v&2cD02e~$}8ux|sEne;#
z|IHSWUmAMu5#Vd2J`Nc85+$d>i1p_M#=9@{`d%`hHt(P!nVD6qT0<tU3+-pq&KVMX
z^zL!gT~nyxuOCb&s-#Tn%#ZrW?el<rH@Wc`dOXid!(M*fT(;34b<Lb#p}mxaE1}6t
z`m9&)?vEQ5-d9e9S@4~}y14X~UGwY>oa4O|B6p`(p;s!YMJt+C(|LF)J#0wk!pkl>
z6@#VUqNqH!kBRpBp6nA$-frT{)(&%wrmIXPuLRuFl_eB;zkQ{!Avbv)&ETa#Q3IJ1
zntTHeESR&kV;KQMmnQ4yj+e_)NEO%@;N(+DYvq%rDb;$ttPjKFg4Jz`3cz@-#i*CR
z-i5tc6PN;kG88#^wd>yQqkgGoUT(0@;v#?5Yl!roDT59XOU~QN?nX7;rTY!5T;nTW
z^FoW~yf>8PJ+$^5Sp@q|sD667*L_jWmo&EfqiD2XeD_DHNS3NhmtSjo?S#|ME9I{C
zCYVu+IIb&e_WQSj6+4FP%KC3V|2;9J*ZX)~OU!DSc;jb-XYRR~HLz_?z3E}+BNsiE
zjxnO)D-M}CV*Sq<+1Dp-oZk6%ja2+`<+8LVFB@&8x~7A(w4r|Ptbr~)0JV)co^R6E
z4=^m`Jo#kljY1f^N7`KIP{BhJ<0t1tZFofm5Jq>Fe)lH05yQeB%w4_4!m9)p*Lo1m
zTVVe7bVmK$$fewzt*Y|_7ejxU+cGsMusrJOzL|Mt>XK7<TGm=O6aRGUnKM;&s1uJ>
zqB93{V>>mYgsI_+=756e=8uU{>{i<Eose;%N;*7mZ6`o`^Lt;Bey9g>gBt0llA<rP
z%&tgRs(u?#ker@>Sfgy3oOJVqK=Gyi=+cXGr^>8dKa6VCdF)SLaab}>N-TW~!c^&X
zT5BnA2}{cmum1@7{Hl?3ONI;0IfsH_AZIGPDITk}OApTnxt7w%#uT`q)L*tCKhN$w
zfD5l|J$AnPJo31)=0hl9?=2uSV*f3$oBeD2?y@QQX7_o3#_+np(9Dfwbk(~0)~>=7
z|J;4$^swvB)$5DSB)vsV9t}N20pGn4gz&w?whsXbBjFE!s#3b*Q!oLw$ZI5p8v7=I
zYu=2;F@X!>#Xnn(Z@pUSePos=WF2M8Jarxu_;9no!1h)gT<~ESChDnwK#=>CdnWaF
z-$n;>5#2&VimrrQ%HUL&6LVNF_HcTU2H9K&>|e3bp0fS@M3yD}I7e+&q}nM5&i5a_
z=Ra=1MENQuc?q~3f4HGMRx)8{)$qK65t;fv)o5tHdue4UvrkUX_IXPZs*N_I28V&)
z$Jn=Tymp}hs~(uw8~0^r_6?Pt1HJND%7(|ivY)C5PHs5t6uAF9P)z=DQ&hb2y&CyA
ztHk}ayF%3WHYZRwIyzH$vx_7h<Z2O(HFaA)qnegzSHs@M)Tn(fRnxIxr0I3bvy(2p
z$A!OK(N#MI_jYtMdwB8;)Q>o}KXr>y`~CM4bH>64MmkB(&l=;-pS*JRWy*V4lQ*B{
zU{d;bPRoGz147wPYj0Q-@44RH@3p!cU4G$f`3rDpY;5;B&yydL{4Tbp8_Df)-uN>T
zqL5f&BstXG>}j8gZ^A3qvYEZD8v&*qpGHA<-{&|wp1XZ#8?up^HyWb0&3&=)`YFVL
zw&sx1ONqQOPbuAN($FV*Bk`f*<Bx*jcf(QAk!`mVJlK{Hj;U%(Awwsx#dBvh9S^kp
z=9^r8NBLHM7|Vgl7nb5558iB#j9Q=ZU=aMeK=66OVBMR>2FA02J=npL_UmHO;EINM
zW6KE5{NvL=UzZFW3SMM*dh(%d2I^p6_mF$A3&<tWJ>G9s`BkW2`HVnSkmEIez=P8f
z2I4-qr?RsYH|}Z_*Vkk_0whQP^OJwE%YMGJXmFaJ$G<q$=(S2RM#)L!xVqjJ@TY*f
z$(EocxVZJXGWjX(m^2(0g_3+Ja(xBVRi?kFd2<FmcH9vD^0X0p)Nz@wNLS9=RCP&L
zMby7udg>Xp?ljp>Q~KkbQHzqXMy+{n)*zg)PBT2Rvz|vi?)eT<DIrIv#;yMifc<XV
z&3>W2Ru;)zv1@N6-hY{szsP*8XcCaofVQz}M^NfDSE)a&V|f!_cs=>b(PMS)c-6fo
z>6F98ey$7=^Yr%Rc#Gxqd%<B~QBx62H!Mg84c9uJfc*=qUQ&AOD`-3ZVnN!y(grrR
zTS?2or=Pl3=64RJ;xem$SvDF?X`z7VzL1nt{3&k>ts+YoGI~4#Y+*l0KGxTZ25)dO
zpmr*_B#DZ(sN=h|j7JNrjAbQ7Q+FaEW|F7qTsGTo`ONl87rc%ty4+IrHty+wpcy*4
zu1pkik<9^zj{-&^EEf*KADlPYx08D7k=)6p+F?v#-n$IWh8$e;QK>?WzNCp>ba-z%
zn3qS+^j5XFakrjTsy@3fNp|uP%c=KGK4$NgQkU17(`x+R)F`~oX~fhkiksd(wD|dv
zn6G|@#&-Ru2!vDh`?dDU!Bz3cCmUS8OyAST^z<pbnRA&j1)iZ93Ikuyoqy(j(!~o^
z?O2Q^rN7x|3KezbzNjI1$v27OOx-<YHx1!iqH%Hpxd$QEVbHKsvta9Xc*3~PHP5Z7
zZ0`2VnD48LE`AQbRIh$||LR_WNJ98hi=z6kOX1u1d5Xm`lXG^5*O9}VB5_p^|GRvY
z6Zzx%M>WF^^s-w9hMRxNcK?(E|0z1AA6NX}%Yly$H&mw>7}Ug$S}tt#tMI55Z_mJx
zAam3zUh<dwz(AwX%Kl!!OK0`q55tpx&PV>naOtSd_rJ?fP5g)?1H<h(I-S4O3(qtB
zsTWogKce%m^}=3%uNU^BAFJfu{*QX$nrkNDo-F}^j#ClR+E1T;<h}lqr!LW#ITaAi
zdJ@DKsVsi!ZoB4vZPl0uo^M!BoHmFNyV;ZDFn=n#D)QzN`ES>63O$eE`;?G#-0<-t
z(^f++YL_w>3N>sk#+2I>Z~ChSkmsAQ1sFhmaWjf<=5Utp8Rx|G?&d@{js`E4>v>jI
zW?4N=bD^cas$CgUqP{_!@e@BTKb+}Vt&^W#s~z>7u267t8!IQGaAxVXz}kV2<*v2Q
zs!ee-5=xD=PGlNvuVMaKwOI?KtIxx&HoD%!ZMK})<9mX4b7vjzzW4ohoXuCfhr6}h
z+_zoqUgFy_Q|$g86ic1_7&M+a{q96?aNoBZko<BS^vvi$xvbq+8ND+THo*Pr-ISo(
zFFuQoejVSiXC_aR8+Yp}@L*&&a^^NZc(l>qn)X*Z3cA~be%KsTJ&}?AtXkz*Jp$JR
zz$GcJjsR+D{DHq0_X1PB`*yEOf!tlnYo8fohYXsh?~`A;f!f*VMZw=03o=OWlJ)=_
zT9SaJzwDa@ss|Fo;ndwEo8Nz>z=FEj238TkmVrZVQRT&&<pO&9Ro8cJb{_Fb(3E?R
zxVa*9YSGS5bo5>_hYF$dU?2<nBU_}#x@92hataResF}FOTyWd5o;Gu(#a$%R{+Y{!
zO_{#Pw0CB0z(kZUt?W~+=R|g2m;rbOmErXrRABhd06D0@Q*JKhl48;PsHULRPCWSg
zTBT1H>7<8;LfWee&`n{S>@bgof>H2~IXwgKVLuS{g~@1VWoauZbNa!k+1xiV%Q8)Y
z6a#qHxP(hBP35J%i>V5!fOh;<E!Kd)%4_k<@1ThvxRnrcBreORW0%d9SeZEe79y|y
za+wE)i6PXm=VKqmtLZVyU{_|STR+6=o45Pwwz?3r2VF?`&pJ{^w`3NWiqor@C2REz
zs%&c3_3h!+m%forElT7HYMq(`OU5CAFUUbQydhDc`rktxE26_XW{?L(J%alrV*@7+
zIq=J-wAyEe7E@_%$8H5S!G7)*j%wD?^4(}2rY<>q*mbxQtgd_)1Ra&2ENYEHCF&Y)
zQ`n#{h7PQQ2jNb<Sq(Hj`0w|$xsV-ihpC>AlC)kO{HN;OP53O4Gn%%nqX%u>SO`z4
zX!2U;hQhlVPt4Z7Tpve@t5H4yGF#n;G4A;_0b1m~GiG?u0@3g536ANNY-~kP5@nwk
zK|XBVX>$XR=C;v>4gPv$>axNR9J(j6(M)vr-YsxHFeDFEA+2^gU@fDYb9Q?s#wEmh
z55a?x4j83-bfX+?vy)8m3rYp=Y;e#EgnQg=Q2`sFS|TktCh~MyiAU>;sYE|=J-l|n
zP0ao4OkH447rnlhW2P)c^{||8D6rREX``0hhUoxSro)0)Robf8T$@37usdomOBkmo
zQvR%QCc2`UZ*zZ|(HPY{-xjbTPj-bML~!hu;LU*o%EG|j0jdqX=R#B&C?}8e(CT)K
z;V+SH^0ia9W+46Q^-J|CKbm(T&3hXMJrd-FQCu3?oz|O!P?{;9S{(^C2%*K}FeL)y
z=F4Yn+|2NMj*hD6Vai5vx}++(JGic{#HJN}3^=oQt$8~U0pEnHQ@NVA+cgCtOJwqt
zBLF;79t?JU3-+7W*$mn4;GJ@m2*6*j&+tqDuk=}>;A;ggxTt6j$ncYN8%MK^raAp+
zI4IqkT84+(6UlMa11YGQ-`PaStOUxBy6&NNz}mb~?rnIu+j-G%DlmKyGB1L!9yk*)
zqOa!%{V6fCRsn_VtV@2?g9j1&D%uVhn|;RXJS--j*{8Janem3W#v5;py!AFik$(=>
z?Ga~M7Uh~@_AYox4WVB~va;3%Nc?y#Mb*ItVeysV7wukC+v5)1VVQxpPi!{-{$=fN
zW$RK8>7L%n+Bs*r6JFsC>cz&@$<ekp3bJTWnzM^L4foezIrVK><M+C$m((P;!4m=v
z=t6~%@revybdCpvs$I7di<;aUhi<7`?nJ@hu0-K{*FlKoP$TIxaWYM-+3ny%wS2uC
z;#t+<XOB%gmlbuF3RDF(PhTyFtm5E7bxkKt?GWddzF50rf^`Hci7Fg5MJm+8Nx_Wf
zdVO%<c&*NyPmR7c;a7*d#c%F((&}Z8gT)cx%$N?a2^r}K?#-wE`jRtGB<Z+eG9oFG
zhZq}q-|k}KCX^N_PpIeJ2wC6Yof_j4P^P{0$Z|2Rplp1JKlJI(w*eui>Zzl@<bBt^
zIL>eN_Q}`MWZ<+89LmX9Z%>t67H;**I2sVR`Pb^B$+7QdQ|W?P9>xfA$uTn&v_Rht
z{OtvvG#VyPh>6s!0DBsyPcKGSG_AHw$=9w=Bg$sI8};>P1Ch8@+D}s&=$k-_>emrd
zfZtqrIRsHzM=~pzt6AFSAO%dAmivHZG5B^0B6#^S;^14$1Y<?>?}h<H6Y(p>%wS{p
zy&DQVzpFM$nr^wrz{q@-e$(T>cnf{Owx)LYBDP17Ro3Elfb;g$Wm+wW5=&FEeMjIH
zYQqF?ehLf>Ev_IIlf~dyp9x8+Gj(?CXRtdM0Yj@mPwOCH`%i%<nLbwbW+MZKgD!SF
zj%k0uqU}qGbDmAKl?jk$YHCZ~Yx)hDAu3DY&aUaa#?}d{P0L#giK$oh@2G~Tzz{A7
z&=Isy2ey~EV;0L=W+)GGHo9>BYtLvyF6-{qcEWV@B{z)3V}dWx1ZqMsChl|`rr&P&
zo!{6Qz&P-HP34H6d$Q=eVW4|k<b(fZ3s^}%7LM3-ifcD#U$UkJ_D~l4Cy2N{VW~%q
z&{k3i+`Fg=4q261+bb&b|C!~BgN`A;74}qM1MTJq&wh0f7cT!RAdai$)7{T&zVA&C
z4r}R}tTE4ouuX$$X2YaF!gicvEU+%b9Q`4amvAni<u%*@Ox^MwXelp-N!uCmhT5?F
zkQ*bEj9reKYC%`p4zUj?4ZWa00!^m}teICVi)|EoT>E-a!>NKooM}UXcBW&oEmbSl
z^#=}RIqu70eu6$UlR(EonUa=b(7@g(HiWJyNtO5~<AQ{i>Tcq4zgTk0wnY)<mWkSO
za+=Ljt`>ILr{jf4g0GF0T-CAjh11YkGT&L@xrzH4NJrCWeitz!C7Je@Z<=R~L#MX?
z>L6I(cCi8bKRR*j+|8i25H%PgJ;VNKFMX71PP(ZL_%hMS6#7B<>lH?80H3mH?ju|>
zNq&X>o_YyHn_yn*ryk2K8uW!aZf#>jC^n!pm=F!2<0oHcR+M_hmZKk<%)1mrJ`@X%
zMkkP3RQRj4uwBPBgMJ}X<Xne?EO563<6Z*{`^Cn^KPYg&dAcmR*JD!tlmJ!(mbVz8
zYz|5xKC|^if{9^K1%Af@2Ua7bq&N^#ka6g7*0q*fx2H=gV9CT~3~>1&xN-ZTTb7Nl
zQEhc2-;iZtdC(O;nD^67Y>^d#tY1wsZXCHAJL&<ikra7;ahwvk7g%fDxQ#Jbwg~JB
zo$iwDjNQiQ6!lyFGBj<BBa_rCf-ER6$IY5s&=@m9vO_Xf(MgRZNUg$o)!6%#a=GoY
zi6J?zKN-`ot!@I1(pye0ZX-QuM3l@HIO7NU5aN%=%3igm>XcAl4F?uwjg+RED=t(R
z*Y5P#K}rldKH*=YDangB1W;ARSN8a}zf6(*enEy(#5{hzHA)sUq4g)7(fKT~nO300
zlia)*U|noQ;F9bYi$#JP+T5@PdgBVb3f8#q8Y`oog1p$~)vU*ALuP~z4wuaujmGOA
ze{KP64*0gbD&BiIKD270Ug$H4KJ9!t41>Ja6KGIe5O}k~uivpdd{GM_hQRHPaRZb$
zu_daVTKtlomId{uF&vkXhdx%z7DW{mep&<36UK0KmP482?5g$nVbw&17i=Y#53Px;
zEceb!$`q)=IH?+g^P)I#01V(!@4e#tfrQo37IYiL6fSB8UQYga8{txb&JE-K8E9NQ
z9HL$_`<-!4p(xYHn$`w(Eg(zagsCLhMvd%cl=_mU8Xyj)!!OdAgz*mwoK@#MTD(bh
zuKhjbauT(nzj~yK00y+RBFLnr67GchM)#d<=K_a*EJ@Bkv!`WIi~M$vEl;@|thtp-
zdgyC)AeE@t7T2se8e(F2ph%-DCyrc$Lbt(NLgoVkaY{A@LI@RngSN36J4@51t?HWc
z_*UN#yK%AN2x!@NPVy8X*{qoKg{H_6Tkw>SOeAOmoMFlDoEsyjS53s9Nmjt}`6V$=
zYUaOM9fo``D-JXd7*58J#?(JQU521Oq1JD9$L}ST`e8&|Lb|jR%LBzucir|rvh<$h
z3T47hQ4^F@qOGy2TZA_uhsCYt@hCqoCIz`af#x^-^zJ>A`R<OvSbV%lSGdtQyAwYG
z$M?(A-Kww^BWBWQ1T0s*#U#VFMyh;#E1wi>U9`uxyy}0clPufyb`w4Z!-%QZAZ`v8
z7x*b}jwHn5g~dFA&KZ><xH-gCHnxt}C3LC`DaWr^eTc{JV45Zf<qfGDvNe}Ai&xL%
z6ru!iDj1#ROo!Vz%PTX!6tdSi;#oW<MZaXTVq<?yF!nV7-!E1h<J3hK?|NH~?K3!L
zth&Gd`0y3{i49H@3JST*v6bTmb50Xg;QsJv+-5MuW;vJj$F<K`ZKoJs>(yxyt`kiN
zBO3yjjwoy+<`RCC3mTOqbt|?&=~=*+gQaA0P%(Zqn-HS1?Zo!~wD;a&O=Vr*FvOrJ
zDxefaF$mI2L{WMYN0BDIBP0kS0s<m6ND^YDqf(U;hAKr60cjy&6huJ?ML?0_fDoxM
z1PBlkzJrc)&wW4hlzE@$E#Gy`KU{}%cGg~dt>0RE?X%a}>!P!)$iBd#!uruDU{yn~
zB+A8Gml4UMJLrdl;%2Pir0EhTv?ch?;p`dHB8U}^(rP_G9y?vMoY0baXpAXsgQAaF
z#vNes?=`QpgqIVEbvqiD$LkNrG-=^WXWrZ}8zAMZysxwRWqEf|TM#ban&}*!+ExC_
zhIni2?U^Oa?b2(g1RanbxYQs&Sll_}Ar0pV?`EfD(U@O$3XJkwWn9K}Gh*9bhAZXy
z4+)t4tRbUSrgrEP4`_?P*MWzW^ZTqO0FSuD)^CYf7S$av3)YIJsQ-+bgOEnAYP>31
zFtR}vwHXTT%a5r0Wt^S5%rk7}iR&iiSkrUIi_2NiiD3=B%4TB=`ROIe>=|QETL`Ci
z8DC~ik&fOr?6In0u`UM_f?fbcpYCW)C||%@EtOP_axN26BFb=2v}Uv{>!w=?0YlZj
z$FR_Iw43s<D*ldLOFxOR0w~f71K|dO`)(upgC1YYf0TY(R2B9Ua%qgwFs4WUnXSi)
z&l<NW!{~Xw_fsodvm*^u6=TE1QFE`{4?W`d#DO^B3gz#9V$w0H&dkKjEd`+#bz^6k
z&KGdohpWy(or>v6(^rwOlB~9wthJ>9jBZdS8P1<lA&}gPIZKhRr#bd`?jy?<)cVtN
z=t-A$sxOTB6GmU-#kBCxsJbQC*uBmhiQ-3lSuKNR0G!L$ez9X&^=sUO?aCDk>6*#e
z=-%<A{@KyEw1K5^W}^^v&~MSH+?t{spe-&xf3~5%J3)ogO}Vw-`fxU*XhyrBp0ngs
z&c4X(d`!Q7;lds-QC5lTwM9HapeyLaWdV3R<$PL188b}*DI|+Ubz|NWO}9<kvE?E4
zl=UHGE1)8MZs%55_0bRep!JYm<}Zsx+W4z#D>8PDG(^Q<7I<TSYC@&L{4Q833<QLQ
zZy2*;NTLPImqO~lqV?u5TG1zt%Ya_)iG<T?&wHHdWkhE&m%?dTGs@@8={c)wL;wYk
zD<jHgYRO+TQLM-jSu5f#sz^YLEs{O;lv0u*gxBI-nE`^y*wMYc)CY-~ea)1m)5y7G
z$3n+-TyiT2DPevZV9}`h)(Yq8h-fb+ZTaCFQ<FKk7tLx(KN^t<ersYoP#v9Yr`vLB
z9euPVD(k$E{$ZQz`O~RFcu7!DgRdGY-ygtp^a#8&YHv03tiwccYKYU4fnb`*zU8GH
z`dhjpgelt6da=-T$Zs(m0!y_}pSDBybnmzdyR@8|&{I0dsa>Em4CfR%43S2N$Ozna
zBMgT?FBJkV%LEX$&jnRgn~X28xA2-TjFhL8*qGb>$5j{?4vFInGDV;07fmsQb*YNw
z6voTj*=tUsym?1lE>zI#^q6i&K@gEFJ0V#w>XSeewY{%1HBmo19t7ZTPKu#j5;<|$
zce?iEWGtc_#G#@2#9Nr;tWo($^c*zIcPhwb@NgFwcZ=CNCnn2oaR@U1QmNEtsXCBg
zSJy)vF9zfz4qn7wUB^5L>l~m$4?lV7m?j0X*}g)iCJ!*s#Ek{`sAV7c%Qn2&b~%`u
z`@}fQHLR0Ki@8k|>0{0%%j#vn0o8NHS*Ym5H=V_92kV&k`%FC(<LnQBI9dsm)ioWI
ze=ao@aVfr&Ileo@y}h3{9RUKLm2<f;yRx36@u4z6VqSaxI@C+_5kH%iBL*nt0x!`0
zI46f1jzPuJApGNjM_;A!*{2c4WJn3KV`o_L)1D5n@@tmhM@reDkwvh=chr@`OIw>`
zL=TT%<aVa1nl55Qb;(~a6S1M5&vtpAMu0?+&7s>Vbo#0uJ%`yc6hIo3w;|^Ni-4vV
zW3F0cD*@cR^sIU`U;fVyb`L^`C8yrpIK@O=)3WxfjUbOCB{1E5;oZzPVx2D*1mvUz
zMn$cNa=RMk%92v1Gun97fS{ASe1D&M*gRoY$(m($I;G>7K$-=@|Gtg<tC(eBc18><
z<5-FA1h3jD33byKk|t9Zsd}-H<hw~=d=fZ(-QGStmY==u|K7tI<!NMFD$EVcbT#Gz
zIM_L+TNu`2izQg#atpzmPO#epbrR=2!_8Md2FnFzxuDSlR*mnup#qn)5#2L7^T*e+
z@|V-MK6;wB1kx){SW*r>q3p7D=-hkM=f~DY-^6Aq@bu~H6<y=`mzElHrK<W(G`_Z0
z=6zXh@n#D1^w8t;eRfDsAUaHxG}q_GIV(bNDk$3BDXQ~yt-)k$Qv~?fxSs#sS2m5@
zCO3Tr_bK9uGrA)Qo~Ul-t5;sa=e7c$Cwy*+if^TcH{uSCU3yxX8v>q|&ZbdYL%RX{
z6lD@o!9mfj&*mN{G-&`|e#r|%!{F!)*^VKBeGeiyBX(j0)Qgpha0J6*W6N|F&_yJI
z2zVdV5Xp7R{$x|ZeP-Mn%-RK!f&+|tem)9I!A>iy#Lbp>r!Ald%Qz8udB$wg3hQ_u
zW6s(%8nKsWeaZL2eL4#ohEhsmt(mA@FyR#vHrHtmeK=cXjqb)izKhraqQlNnSc&qb
zk4{K-rAOSpL}wT^!c!QU(kwH-k*eath_VpA-UPU8%3j$eC*{hs*>LWb7|V|py#!do
zRC$P_RM$ED3M+B3@$*x%T>PUtFWcJ)P41S{>mzBT)LAU8WIA~!bC6ffN&>4=L}8+_
zo|93RY%tsHFN=@2utx?+T7{~h(<1v8G56N&#KvWYQg@;2cHTuO114&}9>Tmf-ztQv
zwm+D7Cff{H*R@hHY%CqdD%1p-h(r=*JcjPBu6RA;@jfkKEHQF-iuL-iH0)WkoS6f{
zUyzx46)v+^+%8K9FT^9Xjm|!bDI@s9Q*s2o&o=Sd%SDb-`3ZK#Sz|^(U|}f|>HjQ1
z{}e4zP*xr#Tsq7rcW8{kxmb8-fdoSXK-_3B_%*=rBPsU+x1Q(g2NDs8^9UwtQVj@h
z`muR{)05xPOckAj>XP7Rb@D@()Nmb~`sYtCw#30l!Ub(3o5HwpY(-!m?OV}2_j+b!
zEi3dkA`ArE?V9n0(Y9z{uL`n2(0^7t#0#}PL|_GwT4^5vxAz7V&Hu!Lc9Xs?DQF3<
z=Mhp@MZRFUwuRJWgYkhpuEa0b*MDVy4%D;P1bqOB>~5jay>Z(oqGj_x==y{4lYC!D
z#uR0vsDw$vb=&t15o4!mii)Gx6ROT+qd_7mLd=Jg%J+p2chy9U!pEp8>zo0S&s26D
z$RE|sa{D#F2x8VD?BdQF$am3s1BQHetb_!Z$i1A1+j=vaD0%vWvA(89jS&f3CWK{0
z#2ydT$R;47LDAdzViL&E4;LP1DaynkEJ4v*_)7W%V%rj5xK!{82|YE3Vj(i0qtb6g
zSR4O@pnrmgbi)$CpqMk+*@#$C_91CD>-u^Y{l*PzV+#Z;sz+#-{c(^UkCJsP(_Qa?
zB|rWac#LgtUwL~p;OB)~a`p<KI#4v<@ZBkfQ8G0c!2{kS6iJh6Ae5ZDcSBn0o?<pe
zPo8D0vR`<Qh<W&m+J!wS@rcxjv6n9d^A%`b2yW2!9V?+vU9y>pY&tl~{usy^Ag5dN
z1#=8#u@tU}23rC_9(n8xft=^kBqsR9M&gxmxp>!ZJQ!R>`2x&Y?%S0zr6s#6BLG)@
zjw0SIfZqyyJFG{TVbB`YR<2Z5VK-GpMgBDOm^jMuOh9DvgdeoGk8P`ilrTQH&))>#
zm&2;Kw>!*fO+VWdq-r_@#0{*ODv`<Prl&kYveCI05-6d(?@l4;zXTYIcZviYtOi7d
zg}i-i#xp*Eka*R8(c_uzP0iza;5*9M8r`)b$g3!sZ<tXMT3wkcGu#g1{W5y0pY+qv
zRTGIW<?Pz0R}DtZY|$ZnMgj-`dTbs+u)u2(^Ng=(_veE80{6K5CAV9G7l3a)?Xfmq
z68B2E{2$B~dHjmuU16!t<U$`wb_D`oWYb6kql-b6FmM}hGfS1N9e&sb{gCuH2B82N
zL`5cAcWTVc>LKVCepV&Wm^lXN@4!k*_M=k9KI~Fa{4tO};a4<y=YgQT!d<L4BK_qn
zhDkTwrJkPRZn=Fs*T}}96_{})+`0K;1_ELh<kO!9gBQS|a80wwao$qVE=+3Al78rI
zE|7!&g)YD-3nFeyao81aw;n`DT;S&OmPPI%@!E%Za^rB5jC}1pZqP1~r%z9g<^3Z2
zFn>uod*##h6pU5RDG&&-E5~FrmnkJdd=FrBe9FntVZmJj27qXOVcUab)1_@?vTH3I
z9a;l#A#Q?|O7&RAH&b3&3hYbcq8~fzsO5VvVxF9I^0KEA55A2@Y4I_Q;mtLszbB8k
z$}#$GsgM*{;ERYk(Ox4jJj7&E<RA|A9I#^A)J5n8+|=`*kLYqkwgGFdZ{8rEr>gw%
z&ztceH@Kfq&a`T@6#haHEbtZ793^@<0Z|h%)?7GzLjKgZw|?2R5x0wtVrueaKvRHv
zpd>;@XPrlWBcGh_BFCPl+r0<y1#-xwlem^r0sIKhd%LijT=y_gvvC+TLzoH89wofR
zIBcXWb~$nL8CwkwB|kd`4ioKyy<pqZ?o^#KYkKiO7(TWAsVe{oCF^cZPB_*Wz(iTZ
zvM?XGtx)krn;<YcCt!4#f(99BUI__qKBy2}K~|3o@5|j1)iW=aBoPrYc9V6qj0YUe
z*IfJ1j=1ZDpn*-(ix<M`MOsgN!O^?|vU*966YBvpqCndU=U=P|rWr<PJT$!P#;XIk
zM^`f?nY&z|0|HOY3u-?ggrfxvkO(-WfGTQ?1akrIaUGTuN}m{3-jvGNEpIOT6>bP*
zO&3GbY1kGoaW5-^%wM_43;5dGJXM5dD(8iAQ6%7nu{cAs2w^F6q{RI3iS!;2{uQ4b
zENRd-4t0$0OtuHWzGWl)k}lKW^!FF6{d|{k7E#}I<EY$n7mw!xJW#W-ilDGSBoKH_
z1*C;Whyt3GcKP{M=v2<nTzD1k;RN});!R(}w^D09vghQ5sm%$QD5pFaAI6&pM2FqY
z#?c;x6HEn{#*=yCfKRf&;eF{1(5Q<oN|O>5&OHPQhui@MwP`f+y^;02M`W@!46l4p
z62+6ii|ekKrJvluJpwKS!*QL;HZQPJ;7bb4Zrl*sA=1Uksrpdj9>lvSlf7a-cph5y
zlf*7yT=^AFwowRa(3-AMa^P86UNvh3E2ig;N}~H=^P-5`Ba!Jw4pXH+gGBoEh-!=A
z%2E<I+)2*wJhM+ZPVjP7#O=P8K;PWd_A!CSXPOkIZkq~RejagqVhVfCv(m0(THuQr
zV$^KDT+~tO-kpg2igQ3`^4vIbJ~4Uvb{lV=aueHrkV@!F-P;Arwfz>U3D`ps8KpI!
zfbYRU(R=e9mgMv$^Z;jIi6Nj?g-n2n01cQYVX}H$1DJE)6FZBAq4Hu+FM->30QgSv
zF~b1f+C2O+!SF&ysURN91DxhUz_N06fn~{Szz!E=N*|t(Fi(!S9bZMZ?t`)=%ssft
zxD{d_+W)OM{_u7=-=y81_IJM7)#b+#%X{<JVC02e&L`V{%7${c++FJ9ixV_Bji5hx
z<>ep^>f;qq@6+fp76XN04c)3s_uxmsRpkQE#?|f}oA2XYUIK0~@Ez3#h|1bPFcl^D
zo&E6$z`e!w+~|0G2hju4p~+&_9KwNHQNpEWCxBSQocOwF@djY{^<br>uj@T}r@287
zK`H24TX-e1+07sxc*5Hgc!WQ<A~N|A?-8iW8$Pax-~Uz2rtZss^V8F-;2pg7k53bY
znsSVu0W!lNJ>M7sccTHP5M7fEc#0kN_Vqo1(m?F63Zi4L_+p&_ldg~8;hJx^B=50@
zgLpT;5<PIPW-%He#j`?kmoig~DK<b*6pa+8Q+x#seu}6y98;@&eTYX0Vsw4>&L!R>
z(k_NiQd%~{$z|dE{`2Yp(B!#cqh`~AmG5CZN>u}}UXKtp{p|YyN8X>O?S()$h0WJ@
z1`CDhXOC=Zd<nU$QXlx>4+aaz9?k7h1TBe3m`-CjX}SPA5nSK@$n$=wYblg3%>)A1
zl8D=PqSKMApT$A)ymF>Ww_Bc_{kGfi^Y-^$XI&-0(~8-wqxnnD+mGu5pFU6C`HBun
zQv&f;@$BUJZqHjKni%a|=2Or7LCuEart*AZ*?><+IMw@>Mw|$CSK5i6=6UZ?te)RG
zxn0imO8exKIN`%b<+xiqveybx;71ZJCiY-&gv)D?$SXI({$n~X`;34Yb(G*#^5(b$
z`0SDF+Wp;Mcg7%|fJ<#F-JkKso^5*ZbQYRR_jnh^C!oFr!@b37fKFyy*vEr+<a<hZ
zx3;LddYI3CUqbCSYK{>Ulow_#ZT%jAt@rm{x$XNlKs8`0wPKwUE4h@+yD00zK2dn$
zyLik5x-e-yvImdc<N_VvA)tPI;i>>W2o%nFaXSJ*0%ned63Z?4&pb-JS0O3u`CN$Y
z+<ZY%>3c3)o0vAeh)$N-)^tNq?r=8yjAvAVwVs3~H=m*m@`3h}(mnvw`=Ys6?dECw
z^DY6??tL3%<%Fv9H-1zFyQnO;W8wvPAc8Lysl51rH`X?Sn$o>?T<eGOTm*n=+LK#x
z_{P@n>I~%jO0uunE)6(332d~H1$1+=mql8VcGzFH*qk|X=76<IG$<UW)Kb~`h$ZZT
zo3gvNEWcd<Dxc}oKUf0H7QiM85RD0eBNK%>k9n6GHMPc2>`v@vOLTz&i#>5cwUqmt
zW#3Ud_)@}L5xCdmXp;3t@~9;DnERQpJI@<)mrRsCR0HJ-0-iEZt<uFBSiT{M?`FbP
zKX&uJzkJ8t8tf(CM6zt}eU+!Lgcyol0Ur&20v3WAd7nVb#5CRGiA+3v^Y}qvnD-*N
zSjLjFkKC7~n<zA9i@w`;lzrcz+lS|(u%r+kuLl7rCH4DqpW77aKkwb;a8a1|ZvYd>
zcepp%^IpZB6e+|T#@0o|t*)mX0*xPOA{{gJxLNm3FQT@n{h@1Cu8@T3*&@pOpuG~v
z>EKssG!4Kpjb9Z>+XArErNh0PI+<x@Yn>U9z)vXmJree9%h13g0wcXC)%G3sm$O4H
z9G^p`rGD7hkM1{34p<d8LM{bm=mebIj5qdO=bk(62DgxjY}&^w0Cl^%Ck^zQf{$$a
z^qxXO%fHZ3vyPy$?_}Ah0KwW1%0L{oIOG~|ph_;~MpFc^xY=^(&b;8JEcOVxpbS5s
zseEe$j!ZS$HwFF|8R07gK562VOt&OHei$)!*0ZC468P<ho4bs!c<tS6U;Ok#sQ9q$
z0KytOMJx3H%F{AUFsqv5ufl&c{};-e|GGHwtPHt3QFbymrS#}<6RbCu@mk@I6oR<3
z^`XTN!rahJzrxV~)`wq2B#`wca;KX31Z0u7bT&tbfP-cg*<ZRvHYrM&MPOt$W7gk2
z7~on}Ksw~u_J=Nqg7!%GP6s|uQ~X^rzscx0-`5OyCBHtf0m(&Q2;GJJCRA(bY_^`K
z4ZDZ`d;4C{EYxWmuaKw&GWxv5-f#V(0L6|*it!pmgx&|sf8!?^XlCQpBl_(}tQ`aY
z(8F(69l}p1RpOEl{{TR>z^8yT#|#;nD4@YR+Xzg-^TPP3Kj`&d1{VeTCi9wz(@zmA
za(F$PO&|l|^X&2bH<y1<;kA}*!q`j42HyBL9e$U_5ve9t+nWq85l1{1-I80<Sq&MX
z;khn<3I7LashJ`uszwKD5F)n_w*fDGJ1y`JK>8=jYW5M#!~@~i930(&hVDs$xBltz
zKQwUUFwn7K#wC%>@%`hr#Jx1YM58jVrhH80qW6eG;=Vy8gNqft<#5c8rm^-BtV%%w
z*PZg)T>h3Zz0F7d2&io3CI(+B_Ybm19l7UQ&(2W3{N4@I&1uU&G&XkId(*6b_X5+B
zO$=k_P9esbKi<wl)o%Sq8vv$5JfXh(na7*gH^m;t_muLEwXS;o4|M+6`j3?UVw1~R
zTx8J#G&M>1cWnGVkRL;gXu8a^PKCnYgCA7CIeUaN&wI1!u*o#_-<$j<yAMFf<MHw8
zzuNR41t9Q6518`8EOK}K__28-7t_5jy+>;E2BM^e^TOtLg#AHx|4a9f5rosrlAr&m
zz(13gi#;eDsj1}o>ie4OdlMiLAk=rtCO-IwL=e@Wp;Vh0L%$z4ejw(*%O&6pSSpv6
zf0g?ol0S?UQOOgEEB!!C?fWwg3GMxadFE>R?mtKq`6uE3j>G>Req1uZE4t(3fhW`W
z{)FYfB9^WOE;zKmL0tKpik|VjM?Le=`1kzDe^<``gL+#Mv0rNsqYys)!D;?MYkx%q
z|M=O_-Jtb#jz$L0<Uj8H^{ab8a6-QL#ytLijfubB;h+5CihTr*>HbnBNpka{wG=Qj
z$4NE*cl2ki*hFcx{@-AH|AV0T-$P2x!DR}XPSonZ^Iufxy=bn{n}!qTDm#yXHu!zV
z9=!+h#a+H7iMs+x;v?USRBdLL{}+?Qxwt@FTz{#Y)L+ad=i-tg{zyKT_WquH?p^dR
z$>-h|z9*l1UjUNNZ%O?~KHpV)(cPKrYn@W<C(*5zUUKS=7goC_+-@VTqRjmvJmjO6
z_XW<Mc_1m^llc0?xkD{};c6$8a~%ZC?JO)H0@my24hbvmR6JfOG0q+S({Z)*h5(Mo
zhY$_z<012n#{)Y3(4R7y3?hw|u|cA~ecPC~Eo)<{|Jl32sRK{uv^N?{t97a!8)`0T
zZ3J?ZysLwqNb~Q4X3Vft^7n(ML$3JL3^S;MfymMFYUjdXpG-R@jX0h3<U<!w=4mAS
zOyTFne$&j2`LObCN9|I$J9`T19TaU>^_*74B6lRSV4SgZ8t!u^^gT7Q0nVJP<DNu@
zcrE)>JI-ir_;L~tiIs=-299oH<sRmLKHI-fYZ!ZlIBTf5d%5i}f4>8Y?$MD+pC<LV
z*A$q9_>>u>`;KZghW8z?)mh(Fldq_PW_tM5QP*Zh)^2nx-CvGPacnR5nd{QMN*(g$
z^lC=6!@tIK5~v3QT&P4P9BoMY>Fj!cd{N(E?dojcBjuTl(#!<dQbuY<*y4T<PN#9a
z*76&Q6RUnNzlLXDK)m)uM(RQA%5lp2ODhVpZlTFWaX`PD?vaV<Y}>*-Po}M5mOprG
z?B$1f5-@gIp`Q#B_7~`@!It0l+Gd4z8M1u(m`ye_Ot}A&VP?1X-i*qjpBXe{%ZzD{
z^Chp2gM)*f4jyI4w{@0|uG7{YVmy}9F+C5VxbFRhS>c^mRKHBFrre%kK6E%y>E+-t
z$VzoTFW#u+OwRuW5~QU(Q|fORznJDOe!oM1`e4A)z1|_egln9Z;>pPNXPwiK-cwoy
zH$&gvvGty7?5P_#p83p{I$C|e@p$GJ0_#9_pJ&6BxjXI~4#cbv!Q_e2q59D+Fx?Q4
zw%Zj^cC`O<RKo=-dH#_1z(7a!L+#~r;J0ssMmbZ%C9rFs-KSFLwdUvb%g*1AA7SX)
z;rv@zW$dPpFh)|j4s}~m&<j22=(rvBk>itzCz)Eq=bv=8s|S}+JAM&Wy-Lm*6HFif
zRNEPZovfG5y7p>+=IM^Uw<nCB^xssk3O0zbe0Co8b6swy_%nvO50&t+KR9zB0M|QT
zU$=<F^_H!^YnbE^wddYP!&^{kA91C9gPfO{FM8t#S>83HtEHkFqvygt@*D`k_9vD*
z1+*mJ#mzsfq?jB>jm#8^uFtFYXKZA9=rD_oakZCvJ^YB1=ONgSlX|UJKN@ZGOC#Yp
z-hP}BXVHAzdJj9Kfqt*1kjz>Q8R4v@`<adEaGsDue%&ZFMqO_j2wCrnSL{jcUw5I|
z4WZrF)(gDLIpYI5O3us}yR4-<#!Rz;p=a$ADC}Z-PvL08$J|a0SNFBgTg)1%e16)R
zm~|>@Bhw19#8|BGz|n}#gzKl*z6`hegmvx`EldcX8PfI`fNA@zae#tX%$vPE1^TmJ
z+3c)f-?c6{Zey>R_UD<Yn|7mBmnhp{P-lPk5*+>N@V2~mT4>+#OPp)Y_3K@y=M6?V
zcLP6YCk#z}a_cW_KlJs(0qE+yRPVNtLDF*Uw(6d3%WECXaC&0-xm4`bX#$R7)0$_D
z3jU~BzVZI>R1ZdFWvF8#$iMzGydO5^!B0Bue62<5<odEo<~oGaHrGI2Uk-d@)rt9q
zFi*b~_+}nQ|H7|3H&#WT|A<_vCLkZ499f}N){HP7SJB%Vv(ioNQ%{n{-@65)Y1?{V
z!AYZZ9FshdA2)L-aA<b>w){A=QO5Ow%o^jt(3R8unAN;dqXSG|6eBa;wjhqY{$^m7
zo(n06Ysaza1MpuL-f;*s>Y1D*wk7Lsxt33Fg5SA0+M(wQzrr~+??_WAHS5gstc*bv
zdMyaG_G+|Xxta>y_|U3euq=xlA}!%I6sqAX1i$*_1~Y@XxZ`(q7bkH>bIe9sDbBGm
z#EZ-*J{`Qmp0Q*-(Qlk6=$$v5<Fng%v%iL@RJ}2K!;kajMM%N&oxX&%m70;{LF%DJ
z_3J~yVP51hV(VzoYjVNz2Y&4I?HXjpAejk;IB?Fl7A#MNIED0DX-&oO>wdDL;cUFf
zU)r2Q8rg4(Fbi+gXV<4q8$=IsKH<>oJVk7xZ1@}>1mojH_}VkNGF|6qHaA{Qx-^jm
zlo#_5|C%TT!PtAj*?e}s8>tY(Ik$n0*0<wW^=d+WaBsQ_`{q03^{>ZA7jth~vC>Z4
zWszpgELnyW()7cSMz!h$Y=<_Po+4Tpy!MfyLLAR|hpZg<<!)wxr56O<O~6fvw-zo7
zHV`+g+O=mi*2@!c!{yF|!Mox?YzDc`!wV8deF??jT4(Mls=R04y<4z6j`SZYq4?%6
zPrJdFT3BScYJOz3D^qtQV{qFXwl$YoCtb}yG<V~sX(-wY^7W~eN<AHFdA%?CIyT(N
z<4}U`y0&`85^(;&TgJ#rU!$MF99zB0jJUGDsB!deIcz<urf`|Gf$i6YvVvuu={@||
zNdhA@{)tYHa`!&mwnIV7&@v>WR`hwi>GkkyebsP3$kKzm&TJ+>#+7P1+o6BOb58ir
zj%Z3$5TGrYOs~I<0rC`@-S<+$0tRNo{O0Rw7H*3QnB9DMP3*0wrK2V)@tVxr+)mI~
zX{X1);P>OM@GSz5>v6>js94Q?iJM2a$)K7@M^~w@cg{#jj7$fp%v3$_B&!r+I6?I(
zrS4PUt&Q`vrH?`>zSrW<xV!*wjZ1BLv|O7&wx2oX+sJdL(Ki>9x71&u1zHjiC|}Y@
zNPHH3o0KIIMl|$4#=m)9-@ft)wmXp%4|4#;2=m9Bn7=snu)}qG?9AH!rP|zOl#<!_
z9(x1)g>uHJkNktahX-7hKzzzazXWB}#DwRmmg_(!dLNjw)(czerWkL_D|N{4VL|W_
zSJr;qN6<?o3eB=s4Ie4X`Owt@5bxUCjc+;$eif7%KU4cLJFQp&n~1Jo%S^>Cv29Oy
zZD7as*_Pt6@v@n|_~IF9;r3}eNJ&EQ{j+Ert;-F?9Cq7XadI)n22CC>k-drntu6_z
z8qi7-j)B0>MFgPdufOg)FHy$S>97-DLZd(D!g`2Ah0$`IX(=GDaKV7pr>`xKZxJ$R
z`Lw!}Ms?eLGsUjDM63|SD#lu<NIzA^C)wCK`%N7_uO5p<vsPBwL>-F$hR*Ps9D_2^
zI+L0prPpAPOZ4>Y9e~fSnqstWdE}v4I*^xoI?1b{_w!YKZHJ;4iq-t!mthdUF6`h5
zD+%HDc(pPDRt(+nT3?JU1sT^LIs+S=nCOO1OA!h!AfE5^u|N0jgPV!2X}n6mg10^g
zExVcK128{N0ctl8bq4h@c@{`vNVf_co0Sl3<5>k!i}TOHWPKk|gP~!e4=Hghg(6k<
zQ7$2TWzthYY`IyVp0Rh%EyIAZa16Dm*cwNd!dzzR#aao~SiaLgb^lc@`q7D*I&%Ar
zsKTh<A-<r4<4UoJjPliZ5w?{o;^An;`FmqwBaHyh6WzYgA18$84Ua95b^X(-jrH0u
zi>_Xa7}L|FT!whYrnZDsC8$BnNDl#h0M$G$pTXqOnPCq_c3Cy+IpHi0`S0<pYM+T&
zSkSm0eBg6y)LOCHB*fD{?b=3Uc<h1kGw7_d4a9KLDz^1N@K88W2R<7*$S-qG;n4)g
zpPhF%Nc#-=8Y*hRV4PF@x;4}@c!o3VR6->cUm1n&j^TF2c0(cR3+x^uZiG|}6W&!a
z&G4U0N6pfy+U51k8EV#!h+}axC&vy?yKO`+2oFbJ8;!!aKbNBnioUFp-0g<+XBXz2
z+MlYL5M6b!{BZ(ZGL9|@nPaGxCt%SH^%IPs0Uu0mnQXYENlSz~Y&pt6N@a{BKX#$q
zw>2qMAF~(}`et0=6*L|~?+t-DI0+43UoMd=H*#yVSz@=bmT18DAmK}i()S3I1scQk
z5#8zDm=@~869Yk6!EG$jP4T|gC!P`W7h0-7oxyx=7(PEij|qQUqhQiqZooLQHZgIZ
zW(S#D2qdu={6DuUcF}asViXD-CdX;)+IY8(&?_|;R(1AFutX*V)dfZGcZw{jB9lU5
z6vrj!DXvKLqi7BqomK8-7f3_nB-*x_jECNQ4yhitg!sK%?a#;6>Wd9{tRCFM(w|in
zgrLO+Hm=we3?(eo#sqhhBb1o;-=ieUBq8G$zpBzoYh}};ADsF~om&V$m|d(hoxijv
zawU?1MArDDf?Prx<$K<vT;XqeBTFJ<KIGu?{D(a}(d)hp+qRy^&}c>GO`<b^7Qi6e
zrO4)<2W5M1iw$TPE8ZiJtq+pIeG^#4k%@i7Tc`Z>#8{?zg0hu*m{~px2;JgcD_sg{
zRdas7tlsmt`VCnNU+jS7b+`4`#-g^_@+6y?*cny5j5A`abDyQ+Fw$A9agKaKLS$kN
ztkY2enZ$BP0UQaC-kYa7T56xtL<AR)Np6i<vldFCjh_u~#@Nlu0<yd0FEHRYYc{b`
zLHX?IKXh5t^KIc{l|;FOV0XP`$wBv)XJC;t5@qCtRoHS|Eq@IANS6MXrtb26xrETj
z@OPq>=ky#8LDMlw`j2c<DT-)J7KQa@uRSlmwzCa$xQr2*+Bck=_ufO&_fpJ)q;hK%
z#~&raGHVIY`%td!w4eXm*yK|8>?=f=5f64{nrJ;XkV4vutiD&5RF;7Z1hp+#9P)f8
z+R>I*OuTaRQ;bGuM8^rKrw-ZE^-&Af8p#A|D@TxdUQnG*+(Umrg8<1k2`!{O6KIkg
zs`yBMgLP{_2Ay~3^(<C(!djAY;SkCd5~VgKJ$7WiW^4_*oO-2}#I#UND8`~!^FEhM
z>dG^`*K|^**ONxK!~kZah)QW?x$6->GmKiIq8D_Oi-BUtYLK}+u*@k5;Y2wC)eRm+
zO^HnW;QKQ2<^rawc0=LSezhq7tv==AfbabhqI(vzs(E2nk-q+*)mZgItzQfWE#_UU
zkg(QbI=<qc0Zb{DpqIQ!l+}uaW>zaO4vM^<m>aW!bLUC97yJtPt(lRFD!7d{-Wm|z
zY)mxM`KasVRO6YI7|CJG!mifjYjrm5?RhE!n!voj5~3*&BgC>Q48~zk#xPGEUF^a5
za}u19j@|{V(g-}hHcBIk;t)+;yYha4UA9lXJRzLCc3D>Xs40G3D=s{50=AdFhqBs|
z*Bg<TK;L?$Y(G9q<Z*pgJ6oD$x1VLCLltEiEmC@m_4FHx!zRE5hS}>QSRHCBRmFe9
zMiAOtbv8iyXfa~F5kkmKf{ThUFH{ao<yr7sS8{#@wo4sMA$Sei1jc@u<c~ZwS8_+Z
zn&jm)0Z21@$Zm|JIRA2tq;wr=h>;nFBG&@#Ok8rn3p1KPdh&}2X44pJtp&SxPR2`B
z=Cn%`kL2z)DmxUe1S7m%PL<BedMTL3i%;svi3Ly>UG>Z$hE-Ksr-VPdJ&9#hflFdu
zpzF5JuI5{4h!{)kjSP*d3Wq4gq&}wfUa8fiyrPS;5xWsC2S5`nw^iyY9XOi4W;*PN
z!sx5RHiBu6c<Z;sb=bA$DC&BVZlMh6E2wDSvE{pfWx#N6rbV$^&mK9d3o1(US|!;x
zzqwK)IURFp4=7G3l6^!3MS_YDYMNQDElkB+=REWV0Qt)|5J1ayCwn5?wjvbRE6J_x
z2l)@DvW%Mz!w)4sfem^9)+Z)y&cVjavRv_9P}8-(R*iAf))N$e@sOdBAZ_ewUHg>~
z_X$d)5X-o9V0HI)dAw~&IIc95WX5BZmn#-kW5KfoDMVn|DNHxIrn{hew@QzH=ZA{z
zvFEdikQhw<${}2{e%^Ri;rM3V4B+TCZ0X7fHi;4vkKWLMkNPi39o>zu#RBTb?yjqa
z#H{(h?NHr*t{Bm*S4Upa!C9ImqfH_Z`=NA3>zF&JGg=;aW=0fZx+Q9%7?J?2fJ#=^
zH)=;SBs<rtt<17_gZvu{T~#L-#>;hB@qVveU=gG2^^9xJwOyMnVhPk-C!_2rkRX2r
zS*^S2(e|AR=A5DSD>eGflx%jjRM!#wd%Ly`xMiY<-xKg1AvtF1N?qz<mbYoPA!s>S
zCez1T+(rU$0HbDyf)MQO@9eXu8AcTL3D{K1$tpziX*B$#_M(86z~xv(UAgtTj=P(L
zQc{xusDiw+Th9s<w^xqV8a?W^e>LJY1Bj7;{&aR*(pjqu3o;Mw74fHdh7xuu<U1$w
zSXQh=roLAH&^)-^dV5|@&N}Rz-f5<%@X_c4(&n!W=>e=bM`FvKedZF&F~^V!SaM6m
zOV6{}f}ja%a%G%-nckj#hS}krP#A8%=h4jV@0_xSQ0tt$M+DFI1A5uv*eDuJl)NoC
zNp(ruJOpp8baeS)efhu|Z>*#_32wV{535toZQ3}y3^c(Yu`24bTxT2+c)&vka^F2E
zgdg0VL>}m&CRc0Mh*W(Atul-c{t|G`z+n3o!N^cS$oX?b$$NZf!t=%g*h%X5WD(i=
z63eMsJ9vPq54_Ek!-ZL5z~yHU{ranMj#7E5EW!my0YpZpD$5`RU-MX6+*@~YH<dEW
zwrpZxrQ>!0H7XR4dc<0Yq9;+p94I7^Y9uBM4g*UfXQ8@kJ9v{of)GKd`)yG$ul0^*
z=86<n;XtuOmkT~JVFg{W)yiHAkLMX<u<_cL4YGd$75(B}69{^_qZ#9n_i}k43&^mP
zR2kER_n|BL&${5jlD?A8*%x-5&Q1Xp+1sp>VjNX^a(I>QArxNG(GwuNSOH*35D|)M
zHL|otaG!;QFuSEa@?O3APZF?!gwQB<fnZFH0qFFOyv-s|drTQ;cf7X~Lk}v4Wq~Tn
z&kVf{*f-3p_LIceku0|bJ$;>ti3RW@xqH!i`h#b-c*`7pjF^RZcI_S%l;|q1hnzkg
z(IkSen#c!=`YsAa0=}1}W@nP!4HA?s>fRw@=JFBnD6jP@tL7%XwhMgWZKl;3r~!o!
z6z{D)?~35U3q&&Sk88@b0Y{`S?#LUFlm>*v+q}`5ljEq$4Px#wXqgEm&mwJ*dnE2{
z&y$3B+A14IG&O?9&>KiS;HLRu0a}cnWA{0rwK%u`*fEKtwFn(*O!_Y!RadeQ&H4%{
zhgK&I?N#x=@br8#sXo$W50qFN!<JSEh?u$n8ED8u#}IBNylPGoW70i$0_<UUJlL|^
zXlQ41zQs$PX2X|Iu6Bq~6N6*!B&P1v#B3GAa@%tdfQV&~TgQ4Ln=aua4;%!z7ZrfA
z<EnuAA!E`<^1#KwPpjj^KsKNev$6*krrCDD*1LDkEktqzJKyN6CSDni%K$P+jy^>U
z-PO~&s=sTynz_VSGuwWP8A+O7+5A;|mu`?mmlHlR-5_1@!o(&;{Y{E>h-TwFbfXS8
z=sdXjjR;>k?@k2?->u0#@l7cdr~^<QMcD;a^xWdKgec(U)i&{)ZHFK(K9W0(00h3u
z_ug*Ouh1g^D3y}xI|1G*63LeC!E1XNWaBrLhWD<y2x16F);{#Zu{XbX5RMUonrM9N
z*C<9PY>f<jM_-gR_W&G`<P(PDlenhGps|bAr`YL7B#@5+pFat*=B6Jj&gPixF)(o~
zt=<FRFYY}m1Z1-3>-DvNZt#^jnpf%b1X!2?7El#<_u7%oMe5~y5!BSJnOBU(9VFs(
zZgBIe={IqDG6H-6gb3Ih;j5Ci1t@{6RPIHuHbB!C_Ij=vOG=og-1>0JOHu;qbxl(F
zCa{>^D|1ig3asF?ix_w))#Rnau~vV+$V7hvcPl{#V9>;)-RXNd0=(7^K<0TtPPr#G
zivTqB0OCm*84@f6M0So=cbt(A=+)S79w-STaQQ~WsNaWMl~>|`fBtbKe3u$GO;&v_
z_B>4vH0Ca(i}^>PgUwb=7q^FY5=t2tP6EGroDSq305zS%M63`*rGl~wplz29{vA>J
zJOJR|$yzCV!c@rX`WNWha{C{*#cl&iARv5rEPd=QT>y)iO3c6Pd9lm$o4h2<cQnJS
zL_$wRAo3$Zi}EZ^W$Odlz9;vCL8YPyR)zFL@Kbr=NR6mKqtpw7;1p{K<*Dw2v74_Q
zrd2+{O6NV^ul@E~!c)+o_m~Ca2mqE%j^0@}u`Zs_tw^cP1~h$2bC+!Z-kP!GaA@;k
ztiHs2YftEI+T-6?2l7qxUh#2B93>)jT6d0GeS7a9(WU5unmS?eZ;!n1I>=2+I6-wf
zFAq>lv25ZzOB#&V15i52n>VS{XS6+5=@a(K#h<qPLEM`lj@oOmS6yna<ZnEx8B4%$
zQB@8lHIADyv{&M0wYk0#O9ckdA!5*hB=zmxXVk=g518;Dj`BZxmYV`~YYw@*IZFE`
zj^%iwx99IF`27OILxlII8f@95oOTy@gL3F{{U$pI6gSHlCAsU5Q233Xf97^mP=Hk$
zbB{twCOoU<fPZ71CE_q!e=imRWSUahmk79$-zr_af&!ynr#9{Ux3pi$eqvY9KWHEO
zP8H!|Dscg*@!<jtWs8+3P=$mWbei`)>P}Nn(x&>hiWqFEd3W!dCIPO1P+Eruy;H2Z
z3DG9bnAi`H{29mp!C?OE5rFfR#-oWru(x-<cXzogg7s`N>$>8P`W=9`^8sa{3g%)T
zLVsibJI4O~-p}s++E~%OzWBSS{&R*ONqxGW+MzjMYxJKH|F_cp{mvtmPgI3;ERQ|>
zJEKSCLG$f06X(*l{d<fboddx_9xflh%=sNT-!1^h&%l;P0Y#~PC;o?vFT!83*M57t
z&$)?nto~=gKRo`=NBjSNV;~0T$c#~ce*d4{;txpvhf2b;M80A;tYC}(qAv0Ed@CIL
zVan5`zgE%Z|Eip?7y|<OSZt&J?}(xv2qA){|9u<zUxGld$|t7032VO^?$3b<_=+BI
zIOmwBM|1dJ?jXApL~2D*E4_XG6!iaBrTsrF$G->oUozx>@|x!=n;vNwqL(NC=S%kd
zJ0|p(5&SouzwP9o9plj{Ah$@4Arl4U7O8&AExK=02ITI1JAiWYyvqM_Zqb(ir*ez9
zxb{E$k%}Y}{5=)vp4z{pBDr7ro{Hpt1xQ5#N<xxu?gYw<4D%(~^r5*v9pA|-2ikE~
zK;TC7%hQ6*+fJq^t0t{IE?s*p9QO+%=gxF-PJH)~eTnZL%Xd5G%6MW<L=`X}oqToj
zjaPM<@Qjq$3ELUPI9Ey6@IFz?T;1c$bvxbA<8KAh_OHA;vF<@;a^^<@<{!p9+M3DA
zE!A}2LhD$Y#32bwW)3H{yxTV3kRMfEL2TTP5qW;7;j(`-F&xofNU`t<TA1m0217S`
zER=LGN_71Cdz9KArC{VJnuiahZD`IqH@dSOMXR!|dZl>K(v~`A8vBqMQ~YJ<%M8rj
zgSzz7<s<c_Pi{V+Azrjhapz;umdn+2%KBrc*xR%#6K|d<k1XpA*w!*3R`H7qVg9qr
z>RNQU1N2oNt7TGn%KOBv#cHXKitomY5}tGPM+0p06U2ji(S=X6zNW;)`$7|H&-%&F
z-~#%~PhRB=KF3W!T{mV%I?21!k1OMx*Kn0Sd}@UBT?P)zPnPXu%lk9P_NOmF9j~qD
ze26jQtk7BFi}%q1WF1_MTfw2=uQ#MyX_XXWBE>$?$hf_#NcqYApf9MbXL;^}pZg7;
z?3FdmIEG`o$(VXQgwLfpYxPPhKmr$J&rh1cFD+ajG)*7Wa_7XlEsQX>KX+JeC(8;A
zpdTwAcQ?82z9pSechbb$%V)#a=g~Vg^bkju5Ii)d*Q&DLE^rVw8g^->xuD3zDxKkJ
zRzz*|SPb)>(H8D}d*B4sdq_fRt|ud=Dj_LBGdHm}L+3ec!_~&8d5IXl)bQCL;99e$
zapjfy2G@^4R#~YNqXDUxOjf7-SljPS1f^3|e2@-K>g@SoHto7!gEIErGH4RB7!|As
z$52(ouLb#(UO@Sl^qvg1I+#{)-Z$KL(0g=a=59|U?QPn|+AzVe5Yn^J_mEvM96mN)
zDeDx_{UDvN!Iak~9`k7hYGYU8(F~<_DER@6y8Keq!=*CErC~YX>QxfO47u_vPmCz_
zbair-jj_YYP|tqolhh2Kw!K~ls(oIs%ldH4tomiizAQuSQeCDq+|z^ob>VrrVo!t;
zDH@VyqD`%Y5InHQ*prp1n6-vfr_m%5X5y>2%Fx;ahmF-rH%gj+D}T?=Q%N<@N6#+y
z+ejTc(D3rty9JYOcdMN07Q2@9^h3ysWtbN0y18cs+9bAX@sQ>U=jLuLc^hnf*kn;|
zr%^H?Gh}M=x>JS%Q}-wbl7&oDfa83pD_b+u6gHM(6gyX^KUt%{+^K10(>se6zRG5?
z8|KOllJ!~82KSXNYxJz2Y}V?md(lG76aM!N6C~?ZLckjA#KI@*)rkdVvy1?nVjWFJ
zmay4Sx<WfSus>uTlB{4htT}eSFUD$OfX(sxMf8wqa-m1oM&16@!nSC$q2YH`LB%?Z
zWK>5Z^|rMG`sID9(kR>}t8pg1fwmG_ynbml?PblQ`wEoC=OG&}8OzYc^v1O}WV2+^
zyU42(NJ>GQu~huf6y4iu0;gY|7UDuGORa%p)^c>MCX!XJa{_zbTBCzLTcqwJ>e?U=
zJfC316y&wdq!Z9nk*er$R=>^N>Vmf7)NL3#A<HmXJxgmYZgf6Np#-vz9rp{{YGT7z
zJZLShFuyImk58p`y0vKGV>`J^El=8JL@T6Wi|<Eun-R`C|GVKXw#AolVLfOonqgl-
zBQ*}5=U2SEzjE&R`T3OMjt}erna<Af1$R+j`ndHRD{+;7o_Vc!nbm)nIW$GK3aW&5
zT%U@Hv8>+6amO*%Ao<e~>}^E<>9``TuQzKjffEzx@R^Yo2JT*QI$Y~SoX$y`np#iy
z8FA6P>Fk(ooL{8}@}^T_ps=q!Ggd)gGx(6CwN7jCW6K?HWY4`GFBWd!Ix)`QNQ)dV
zq%sDl>$i|_OPr$mr5<hBj9Q1{-K@&Hg{Tel;@zW5v!^i6jrx4nMiz$^^(X1~`5^D7
z^Q_rT%(b2538|v`IZ^!);7|&ypk@NR>f@J${OQ-u%*^ATZ(9$SQ3ynh@l_p|M(9NK
zYHe;AYe+Zwyk*ehhk@NF^y%Vw&Cit|=QRP|A@IPN%Fgz2Mg5wC1d7|N51;7bYEW@}
z!<8O*ysowlJoE)KYyHh?krtL<Q$tYIv<cbwE~G05*D>zKfAA`O#d^46q3BdXscBy(
zETGldZ#J=b?{4f=3pp*RSUUkdP1VQgCQzIM7_I!KeJQx~wH`tJ{OM*s3VF6e$}*^D
z{Q!jXE73YrXX;Kvf%Ah_^-dTgD%CQmZ7w&9^lNTWeS3cU*32ZbN`p%H*I#(htM8`D
z&#h7KS_Re2PN9_H<nXL<+e7qBAVXsdQ`_~6`Xemb@8tD#PZZ=t`1@&g0m}c^qP?Uw
zs1NOzk)pZiZ=r8#1HzH1(qDYHw)KgPd-DI`Q1hwaj+KqtE4QG<(!i28`vAB;3_j)6
z;A@|Acyqm-ExNsyyqF`}Zrz~tl{G>o(>V3Sqr0|WL2$FMTU-xBH#zRmNcG;H@0hUX
zso(W??H7&#E1EAsz)O1r{V8M(Vz128Lp)=|BsY({vCEl0g)qFyDTz+~W|IpY{RR?X
zOQ*nNTF^|e6((1};IwzU!GRVpr5jCaz)ryFX6v1M{n@dDwF7DT<6=g<fn)@lkkad|
z{j%)&`H!4Fsjft@rRu}8ERqFQ-}{SKL;RYJX|^`FCt~oipxbWm21+in`7(6R`!ujC
zuF%!OA33%PZ)>lp5o^{8n_2H@xX$Cck>@Hl5wK@UE@3V1(y#?OH+W(D;#T|`dY^N;
zkC<Kauhsil@V*BQ&99-)W)3e%_48@fe!0TlU2k>HrI9x;XFx3{*a8mq$L0_tE6T?)
z=kI6gVKp-)bA{TMj9VUn!}1Ucmu%p&<KFw_JB{T}I-b#Ik#iT$vVFg-j~t^^2Yq<o
zU-^6niCUVbwK{r9<!#rCrLe}=g=m!NrI_vcl4rlJ-MUq)x}%v@Ty@?j$g@>4WnR&9
zKoX9nIlYJsco4lFneVl9@m;w`twa|ec$(_Y9K1L~mF(0XLP@>~F!l85Hp{qaP-xjG
z2X!Avv_Qr$?NH$j%msFPO#Kgy!H(zKno|;(9(0nDt3TEi*J_zFYLE)#TR1@%8pIC7
z_4aY9q`vy%g;-ceW%U-#s6ZV`Ge)}8q}vd<fSJ(diF^ayW_6yq>(_KWK{Hko%NGV4
zP8#XHVpV9!`fo&0Q#CyD8G8esFc0jP#K$2dguM+|?^Jq_dI#;AL+2A17BJ}jGMjST
zZYvnHNKTY@@e)21p1l0{iKK~wt4?dIVoyHHpeiOnP&JF;{p49`M*Mmtun{})-U3}m
zQo}FF3$RF76wptoZ`3F>p`;^+!{jnzCToaxquy!?i&oIurykqITaYuLCfV9Nma3s_
z&atu-i;aUaS?7ThJ(_#a^7Tbzh+<(o3$Dj=9|7GHiOIOsFXN1ZpDZ)bS7MK7%u9JU
zJ65jM8RiGxFHO~VJkoUlJWV!>XI@aQb{L?T4CzY(8)p8UZQfUsMlBkveeC8vPUrXD
z5qqtn>nO2&YOwNT$J%JaLA$&dCaQukv*g{f1kWXAb$`mq4HHjR|6+R*ze~~_XPUq=
zT4lDjC;|I*%K=p}fw^yjtw$a^+cDG^L!C1VgBt<)+OY%JTPchpdn6m3Z=lbLv?!M=
zgBPAbe~tl%RU)Fv8+*w7-d~9iqV-;JpsiV#1Q;`(jAMimWpf-LZ_BmtXGEH$ct!WL
zwu0X~Qy4x=PJmU%K^sqHG;2nORopQsd@a}V4KCuz1B&$$U0Zm@(6B5s*Jybaqlyd&
z3+`G}e*_=*lZ0?dB#t62(sy}OAKUzxTvn0G(NaJS(5q#;_Ty{I^KvXS2g@_NtwY>|
zBLmDiI?@a3{XyqN@j-%-Ot+#hM>>u3a%ry<stjDL_C&H#IJj;2e9!A_B-!Hj?uA^;
zwmzY~U+#k{aO1=?^%k?Uk~xIV=T-!S8ITunu%aZuX<*q|+)f)beAtB?a<^XfdI=&t
ze$<%}f^0Azsm&CO#E4O#-6t2MhULn98qdy6?EMl48phzJ;x0`Rr+|%{{BE|SPraJ0
zsmm_>S?a~O=wsj9Z&R4iRk&k9ZT)YrjZYeK)5=f2NgdA7mj9FyxaVBZ=GyK{QUvQ<
zVXd{yR}!eRoAvDVk*+78PWPIH0~vq)oyMOnm%GNn>rA<@A!FB{{`oag_GxYkLwrxs
zpS?!^pXvH<%KwMw;lBm`x4{4DCFh?xgug-ZH%R^l$=_D+w-x+t1%F$?-&XLKtN?#V
i{_Dn63i)SXZ-(=21FGpCCD)23|H1|o`2WuvQ2!rOV{6R-
literal 0
HcmV?d00001
diff --git a/mpddata/MpdFillDstTask.cxx b/mpddata/MpdFillDstTask.cxx
index c52653c..fd057ca 100644
--- a/mpddata/MpdFillDstTask.cxx
+++ b/mpddata/MpdFillDstTask.cxx
@@ -34,231 +34,230 @@ using std::cout;
using std::endl;
// ----- constructor with names ------------------------------------
+
MpdFillDstTask::MpdFillDstTask(const char *name, const char *title)
- : FairTask(name)
-{
+: FairTask(name) {
fEvent = NULL;
fZdcSkeletonesSaved = 0;
}
// ----- Destructor -----------------------------------------------
-MpdFillDstTask::~MpdFillDstTask()
-{
- if (fEvent) delete fEvent;
+
+MpdFillDstTask::~MpdFillDstTask() {
+ if (fEvent) delete fEvent;
}
// -------------------------------------------------------------------
-InitStatus MpdFillDstTask::Init()
-{
- fEvent = new MpdEvent();
-
- FairRootManager *manager = FairRootManager::Instance();
-
- fKFTracks = (TClonesArray *) manager->GetObject("TpcKalmanTrack");
- fKFEctTracks = (TClonesArray *) manager->GetObject("EctTrack");
- fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
- fMCEventHeader = (FairMCEventHeader*) manager->GetObject("MCEventHeader.");
- fTofMatching = (TClonesArray *) manager->GetObject("TOFMatching");
- fEtofMatching = (TClonesArray *) manager->GetObject("ETOFMatching");
- fVertex = (TClonesArray *) manager->GetObject("Vertex"); //AZ
-
- fELossZdc1Value = (TClonesArray *) manager->GetObject("ELossZdc1Value"); //EL
- fELossZdc2Value = (TClonesArray *) manager->GetObject("ELossZdc2Value"); //EL
- fELossZdc1Histo = (TClonesArray *) manager->GetObject("ELossZdc1Histo"); //EL
- fELossZdc2Histo = (TClonesArray *) manager->GetObject("ELossZdc2Histo"); //EL
- fHistZdc1En = (TH2F *) manager->GetObject("HistZdc1En"); //EL
- fHistZdc2En = (TH2F *) manager->GetObject("HistZdc2En"); //EL
- fZdcSkeletonesSaved=0;
-
- FairRootManager::Instance()->Register("MCEventHeader.","MC", fMCEventHeader, kTRUE);
- FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
- //FairRootManager::Instance()->Register("MpdEvent","MpdEvents", fEvents, kTRUE);
- FairRootManager::Instance()->Register("MPDEvent.","MpdEvent", fEvent, kTRUE);
- // FairRootManager::Instance()->Register("EZdc1","EZdc1", fELossZdc1Histo, kTRUE);
- // FairRootManager::Instance()->Register("EZdc2","EZdc2", fELossZdc2Histo, kTRUE);
-
- //fhTrackMotherId = new TH1F("TrackMotherId","",1000,-5,995);
- //fhTrackPrimaryPDG = new TH1F("TrackPrimaryPDG","",1000,-500,500);
- //fhTrackVertex = new TH2F("TrackVertex","",1000,-500,500,1000,-500,500);
- //fhTruthVertex = new TH2F("TruthVertex","",1000,-500,500,1000,-500,500);
-
- return kSUCCESS;
+
+InitStatus MpdFillDstTask::Init() {
+ fEvent = new MpdEvent();
+
+ FairRootManager *manager = FairRootManager::Instance();
+
+ fKFTracks = (TClonesArray *) manager->GetObject("TpcKalmanTrack");
+ fKFEctTracks = (TClonesArray *) manager->GetObject("EctTrack");
+ fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
+ fMCEventHeader = (FairMCEventHeader*) manager->GetObject("MCEventHeader.");
+ fTofMatching = (TClonesArray *) manager->GetObject("TOFMatching");
+ fEtofMatching = (TClonesArray *) manager->GetObject("ETOFMatching");
+ fVertex = (TClonesArray *) manager->GetObject("Vertex"); //AZ
+
+ fELossZdc1Value = (TClonesArray *) manager->GetObject("ELossZdc1Value"); //EL
+ fELossZdc2Value = (TClonesArray *) manager->GetObject("ELossZdc2Value"); //EL
+ fELossZdc1Histo = (TClonesArray *) manager->GetObject("ELossZdc1Histo"); //EL
+ fELossZdc2Histo = (TClonesArray *) manager->GetObject("ELossZdc2Histo"); //EL
+ fHistZdc1En = (TH2F *) manager->GetObject("HistZdc1En"); //EL
+ fHistZdc2En = (TH2F *) manager->GetObject("HistZdc2En"); //EL
+ fZdcSkeletonesSaved = 0;
+
+ FairRootManager::Instance()->Register("MCEventHeader.", "MC", fMCEventHeader, kTRUE);
+ FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
+ //FairRootManager::Instance()->Register("MpdEvent","MpdEvents", fEvents, kTRUE);
+ FairRootManager::Instance()->Register("MPDEvent.", "MpdEvent", fEvent, kTRUE);
+ // FairRootManager::Instance()->Register("EZdc1","EZdc1", fELossZdc1Histo, kTRUE);
+ // FairRootManager::Instance()->Register("EZdc2","EZdc2", fELossZdc2Histo, kTRUE);
+
+ //fhTrackMotherId = new TH1F("TrackMotherId","",1000,-5,995);
+ //fhTrackPrimaryPDG = new TH1F("TrackPrimaryPDG","",1000,-500,500);
+ //fhTrackVertex = new TH2F("TrackVertex","",1000,-500,500,1000,-500,500);
+ //fhTruthVertex = new TH2F("TruthVertex","",1000,-500,500,1000,-500,500);
+
+ return kSUCCESS;
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Exec(Option_t * option)
-{
- Reset(); // Clear previos event information
-
- if (!fZdcSkeletonesSaved) { // empty skeletones saved only once
- FairRootManager* ioman = FairRootManager::Instance();
- if (fHistZdc1En) {
- fHistZdc1En->SetDirectory((TFile*)ioman->GetOutFile());
- fHistZdc2En->SetDirectory((TFile*)ioman->GetOutFile());
- fHistZdc1En->Write();
- fHistZdc2En->Write();
- fZdcSkeletonesSaved=1;
+
+void MpdFillDstTask::Exec(Option_t * option) {
+ Reset(); // Clear previos event information
+
+ if (!fZdcSkeletonesSaved) { // empty skeletones saved only once
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (fHistZdc1En) {
+ fHistZdc1En->SetDirectory((TFile*) ioman->GetOutFile());
+ fHistZdc2En->SetDirectory((TFile*) ioman->GetOutFile());
+ fHistZdc1En->Write();
+ fHistZdc2En->Write();
+ fZdcSkeletonesSaved = 1;
+ }
+ }
+
+ Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
+ Int_t nEctReco = fKFEctTracks ? fKFEctTracks->GetEntriesFast() : 0;
+ cout << "\n-I- [MpdFillDstTask::Exec] " << nReco + nEctReco << " reconstruced tracks to write" << endl;
+
+ FairRunAna *fRun = FairRunAna::Instance();
+ fEvent->SetRunInfoRunId(fRun->GetRunId());
+
+ FairField *fPar = fRun->GetField();
+ fEvent->SetRunInfoMagneticFieldZ(fPar->GetType());
+
+ if (fVertex) {
+ MpdVertex *vertex = (MpdVertex*) fVertex->UncheckedAt(0);
+ fEvent->SetPrimaryVerticesX(vertex->GetX());
+ fEvent->SetPrimaryVerticesY(vertex->GetY());
+ fEvent->SetPrimaryVerticesZ(vertex->GetZ());
+ } else {
+ // Vertex info is not available
+ fEvent->SetPrimaryVerticesX(0.);
+ fEvent->SetPrimaryVerticesY(0.);
+ fEvent->SetPrimaryVerticesZ(0.);
+ }
+
+ // check clone track into ECT
+ typedef std::set<Int_t> trackSet;
+ trackSet EctTrackSet;
+ for (Int_t index = 0; index < nEctReco; index++) // cycle by ECT KF tracks
+ {
+ MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fKFEctTracks->UncheckedAt(index);
+ if (track->IsFromTpc()) EctTrackSet.insert(trackSet::value_type(track->GetTpcIndex()));
+ }
+
+ Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
+ MpdTofMatchingData *pMatchingData;
+ bool matchingDataExist;
+
+ //AZ MpdParticleIdentification *identificator = new MpdParticleIdentification();
+ MpdParticleIdentification identificator;
+
+ for (Int_t i = 0; i < nReco; i++) {
+ // check clone track into ECT
+ if (nEctReco > 0) {
+ trackSet::iterator it = EctTrackSet.find(i);
+ if (it != EctTrackSet.end()) continue;
+ }
+
+ MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
+
+ MpdTrack *track = fEvent->AddGlobalTrack();
+ track->SetID(kftrack->GetTrackID());
+ track->SetNofHits(kftrack->GetNofHits());
+ track->SetdEdXTPC(kftrack->GetPartID());
+ Float_t Ppi, Pk, Pe, Pp;
+
+ if (!identificator.GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 0)) { //0 - equal bayesian coefficients
+ track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2));
+ }
+ matchingDataExist = false;
+ for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++) {
+ pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
+ if (pMatchingData->GetKFTrackIndex() == i) {
+ matchingDataExist = true;
+ break;
+ } // first matching
+ }
+
+ if (matchingDataExist) {
+ track->SetTofBeta(pMatchingData->GetBeta());
+ track->SetTofMass2(pMatchingData->GetMass2());
+ track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
+
+ if (!identificator.GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetBeta(), Ppi, Pk, Pp, Pe, 0)) {
+ track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1));
+ }
+ }
+ Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
+ Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
+ Float_t combProbs[4]; //probabilities combined from TOF & TPC
+ identificator.GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
+ Ppi = combProbs[0];
+ Pk = combProbs[1];
+ Pp = combProbs[2];
+ Pe = combProbs[3];
+ track->SetCombPidProb(Pe, Ppi, Pk, Pp);
+
+ if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
+ else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
+
+ track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
+ track->SetPhi(kftrack->GetParam(2)); // Phi
+
+ TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
+ track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
+ track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
+ track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
+
+ track->SetChi2(kftrack->GetChi2());
+
+ Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
+ track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
+ track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
+ track->SetFirstPointZ(kftrack->GetParam(1));
+
+ track->SetLastPointX(0.); // AZ - currently not available
+ track->SetLastPointY(0.); // AZ - currently not available
+ track->SetLastPointZ(0.); // AZ - currently not available
+ }
+
+
+ Int_t nEMatching = fEtofMatching ? fEtofMatching->GetEntries() : 0;
+ MpdEtofMatchingData *pEMatchingData;
+
+ for (Int_t i = 0; i < nEctReco; i++) {
+ MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFEctTracks->UncheckedAt(i);
+
+ MpdTrack *track = fEvent->AddGlobalTrack();
+ track->SetID(kftrack->GetTrackID());
+ track->SetNofHits(kftrack->GetNofHits());
+ matchingDataExist = false;
+ for (Int_t etofIndex = 0; etofIndex < nEMatching; etofIndex++) {
+ pEMatchingData = (MpdEtofMatchingData*) fEtofMatching->UncheckedAt(etofIndex);
+ if (pEMatchingData->GetKFTrackIndex() == i) {
+ matchingDataExist = true;
+ break;
+ } // first matching
+ }
+
+ if (matchingDataExist) {
+ track->SetTofMass2(pEMatchingData->GetMass2());
+ track->SetTofHitIndex(pEMatchingData->GetTofHitIndex());
+ }
+
+ if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
+ else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
+
+ track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
+ track->SetPhi(kftrack->GetParam(2)); // Phi
+
+ TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
+ track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
+ track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
+ track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
+
+ track->SetChi2(kftrack->GetChi2());
+
+ Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
+ track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
+ track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
+ track->SetFirstPointZ(kftrack->GetParam(1));
+
+ track->SetLastPointX(0.); // AZ - currently not available
+ track->SetLastPointY(0.); // AZ - currently not available
+ track->SetLastPointZ(0.); // AZ - currently not available
}
- }
-
- Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
- Int_t nEctReco = fKFEctTracks ? fKFEctTracks->GetEntriesFast() : 0;
- cout<<"\n-I- [MpdFillDstTask::Exec] "<< nReco + nEctReco<< " reconstruced tracks to write" <<endl;
-
- FairRunAna *fRun = FairRunAna::Instance();
- fEvent->SetRunInfoRunId(fRun->GetRunId());
-
- FairField *fPar = fRun->GetField();
- fEvent->SetRunInfoMagneticFieldZ(fPar->GetType());
-
- if (fVertex)
- {
- MpdVertex *vertex = (MpdVertex*) fVertex->UncheckedAt(0);
- fEvent->SetPrimaryVerticesX(vertex->GetX());
- fEvent->SetPrimaryVerticesY(vertex->GetY());
- fEvent->SetPrimaryVerticesZ(vertex->GetZ());
- } else {
- // Vertex info is not available
- fEvent->SetPrimaryVerticesX(0.);
- fEvent->SetPrimaryVerticesY(0.);
- fEvent->SetPrimaryVerticesZ(0.);
- }
-
- // check clone track into ECT
- typedef std::set<Int_t> trackSet;
- trackSet EctTrackSet;
- for(Int_t index = 0; index < nEctReco; index++) // cycle by ECT KF tracks
- {
- MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fKFEctTracks->UncheckedAt(index);
- if(track->IsFromTpc()) EctTrackSet.insert(trackSet::value_type(track->GetTpcIndex()));
- }
-
- Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
- MpdTofMatchingData *pMatchingData;
- bool matchingDataExist;
-
- //AZ MpdParticleIdentification *identificator = new MpdParticleIdentification();
- MpdParticleIdentification identificator;
-
- for(Int_t i = 0; i < nReco; i++)
- {
- // check clone track into ECT
- if(nEctReco>0)
- {
- trackSet::iterator it = EctTrackSet.find(i);
- if(it != EctTrackSet.end()) continue;
- }
-
- MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
-
- MpdTrack *track = fEvent->AddGlobalTrack();
- track->SetID(kftrack->GetTrackID());
- track->SetNofHits(kftrack->GetNofHits());
- track->SetdEdXTPC(kftrack->GetPartID());
- Float_t Ppi, Pk, Pe, Pp;
-
- if(!identificator.GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 1)) {
- track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2));
- }
- matchingDataExist = false;
- for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++)
- {
- pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
- if(pMatchingData->GetKFTrackIndex() == i){ matchingDataExist = true; break;} // first matching
- }
-
- if(matchingDataExist)
- {
- track->SetTofBeta(pMatchingData->GetBeta());
- track->SetTofMass2(pMatchingData->GetMass2());
- track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
-
- if(!identificator.GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetBeta(), Ppi, Pk, Pp, Pe, 0)) {
- track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1));
- }
- }
- Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
- Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
- Float_t combProbs[4]; //probabilities combined from TOF & TPC
- identificator.GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
- Ppi = combProbs[0];
- Pk = combProbs[1];
- Pp = combProbs[2];
- Pe = combProbs[3];
- track->SetCombPidProb(Pe, Ppi, Pk, Pp);
-
- if ( kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
- else track->SetPt(1./kftrack->GetParam(4)); /*signed Pt*/
-
- track->SetTheta(TMath::PiOver2()-kftrack->GetParam(3)); // Theta: angle from beam line
- track->SetPhi(kftrack->GetParam(2)); // Phi
-
- TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
- track->SetPtError(TMath::Sqrt(Cov(4,4))); // Pt error
- track->SetThetaError(TMath::Sqrt(Cov(3,3))); // Theta error
- track->SetPhiError(TMath::Sqrt(Cov(2,2))); // Phi error
-
- track->SetChi2(kftrack->GetChi2());
-
- Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
- track->SetFirstPointX(kftrack->GetPosNew()*TMath::Cos(phi)); // closest to beam line
- track->SetFirstPointY(kftrack->GetPosNew()*TMath::Sin(phi));
- track->SetFirstPointZ(kftrack->GetParam(1));
-
- track->SetLastPointX(0.); // AZ - currently not available
- track->SetLastPointY(0.); // AZ - currently not available
- track->SetLastPointZ(0.); // AZ - currently not available
- }
-
-
- Int_t nEMatching = fEtofMatching ? fEtofMatching->GetEntries() : 0;
- MpdEtofMatchingData *pEMatchingData;
-
- for (Int_t i = 0; i < nEctReco; i++)
- {
- MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFEctTracks->UncheckedAt(i);
-
- MpdTrack *track = fEvent->AddGlobalTrack();
- track->SetID(kftrack->GetTrackID());
- track->SetNofHits(kftrack->GetNofHits());
- matchingDataExist = false;
- for (Int_t etofIndex = 0; etofIndex < nEMatching; etofIndex++)
- {
- pEMatchingData = (MpdEtofMatchingData*) fEtofMatching->UncheckedAt(etofIndex);
- if(pEMatchingData->GetKFTrackIndex() == i){ matchingDataExist = true; break;} // first matching
- }
-
- if(matchingDataExist)
- {
- track->SetTofMass2(pEMatchingData->GetMass2());
- track->SetTofHitIndex(pEMatchingData->GetTofHitIndex());
- }
-
- if ( kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
- else track->SetPt(1./kftrack->GetParam(4)); /*signed Pt*/
-
- track->SetTheta(TMath::PiOver2()-kftrack->GetParam(3)); // Theta: angle from beam line
- track->SetPhi(kftrack->GetParam(2)); // Phi
-
- TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
- track->SetPtError(TMath::Sqrt(Cov(4,4))); // Pt error
- track->SetThetaError(TMath::Sqrt(Cov(3,3))); // Theta error
- track->SetPhiError(TMath::Sqrt(Cov(2,2))); // Phi error
-
- track->SetChi2(kftrack->GetChi2());
-
- Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
- track->SetFirstPointX(kftrack->GetPosNew()*TMath::Cos(phi)); // closest to beam line
- track->SetFirstPointY(kftrack->GetPosNew()*TMath::Sin(phi));
- track->SetFirstPointZ(kftrack->GetParam(1));
-
- track->SetLastPointX(0.); // AZ - currently not available
- track->SetLastPointY(0.); // AZ - currently not available
- track->SetLastPointZ(0.); // AZ - currently not available
- }
}
// -------------------------------------------------------------------
//Delete MC tracks being outside the MPD
-void MpdFillDstTask::CleanMC(){
- cout<<"-I- [MpdFillDstTask::Exec] Cleaning from outer decays..."<<flush;
+
+void MpdFillDstTask::CleanMC() {
+ cout << "-I- [MpdFillDstTask::Exec] Cleaning from outer decays..." << flush;
Int_t nMC = fMCTracks->GetEntriesFast(), motherId;
FairMCTrack *track;
@@ -272,24 +271,21 @@ void MpdFillDstTask::CleanMC(){
vector<int> removedMother;
Int_t i, j;
- for (i = 0; i < nMC; i++)
- {
+ for (i = 0; i < nMC; i++) {
track = (FairMCTrack*) fMCTracks->UncheckedAt(i);
//filling hostograms
motherId = track->GetMotherId();
//fhTrackMotherId->Fill(motherId);
- if (motherId == -1){
+ if (motherId == -1) {
//fhTrackPrimaryPDG->Fill(track->GetPdgCode());
- }
- else{
+ } else {
//fhTrackVertex->Fill(track->GetStartX(),track->GetStartY());
//if decay was out of MPD
- if (sqrt(track->GetStartX()*track->GetStartX()+track->GetStartY()*track->GetStartY()) > rMax){
+ if (sqrt(track->GetStartX() * track->GetStartX() + track->GetStartY() * track->GetStartY()) > rMax) {
fMCTracks->Remove(track);
removedMother.push_back(i);
- }
- else{
+ } else {
//motherId < childId (i) REQUIRED
/*for (j = 0; j < cntRemovedMother; j++){
if (motherId == removedMother[j]){
@@ -320,7 +316,7 @@ void MpdFillDstTask::CleanMC(){
cntRemovedMother++;
}*/
//binary search
- if (binary_search (removedMother.begin(), removedMother.end(), motherId)){
+ if (binary_search(removedMother.begin(), removedMother.end(), motherId)) {
fMCTracks->Remove(track);
removedMother.push_back(i);
}
@@ -331,35 +327,36 @@ void MpdFillDstTask::CleanMC(){
}
fMCTracks->Compress();
- cout<<endl;
+ cout << endl;
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Reset()
-{
+
+void MpdFillDstTask::Reset() {
fEvent->Reset();
}
// -------------------------------------------------------------------
-void MpdFillDstTask::Finish()
-{
- //cout<<"\n-I- [MpdFillDstTask::Finish] "<< endl;
- //fEvents->Dump();
- //cout << "\n";
-
- /* !this code gives some problems! -- !add extra event and so segfault!
- FairRootManager *fManager = FairRootManager::Instance();
- fManager->Fill();
- */
- //fhTrackMotherId->Write();
- //fhTrackPrimaryPDG->Write();
- //fhTrackVertex->Write();
- //fhTruthVertex->Write();
-
- delete fEvent;
- fEvent = NULL;
+
+void MpdFillDstTask::Finish() {
+ //cout<<"\n-I- [MpdFillDstTask::Finish] "<< endl;
+ //fEvents->Dump();
+ //cout << "\n";
+
+ /* !this code gives some problems! -- !add extra event and so segfault!
+ FairRootManager *fManager = FairRootManager::Instance();
+ fManager->Fill();
+ */
+ //fhTrackMotherId->Write();
+ //fhTrackPrimaryPDG->Write();
+ //fhTrackVertex->Write();
+ //fhTruthVertex->Write();
+
+ delete fEvent;
+ fEvent = NULL;
}
// -------------------------------------------------------------------
+
/*MpdEvent *MpdFillDstTask::AddEvent(Option_t * option)
{
TClonesArray &events = *fEvents;
@@ -368,7 +365,7 @@ void MpdFillDstTask::Finish()
return event;
}*/
-MpdTrack *MpdFillDstTask::AddPrimaryTrack(){
+MpdTrack *MpdFillDstTask::AddPrimaryTrack() {
return NULL;
}
// -------------------------------------------------------------------
diff --git a/tpc/CMakeLists.txt b/tpc/CMakeLists.txt
index bfd935b..63dbdc2 100644
--- a/tpc/CMakeLists.txt
+++ b/tpc/CMakeLists.txt
@@ -50,6 +50,7 @@ MpdTpcHitProducer.cxx
MpdTpcSectorGeo.cxx
MpdTpcDigitizerAZ.cxx
MpdTpcClusterFinderAZ.cxx
+MpdTPCpid.cxx
)
# fill list of header files from list of source files
diff --git a/tpc/MpdParticleIdentification.cxx b/tpc/MpdParticleIdentification.cxx
index 34439f1..554a64a 100644
--- a/tpc/MpdParticleIdentification.cxx
+++ b/tpc/MpdParticleIdentification.cxx
@@ -32,361 +32,6 @@ MpdParticleIdentification::MpdParticleIdentification() {
MpdParticleIdentification::~MpdParticleIdentification() {
}
-Int_t MpdParticleIdentification::GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
-
- /*Parameters for fitting*/
- //AZ Float_t *ProtonPar, *PionPar, *KaonPar, *ElectronPar;
- /************************/
-
- dedx *= 1000000; // converting to keV/cm
- const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
- Float_t ProtonPar[9], PionPar[9], KaonPar[9], ElectronPar[9], resultProb[Ntypes]; //AZ
- const Int_t Nintervals = 10;
- Float_t sigmas[Ntypes][Nintervals] = {//array of sigmas for different particles and different intervals of momentum. They were got from 100000 BOX events
- {0.113, 0.107, 0.100, 0.107, 0.097, 0.111, 0.122, 0.105, 0.115, 0.118},
- {0.407, 0.276, 0.260, 0.159, 0.185, 0.154, 0.122, 0.105, 0.115, 0.118},
- {0.507, 0.488, 0.394, 0.330, 0.268, 0.244, 0.260, 0.172, 0.193, 0.118},
- {0.163, 0.160, 0.149, 0.159, 0.172, 0.205, 0.260, 0.172, 0.193, 0.156}
- };
-
- Float_t sigma = 1.0 / TMath::Sqrt(nHits); //for gaussian
- Float_t sum = 0.0; // for normalizing
-
- /*A priori coefficients for Bayes approach */
- Float_t bayesAprioriCoefficients[Ntypes];
- /*A "measured" probabilities */
- Float_t gausProb[Ntypes];
-
- Float_t sigPi, sigPr, sigKa, sigEl; //sigmas for gaussians
-
- switch (method) {
-
- case 0: //bethe-bloch approximation with "standard" sigmas and equal bayesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- sigma = 0.07 * dedx;
- sigPi = sigPr = sigKa = sigEl = sigma;
- bayesAprioriCoefficients[0] = 1.0;
- bayesAprioriCoefficients[1] = 1.0;
- bayesAprioriCoefficients[2] = 1.0;
- bayesAprioriCoefficients[3] = 1.0;
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 1: //bethe-bloch approximation with special different sigmas and byesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- if (P < 0.3) {
- bayesAprioriCoefficients[0] = 0.28;
- bayesAprioriCoefficients[1] = 0.25;
- bayesAprioriCoefficients[2] = 0.26;
- bayesAprioriCoefficients[3] = 0.21;
- sigPi = sigmas[0][0];
- sigKa = sigmas[1][0];
- sigPr = sigmas[2][0];
- sigEl = sigmas[3][0];
- } else if ((P >= 0.3) && (P < 0.4)) {
- bayesAprioriCoefficients[0] = 0.91;
- bayesAprioriCoefficients[1] = 0.02;
- bayesAprioriCoefficients[2] = 0.06;
- bayesAprioriCoefficients[3] = 0.01;
- sigPi = sigmas[0][1];
- sigKa = sigmas[1][1];
- sigPr = sigmas[2][1];
- sigEl = sigmas[3][1];
- } else if ((P >= 0.4) && (P < 0.5)) {
- bayesAprioriCoefficients[0] = 0.70;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.13;
- bayesAprioriCoefficients[3] = 0.06;
- sigPi = sigmas[0][2];
- sigKa = sigmas[1][2];
- sigPr = sigmas[2][2];
- sigEl = sigmas[3][2];
- } else if ((P >= 0.5) && (P < 0.6)) {
- bayesAprioriCoefficients[0] = 0.39;
- bayesAprioriCoefficients[1] = 0.19;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.10;
- sigPi = sigmas[0][3];
- sigKa = sigmas[1][3];
- sigPr = sigmas[2][3];
- sigEl = sigmas[3][3];
- } else if ((P >= 0.6) && (P < 0.7)) {
- bayesAprioriCoefficients[0] = 0.41;
- bayesAprioriCoefficients[1] = 0.17;
- bayesAprioriCoefficients[2] = 0.37;
- bayesAprioriCoefficients[3] = 0.5;
- sigPi = sigmas[0][4];
- sigKa = sigmas[1][4];
- sigPr = sigmas[2][4];
- sigEl = sigmas[3][4];
- } else if ((P >= 0.7) && (P < 0.8)) {
- bayesAprioriCoefficients[0] = 0.43;
- bayesAprioriCoefficients[1] = 0.16;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.10;
- sigPi = sigmas[0][5];
- sigKa = sigmas[1][5];
- sigPr = sigmas[2][5];
- sigEl = sigmas[3][5];
- } else if ((P >= 0.8) && (P < 0.9)) {
- bayesAprioriCoefficients[0] = 0.44;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.43;
- bayesAprioriCoefficients[3] = 0.02;
- sigPi = sigmas[0][6];
- sigKa = sigmas[1][6];
- sigPr = sigmas[2][6];
- sigEl = sigmas[3][6];
- } else if ((P >= 0.9) && (P < 1.0)) {
- bayesAprioriCoefficients[0] = 0.36;
- bayesAprioriCoefficients[1] = 0.21;
- bayesAprioriCoefficients[2] = 0.36;
- bayesAprioriCoefficients[3] = 0.07;
- sigPi = sigmas[0][7];
- sigKa = sigmas[1][7];
- sigPr = sigmas[2][7];
- sigEl = sigmas[3][7];
- } else if ((P >= 1.0) && (P < 1.2)) {
- bayesAprioriCoefficients[0] = 0.32;
- bayesAprioriCoefficients[1] = 0.32;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.04;
- sigPi = sigmas[0][8];
- sigKa = sigmas[1][8];
- sigPr = sigmas[2][8];
- sigEl = sigmas[3][8];
- } else if (P >= 1.2) {
- bayesAprioriCoefficients[0] = 0.34;
- bayesAprioriCoefficients[1] = 0.27;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.08;
- sigPi = sigmas[0][9];
- sigKa = sigmas[1][9];
- sigPr = sigmas[2][9];
- sigEl = sigmas[3][9];
- }
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 2: //bethe-bloch approximation with "standard" sigmas and different byesian coefficients
-
- //parameters were got from Bethe-Bloch approximation for 100000 BOX events
- /*AZ
- ProtonPar = new Float_t[4];
- PionPar = new Float_t[4];
- KaonPar = new Float_t[4];
- ElectronPar = new Float_t[4];
- */
- ProtonPar[0] = -3.957;
- ProtonPar[1] = 3.525;
- ProtonPar[2] = 16.468;
- ProtonPar[3] = 0.815;
- ProtonPar[4] = 5.247;
- PionPar[0] = -0.739;
- PionPar[1] = 7.373;
- PionPar[2] = 3904.790;
- PionPar[3] = 0.274;
- PionPar[4] = 5.497;
- KaonPar[0] = -2.590;
- KaonPar[1] = 4.918;
- KaonPar[2] = 79.722;
- KaonPar[3] = 0.357;
- KaonPar[4] = 4.511;
- ElectronPar[0] = -1.552;
- ElectronPar[1] = 1.748;
- ElectronPar[2] = 7.425;
- ElectronPar[3] = 0.980;
- ElectronPar[4] = 1.604;
-
- sigma = 0.07 * dedx;
- sigPi = sigPr = sigKa = sigEl = sigma;
-
- if (P < 0.3) {
- bayesAprioriCoefficients[0] = 0.28;
- bayesAprioriCoefficients[1] = 0.25;
- bayesAprioriCoefficients[2] = 0.26;
- bayesAprioriCoefficients[3] = 0.21;
- } else if ((P >= 0.3) && (P < 0.4)) {
- bayesAprioriCoefficients[0] = 0.91;
- bayesAprioriCoefficients[1] = 0.02;
- bayesAprioriCoefficients[2] = 0.06;
- bayesAprioriCoefficients[3] = 0.01;
- } else if ((P >= 0.4) && (P < 0.5)) {
- bayesAprioriCoefficients[0] = 0.70;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.13;
- bayesAprioriCoefficients[3] = 0.06;
- } else if ((P >= 0.5) && (P < 0.6)) {
- bayesAprioriCoefficients[0] = 0.39;
- bayesAprioriCoefficients[1] = 0.19;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.10;
- } else if ((P >= 0.6) && (P < 0.7)) {
- bayesAprioriCoefficients[0] = 0.41;
- bayesAprioriCoefficients[1] = 0.17;
- bayesAprioriCoefficients[2] = 0.37;
- bayesAprioriCoefficients[3] = 0.5;
- } else if ((P >= 0.7) && (P < 0.8)) {
- bayesAprioriCoefficients[0] = 0.43;
- bayesAprioriCoefficients[1] = 0.16;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.10;
- } else if ((P >= 0.8) && (P < 0.9)) {
- bayesAprioriCoefficients[0] = 0.44;
- bayesAprioriCoefficients[1] = 0.11;
- bayesAprioriCoefficients[2] = 0.43;
- bayesAprioriCoefficients[3] = 0.02;
- } else if ((P >= 0.9) && (P < 1.0)) {
- bayesAprioriCoefficients[0] = 0.36;
- bayesAprioriCoefficients[1] = 0.21;
- bayesAprioriCoefficients[2] = 0.36;
- bayesAprioriCoefficients[3] = 0.07;
- } else if ((P >= 1.0) && (P < 1.2)) {
- bayesAprioriCoefficients[0] = 0.32;
- bayesAprioriCoefficients[1] = 0.32;
- bayesAprioriCoefficients[2] = 0.32;
- bayesAprioriCoefficients[3] = 0.04;
- } else if (P >= 1.2) {
- bayesAprioriCoefficients[0] = 0.34;
- bayesAprioriCoefficients[1] = 0.27;
- bayesAprioriCoefficients[2] = 0.31;
- bayesAprioriCoefficients[3] = 0.08;
- }
-
- gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
- break;
-
- case 3: //parabolic approximation
-
- //parameters were got from parabolic approximation function for 2000 UrQMD events
- ProtonPar[0] = 0.031;
- ProtonPar[1] = -0.124;
- ProtonPar[2] = 1.412;
- PionPar[0] = 0.230;
- PionPar[1] = 0.088;
- PionPar[2] = 1.072;
- KaonPar[0] = 0.676;
- KaonPar[1] = 0.621;
- KaonPar[2] = 0.831;
- ElectronPar[0] = 0.000;
- ElectronPar[1] = -0.018;
- ElectronPar[2] = 2.055;
-
- Float_t invP = 1.0 / P;
-
- gausProb[0] = Gaus(dedx, ParabolicFunction(invP, PionPar), sigPi, kTRUE);
- gausProb[1] = Gaus(dedx, ParabolicFunction(invP, KaonPar), sigKa, kTRUE);
- gausProb[2] = Gaus(dedx, ParabolicFunction(invP, ProtonPar), sigPr, kTRUE);
- gausProb[3] = Gaus(dedx, ParabolicFunction(invP, ElectronPar), sigEl, kTRUE);
-
- sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
- if (sum == 0.) return 1;
- gausProb[0] /= sum;
- gausProb[1] /= sum;
- gausProb[2] /= sum;
- gausProb[3] /= sum;
-
- break;
- }
- //AZ Float_t *resultProb = new Float_t[Ntypes];
- BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
-
- Ppi = resultProb[0];
- Pk = resultProb[1];
- Pp = resultProb[2];
- Pe = resultProb[3];
-
- return 0;
-}
-
Int_t MpdParticleIdentification::GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
const Float_t m2_proton = 0.938 * 0.938; //square of proton mass (GeV)
diff --git a/tpc/MpdParticleIdentification.h b/tpc/MpdParticleIdentification.h
index abdb605..c3d2fea 100644
--- a/tpc/MpdParticleIdentification.h
+++ b/tpc/MpdParticleIdentification.h
@@ -19,15 +19,18 @@
// Base Class Headers ----------------
#include "TSystem.h"
+#include "MpdTPCpid.h"
-class MpdParticleIdentification {
+
+class MpdParticleIdentification : public MpdTPCpid {
public:
// Constructors/Destructors ---------
MpdParticleIdentification();
- ~MpdParticleIdentification();
+ virtual ~MpdParticleIdentification();
- Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ // Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ // The function is realized now in the MpdTPCpid-class with corrected probability coefficients
Int_t GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
Int_t GetCombinedProbs(Float_t *tofProbs, Float_t *tpcProbs, Float_t *resultProbs, Int_t N);
diff --git a/tpc/MpdTPCpid.cxx b/tpc/MpdTPCpid.cxx
new file mode 100644
index 0000000..f7005ab
--- /dev/null
+++ b/tpc/MpdTPCpid.cxx
@@ -0,0 +1,285 @@
+//-----------------------------------------------------------
+// File and Version Information:
+// $Id$
+//
+// Description:
+// Implementation of class ParticleIdentification
+// see ParticleIdentification.h for details
+//
+// Environment:
+// Software developed for MPD at NICA.
+//
+// Author List:
+// Gyulnara Eyyubova
+//
+//-----------------------------------------------------------
+
+// C/C++ Headers ----------------------
+#include <iostream>
+#include <TMath.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TH2.h>
+// This Class' Header ------------------
+#include "MpdTPCpid.h"
+
+using namespace std;
+using namespace TMath;
+
+// Class Member definitions -----------
+
+MpdTPCpid::MpdTPCpid() :
+ProtonPar(),
+PionPar(),
+KaonPar(),
+ElectronPar(),
+sigmasPi(),
+sigmasPr(),
+sigmasKa(),
+sigmasEl(),
+fCoefficients(NULL) {
+ ReadTPCResponse();
+}
+
+MpdTPCpid::~MpdTPCpid() {
+ fCoefficients->Delete();
+}
+
+Int_t MpdTPCpid::GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
+
+
+ const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
+ Float_t resultProb[Ntypes];
+ const Int_t Nintervals = 10;
+
+ Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
+
+ //default sigmas
+ Float_t sigma = 0.07 * dedx;
+
+ Float_t sum = 0.0; // for normalizing
+
+ /*A priori coefficients for Bayes approach */
+ Float_t bayesAprioriCoefficients[Ntypes];
+
+ /*A "measured" probabilities */
+ Float_t gausProb[Ntypes];
+ for (Int_t i = 0; i < Ntypes; i++) gausProb[i] = 0;
+
+ Int_t i_p = 0; // momentum interval
+
+ switch (method) {
+
+ case 0: //bethe-bloch approximation equal bayesian coefficients
+
+ if (P > 3) i_p = 9;
+ for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
+
+
+ bayesAprioriCoefficients[0] = 1.0;
+ bayesAprioriCoefficients[1] = 1.0;
+ bayesAprioriCoefficients[2] = 1.0;
+ bayesAprioriCoefficients[3] = 1.0;
+
+ gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE); // kTRUE = normilized by sqrt(2*Pi)*sigma.
+ gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
+ gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
+ gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
+
+ sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+ if (sum == 0.) return 1;
+ gausProb[0] /= sum;
+ gausProb[1] /= sum;
+ gausProb[2] /= sum;
+ gausProb[3] /= sum;
+ break;
+
+ case 1: //bethe-bloch approximation with special byesian coefficients
+
+ if (P > 3) i_p = 9;
+ for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
+
+ if (P < 0.3) {
+ bayesAprioriCoefficients[0] = 0.28;
+ bayesAprioriCoefficients[1] = 0.25;
+ bayesAprioriCoefficients[2] = 0.26;
+ bayesAprioriCoefficients[3] = 0.21;
+ } else if ((P >= 0.3) && (P < 0.4)) {
+ bayesAprioriCoefficients[0] = 0.91;
+ bayesAprioriCoefficients[1] = 0.02;
+ bayesAprioriCoefficients[2] = 0.06;
+ bayesAprioriCoefficients[3] = 0.01;
+ } else if ((P >= 0.4) && (P < 0.5)) {
+ bayesAprioriCoefficients[0] = 0.70;
+ bayesAprioriCoefficients[1] = 0.11;
+ bayesAprioriCoefficients[2] = 0.13;
+ bayesAprioriCoefficients[3] = 0.06;
+ } else if ((P >= 0.5) && (P < 0.6)) {
+ bayesAprioriCoefficients[0] = 0.39;
+ bayesAprioriCoefficients[1] = 0.19;
+ bayesAprioriCoefficients[2] = 0.32;
+ bayesAprioriCoefficients[3] = 0.10;
+ } else if ((P >= 0.6) && (P < 0.7)) {
+ bayesAprioriCoefficients[0] = 0.41;
+ bayesAprioriCoefficients[1] = 0.17;
+ bayesAprioriCoefficients[2] = 0.37;
+ bayesAprioriCoefficients[3] = 0.5;
+ } else if ((P >= 0.7) && (P < 0.8)) {
+ bayesAprioriCoefficients[0] = 0.43;
+ bayesAprioriCoefficients[1] = 0.16;
+ bayesAprioriCoefficients[2] = 0.31;
+ bayesAprioriCoefficients[3] = 0.10;
+ } else if ((P >= 0.8) && (P < 0.9)) {
+ bayesAprioriCoefficients[0] = 0.44;
+ bayesAprioriCoefficients[1] = 0.11;
+ bayesAprioriCoefficients[2] = 0.43;
+ bayesAprioriCoefficients[3] = 0.02;
+ } else if ((P >= 0.9) && (P < 1.0)) {
+ bayesAprioriCoefficients[0] = 0.36;
+ bayesAprioriCoefficients[1] = 0.21;
+ bayesAprioriCoefficients[2] = 0.36;
+ bayesAprioriCoefficients[3] = 0.07;
+ } else if ((P >= 1.0) && (P < 1.2)) {
+ bayesAprioriCoefficients[0] = 0.32;
+ bayesAprioriCoefficients[1] = 0.32;
+ bayesAprioriCoefficients[2] = 0.32;
+ bayesAprioriCoefficients[3] = 0.04;
+ } else if (P >= 1.2) {
+ bayesAprioriCoefficients[0] = 0.34;
+ bayesAprioriCoefficients[1] = 0.27;
+ bayesAprioriCoefficients[2] = 0.31;
+ bayesAprioriCoefficients[3] = 0.08;
+ }
+
+ gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE);
+ gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
+ gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
+ gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
+
+ sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+ if (sum == 0.) return 1;
+ gausProb[0] /= sum;
+ gausProb[1] /= sum;
+ gausProb[2] /= sum;
+ gausProb[3] /= sum;
+ break;
+
+
+
+ }
+
+
+ BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
+
+ Ppi = resultProb[0];
+ Pk = resultProb[1];
+ Pp = resultProb[2];
+ Pe = resultProb[3];
+
+ return 0;
+}
+
+Float_t MpdTPCpid::BetheBlochFunction(Float_t x, Float_t *p) {
+ Float_t b = 1 / (x / Sqrt(1 + x * x));
+ return p[0] / Power(b, p[3]) * (p[1] - Power(b, p[3]) - Log(p[2] + Power(1 / x, p[4])));
+}
+
+/*Bayes function for probabilities calculation*/
+Int_t MpdTPCpid::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
+
+ /*measProb - measured probabilities from TPC/TOF/etc*/
+ /*aprioriProb - a priori probabilities from some magic*/
+ /*bayesProb - result bayes probabilities */
+ /*N - number of types {pion, kaon, proton, electron} */
+
+ Float_t sumBayes = 0.0;
+ for (Int_t i = 0; i < N; ++i) {
+ sumBayes += measProb[i] * aprioriProb[i];
+ }
+
+ if (sumBayes == 0.) return 1;
+
+ for (Int_t i = 0; i < N; ++i) {
+ bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
+ }
+
+ return 0;
+}
+
+void MpdTPCpid::ReadTPCResponse() {
+ TString c;
+ /* Accessing the DeDx response, obtained with VHLEE generator */
+ fCoefficients = new TFile("$VMCWORKDIR/input/MPDTPCPidResponseVhelle.root");
+
+ if (fCoefficients->IsZombie()) {
+ printf("File MPDTPCPidResponseVhelle.root does not exist.\n");
+ }
+
+ TH2D *hdedx_pi = (TH2D*) fCoefficients->Get("dedx_pi");
+ TH2D *hdedx_pr = (TH2D*) fCoefficients->Get("dedx_pr");
+ TH2D *hdedx_ka = (TH2D*) fCoefficients->Get("dedx_ka");
+ TH2D *hdedx_el = (TH2D*) fCoefficients->Get("dedx_el");
+
+ TF1 *fBetheBlochPion = hdedx_pi->GetFunction("BetheBlochALEPH");
+ PionPar[0] = fBetheBlochPion->GetParameter(0);
+ PionPar[1] = fBetheBlochPion->GetParameter(1);
+ PionPar[2] = fBetheBlochPion->GetParameter(2);
+ PionPar[3] = fBetheBlochPion->GetParameter(3);
+ PionPar[4] = fBetheBlochPion->GetParameter(4);
+
+ TF1 *fBetheBlochProton = hdedx_pr->GetFunction("BetheBlochALEPH");
+ ProtonPar[0] = fBetheBlochProton->GetParameter(0);
+ ProtonPar[1] = fBetheBlochProton->GetParameter(1);
+ ProtonPar[2] = fBetheBlochProton->GetParameter(2);
+ ProtonPar[3] = fBetheBlochProton->GetParameter(3);
+ ProtonPar[4] = fBetheBlochProton->GetParameter(4);
+
+ TF1 *fBetheBlochKaon = hdedx_ka->GetFunction("BetheBlochALEPH");
+ KaonPar[0] = fBetheBlochKaon->GetParameter(0);
+ KaonPar[1] = fBetheBlochKaon->GetParameter(1);
+ KaonPar[2] = fBetheBlochKaon->GetParameter(2);
+ KaonPar[3] = fBetheBlochKaon->GetParameter(3);
+ KaonPar[4] = fBetheBlochKaon->GetParameter(4);
+
+ TF1 *fBetheBlochElectron = hdedx_el->GetFunction("BetheBlochALEPH");
+ ElectronPar[0] = fBetheBlochElectron->GetParameter(0);
+ ElectronPar[1] = fBetheBlochElectron->GetParameter(1);
+ ElectronPar[2] = fBetheBlochElectron->GetParameter(2);
+ ElectronPar[3] = fBetheBlochElectron->GetParameter(3);
+ ElectronPar[4] = fBetheBlochElectron->GetParameter(4);
+
+
+ const Int_t Nintervals = 10;
+ Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
+
+ for (Int_t k = 0; k < Nintervals; k++) {
+
+ Int_t p_low_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k]);
+ Int_t p_up_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k + 1]);
+
+ c = Form("hpi%i", k);
+ TH1D *hpi = (TH1D*) hdedx_pi->ProjectionY(c, p_low_bin, p_up_bin);
+ hpi->Fit("gaus");
+ sigmasPi[k] = hpi->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hpr%i", k);
+ TH1D *hpr = (TH1D*) hdedx_pr->ProjectionY("hpr", p_low_bin, p_up_bin);
+ hpr->Fit("gaus");
+ sigmasPr[k] = hpr->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hka%i", k);
+ TH1D *hka = (TH1D*) hdedx_ka->ProjectionY("hka", p_low_bin, p_up_bin);
+ hka->Fit("gaus");
+ sigmasKa[k] = hka->GetFunction("gaus")->GetParameter(2);
+
+ c = Form("hel%i", k);
+ TH1D *hel = (TH1D*) hdedx_el->ProjectionY("hel", p_low_bin, p_up_bin);
+ hel->Fit("gaus");
+ sigmasEl[k] = hel->GetFunction("gaus")->GetParameter(2);
+ }
+
+ // f->Close();
+
+}
+
+ClassImp(MpdTPCpid);
diff --git a/tpc/MpdTPCpid.h b/tpc/MpdTPCpid.h
new file mode 100644
index 0000000..db0d3fd
--- /dev/null
+++ b/tpc/MpdTPCpid.h
@@ -0,0 +1,56 @@
+//-----------------------------------------------------------
+// File and Version Information:
+// $Id$
+//
+// Description:
+// class for particle identification
+//
+//
+// Environment:
+// Software developed for MPD at NICA.
+//
+// Author List:
+// Gyulnara Eyyubova
+//
+//-----------------------------------------------------------
+
+#ifndef MpdTPCpid_HH
+#define MpdTPCpid_HH
+
+// Base Class Headers ----------------
+#include "TSystem.h"
+#include <TGraph.h>
+#include <TFile.h>
+
+class MpdTPCpid {
+public:
+
+ // Constructors/Destructors ---------
+ MpdTPCpid();
+ virtual ~MpdTPCpid();
+
+ Int_t GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
+ Float_t BetheBlochFunction(Float_t x, Float_t *p);
+ Int_t BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N);
+ void ReadTPCResponse();
+
+private:
+
+ Float_t ProtonPar[5];
+ Float_t PionPar[5];
+ Float_t KaonPar[5];
+ Float_t ElectronPar[5];
+ //const Int_t Nintervals = 10;
+ Float_t sigmasPi[10];
+ Float_t sigmasPr[10];
+ Float_t sigmasKa[10];
+ Float_t sigmasEl[10];
+
+ TFile* fCoefficients;
+
+public:
+ ClassDef(MpdTPCpid, 1)
+
+};
+
+#endif
diff --git a/tpc/tpcLinkDef.h b/tpc/tpcLinkDef.h
index badcf27..9c71ad5 100644
--- a/tpc/tpcLinkDef.h
+++ b/tpc/tpcLinkDef.h
@@ -27,6 +27,7 @@
#pragma link C++ class MpdTpcSectorGeo+;
#pragma link C++ class MpdTpcDigitizerAZ+;
#pragma link C++ class MpdTpcClusterFinderAZ+;
+#pragma link C++ class MpdTPCpid+;
#endif
--
GitLab