From c61683ef13909c651d3981d795e00a5a9bca53a6 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Wed, 30 Mar 2022 11:17:02 +0100
Subject: [PATCH] Add scripts to generator events from Foresee and a sampler to
read these in to the DIF + config for this. Add ability to read HepMC input
events. Add Validation alg
---
.../python/FaserParticleGunConfig.py | 47 ++
Generators/ForeseeGenerator/CMakeLists.txt | 12 +
.../events_14TeV_m0.1GeV_c1e-05to_11_-11.npy | Bin 0 -> 59240 bytes
.../ForeseeGenerator/python/ForeseeSampler.py | 433 ++++++++++++++++++
.../ForeseeGenerator/python/Validate.py | 198 ++++++++
.../share/generate_forsee_events.py | 256 +++++++++++
Generators/HEPMCReader/CMakeLists.txt | 9 +
.../HEPMCReader/python/HepMCReaderConfig.py | 24 +
.../test/G4FaserAlgConfigNew_Test.py | 18 +-
9 files changed, 994 insertions(+), 3 deletions(-)
create mode 100644 Generators/ForeseeGenerator/CMakeLists.txt
create mode 100644 Generators/ForeseeGenerator/data/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy
create mode 100644 Generators/ForeseeGenerator/python/ForeseeSampler.py
create mode 100644 Generators/ForeseeGenerator/python/Validate.py
create mode 100644 Generators/ForeseeGenerator/share/generate_forsee_events.py
create mode 100644 Generators/HEPMCReader/CMakeLists.txt
create mode 100644 Generators/HEPMCReader/python/HepMCReaderConfig.py
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index 28bb6147..89b64e1c 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -140,6 +140,51 @@ def FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs) :
return cfg
+
+def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
+ # Supported keyword arguments:
+ # model_path (detault: $PWD)
+ # model_name (default: DarkPhoton)
+ # mother_mass (default: 0.01 GeV)
+ # com_energy (default: 14 TeV)
+ # daughter1_pid (default: 11)
+ # daughter2_pid (default: -11)
+ # mother_pid (default: none)
+ # mother_pos (default: CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0))
+ #
+ # Note that ALL of these can be samplers themselves - either the simple, "literal" variety or a sampler object configured by the caller
+ #
+
+ cfg = FaserParticleGunCommonCfg(ConfigFlags, **kwargs)
+
+ pg = cfg.getPrimary()
+
+ from DIFGenerator.DIFSampler import CylinderSampler
+
+ from DIFGenerator import DIFSampler
+ pg.sampler = DIFSampler(
+ daughter1_pid = kwargs.get("daughter1_pid", 11),
+ daughter2_pid = kwargs.get("daughter2_pid", -11),
+ mother_pid = kwargs.get("mother_pid", None),
+ mother_pos = kwargs.get("mother_pos", CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0))
+ )
+
+ from ForeseeGenerator.ForeseeSampler import ForeseeNumpySampler
+ mother_mom = ForeseeNumpySampler(
+ model_path = kwargs.get("model_path", "."),
+ model_name = kwargs.get("model_name", "DarkPhoton"),
+ com_energy = kwargs.get("com_energy", "14"),
+ mother_mass = kwargs.get("mother_mass", 0.01),
+ daughter1_pid = kwargs.get("daughter1_pid", 11),
+ daughter2_pid = kwargs.get("daughter2_pid", -11),
+ )
+
+ pg.sampler.mother_sampler = PG.ParticleSampler(pid=kwargs.get("mother_pid", None),mom=mother_mom,n=1,pos=pg.sampler.mother_sampler.pos)
+ pg.sampler.mother_sampler.mass_override = False
+
+
+ return cfg
+
def FaserParticleGunCfg(ConfigFlags) :
generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
kwargs = ConfigFlags.Sim.Gun
@@ -149,5 +194,7 @@ def FaserParticleGunCfg(ConfigFlags) :
return FaserParticleGunCosmicsCfg(ConfigFlags, **kwargs)
elif generator == "DecayInFlight" :
return FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs)
+ elif generator == "Foresee" :
+ return FaserParticleGunForeseeCfg(ConfigFlags, **kwargs)
else :
return FaserParticleGunSingleParticleCfg(ConfigFlags, **kwargs )
diff --git a/Generators/ForeseeGenerator/CMakeLists.txt b/Generators/ForeseeGenerator/CMakeLists.txt
new file mode 100644
index 00000000..e2e5137b
--- /dev/null
+++ b/Generators/ForeseeGenerator/CMakeLists.txt
@@ -0,0 +1,12 @@
+################################################################################
+# Package: ForeseeGenerator
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ForeseeGenerator )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+atlas_install_joboptions( share/*.py )
+
diff --git a/Generators/ForeseeGenerator/data/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy b/Generators/ForeseeGenerator/data/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy
new file mode 100644
index 0000000000000000000000000000000000000000..36320c3f1c9cced47c9307b48529b444a2384be0
GIT binary patch
literal 59240
zcmeFad0bBI`ZtV<3aLbrR77)<CSBJ$k0uQYjT$6MrBIrcu|Y{HLlKceWN1Xm5KTfP
zDWQ~-AtWVI8ob><9iHb|`?=oSd++<+`~AH9^Zgv-IjnQ7b)LufT8Uk3xyZz3<Ou%}
zn<Sl_eO7r(YKTcPUDPGz#3WtTd;5Ajc&u3O?d1H2Ue95jkMp43$JN2pc~F*`Cnu()
ztU6CxPHc<VKlhJYq+D>c{K*Ac-?yJF++xDeN-BGt)5Vypg?hmp*O8MYP5(iUc0`?+
z{Z=2qNb5^W-G}^pleC1FH^hhCnDsNS!^jDWp>wtVAuqN|w_`qRe>B(l%i|7O{_F?c
zAN{V|caSL#vMl}Rd2zhFS!A*0c@>uZrN^o6-+s_jPyA;)|7<+Y(?8n}jHi$L&-MrB
z7xVqI{lPqOUV@i5EE}oy@9LMDOJsAq3E$84#kgf&d(BJu=Q{gW+rfNzU2z`%)x0<#
znd1b~*M^MM!gl`M`TweZ<Gf?P|EhL;HD%lfo^oi>=jpF%58J`<{Z;K?eta(e)$@q=
z9esR%SNrtqh~xjO$NirgZ}HSKT(e?(HF4ZH-v9J`|6T3ldc*Pmr>?_)OFwWu{I~Vz
zKYe}vTj%A!?K=JL&yRBi?&op;^n3mMIH%FiA;0J6|2$Xzo}WK=j{UuU{=d0?=yezR
z^O7#Je^}MyGH2O8)L-!Wqn%Y8UC-)1!#sFyrpIA>xIe`@4&4v+!;O=kXE=ZCAJ#{H
zwH-XKv#N{y*>?YI-3$AT=Tuhx_|@@XJ{%WbFM8fTn-BLdm=EXqpB*o~KmTl=zxw<+
z+<gAC{l@tm?mGRm`Efn{Rr!ZI592!dcU@=PC*nEuuO1((bJgFq-thVNceRJ}i{l>d
z`u<hx3hxW-2VU>Ldp@uq*gnq3-*r9yQ{$!2JC6JR@BN|QAM|zcchAp%%XPx%Jnob5
zxj=t!(eE3K!~G`uqm27ToL9_Cf1jWq?mux}usv)C<FKFfIJ6J9uNm&W^Jn)}n2+Ai
z;rc=EKdwidr$3vIKJTpV`+rw{yk6LEyv}&ta9;l1?csfl&-cHpUED{~*XLg~&cAvd
zaXdJmcpY)w{9Wzid}8~5b-OrTTtC=8j{EOgPxOA^`2X&H(ANQdUa(((*SrkZpa0bS
z;5xzi_)qnN)%E#L^@qNW|6AticfWqEemdO!_Mcwwzpvx(>-fj-`TpJe$lvS!$L~u2
z?0M(+{rLNS{Cz);)cU>8`M>`j<o9#+_jC1+-+%vAb&=oe>-YNly}o)if7b#2-+y=a
z`#JIZIq~~B@%uS}d?B45L;mXbbK>`N;`ejn_jBU+bK=K$;s5II4t_sZ|GD4!|9%hs
zeh>YA5B+`*{eBPqeh>Zce{bRU?>c_`eVO0y$KUVA-|xrY@5kTo$A9np@xSf&9MAU#
z^?tvXP0VDDUFGVzK<YN-Eo^noCABo<Z}RPFlekO}=UVh;EmOTlMo!otbTst_88L)s
zGPi~(7Zno3vox8me{}H<?}yf6QbR*l`YED{hL+1p38D~9&bF4@kU8lNp<@;cJ@2df
zCYW`P(2;K#{kh;cy>)RFBz*|eSjecSV@9*t{`HqC31Sva&K#1{hY)4-N4@%3-RqAV
ztBBeVqWotO<JV`nhLtX=A)Z6XN=7{^%uC}}4jlcIT>eG8GUUf9?(V))Zv}y8q>P4i
z{ctkIqaWgqCz)9bf}fL}Ll{BB@F5xf5Ydk=)AhT~awuk}HIc|6M467m+0mcwm)w2i
z#+~YB!U~TK*%5!$PesS+8)&{Jh{{7U9aU+5bbpi)>G2l>Vq}cJwve15%pFEXKYASM
zFVW)Zep*@X%O%-b38FSlh9OyV2vJ6VL>*cj`Xe?@eRfI4<0E-KgiSPjF(jiOD{P?U
zq36T6!dWScdyaLG;jrKr{f5g!x7SLn<}JS6NuJV>u16UW{pm955z+pHmS5jb?nI$d
z7cm$@R<a?@5B(8`v(xi7kK?tJYUw6TLx?isiy?jE5TeWqF^--O{posn!-I!qKJ^lX
zA(Uf*vO{|GM?cgfq8^cz9@}%06SWjx_MJEkA<B+JG9A&6u16W|tA^qc(T+&hyB6qu
z(2pMph^vQW#9!6Z{phkAtsT03gTHQC@ijJJg#j$=bp4-|S+&QV_<W<$tkJ+V45A+*
z-Jh=K93~$<{-%^f#`ve4fLKbCONL}bR(`0b<DZS|JD?<~-ZU0iVgHbwj&%P47I}JU
z@fe5s(2nsN<a2@r4JJU~5c<<(v<D37S;@oc(T@31rsKyrS<b{LKKNB^qxqvuM^^dJ
zkClDS^0B+`W=#ZzA!H@f_2@qwVjLpfpRTtJdbPP;a55~TAzg1nv!g%VkConTD4y<*
zc~Kv->bhdzohd-aV46OZCe!U0PuHW2I9#4X(=LBo@<138hVTeYM*CqFc8o(ml<A1^
zbUntQ-K;Z3O;kk`5G`mj>dk3-lo8Pn<A>9upD8WxDMQgC-NE9JI)tob)FYmz`46Z6
z)%=_twmG#aGl7eSsHfwYA-`YMqu+4(F|TZN!O0EWvq5GEhm+}k=#My@9sLl|A5k!G
z<<rh9QZQu*g@zFIh?7|8(T+Ho=Eo`y{jgm`tk=-%Bv|*L*Bh`-fpr8{^#S^M9_s+~
zb2*-y>E~QLKjS%-j(GlLbq>UHA^jXjKfj@!er}`7cuvFf7$TmR=;tcb<2en_Sw;V*
z^Aet4@H~Tv=NYu4jCMpke<0$y1(B5=?R3O*1uLYVH&8}@`gw(ZUZBf(u3&|DKEXKL
z?<3;*01?jtbQ%2-@f?6i-?y`p|EwPO;q?7G>i;a#_t7X1w;%u2{W*Of{i`y4|4xrX
zJMMQ8alecEVY-ZZ+#h2cU5_%_5pf@j`(#AaBhvNs{V?5*afs;8YQIX~&!SAnUyWzw
zr^sSIOy9rC53_$pfAqusE+X1-zln(ZR@{%GjEH`Ss2|Rs?uYS+xKBl-?@Lh*VzJ*v
zJ3SutxW7e2J>4Jone_eVaQj_War8VWqd)FvQO5lz?lVzFL_O|XEou8uw4*;FeZPwF
zzq+5qeJS1lSM?Z&h;b)r`(OILl9i16O0*-Q9`~iFM?^g$>RBP~9}#i?NS9HMC^xjv
zMMOO!#-p8;jQ)tY-~4A0_nX7T<31Jr=y450JR>6RUlEaaM1GNud_(-B;t=nM{3FW9
zKMscjJU{qK<RMw%=^@^d&R^1b%EvU`679%OvO+p<$;uD)xUc_JWW{Dp&1Sn!+3
zgAQj$eiY9Yh<NTmew2=AXNBlbw<GU~i02rT=}5N^Cu2PNss2em75P(E=PK0C|HZs`
z4#V>oBI*(8dXy2-j}`tZKTG$+^B|qKMLqgo9pZKIyokt(ho#%m5Ahz2|3y8<A)+1;
z`B_$o=U6%-?@QOy^PnB|c<x1>77_Ks;c))L*`tPdW>!2e>hU~|NY|rGN4h`Vew+p0
zjOTkgUp#LZ{utvB>GAY>0?LROhp0Bx9`ekH^g09T5wY$-mr+lTM?Ia7Mm-|tV^zPP
z>*+GqIfg@e9M(hVh;<RHgJ4|*5$hr-)BWgn^h3n_=!e+Mf)7WYxp=4!^Q-lk8#I56
zW5s{7(qsE{#JbF{BGz|UA?86l%5+3K)`<{Nk98nKtozVq)FYyu&V!>K{Sb%K)9qL{
zLZsJ|Ueb7UR(`0bBfTC)m$5z6V?5TQu<nF)DAW&!bbpk8H6JTKtY0Bwy$fYT<lV6j
zwvxubqaG3CFdri75m8T%V^tT!`r2?fTs$i~>RDkRjbHbt$*k&fbp3GhzpIXi^*lOa
zeGie99ral6Lqz>>NUz(m;?Yr$bv=}ai$^>9V?62+QIClAJw&VrqKx%EI-;HKk9tH_
zaj0kIhxI>1tP>)l9_xcB)9rLU-5>3AJ^nrb?TGaHAgelK`w&l$c6uDj|EeAH(Cde+
z<X^30VjU6@^I<&df3=Q@{5sYfvA&2hy-rED4<}<D^h3nFC?jIM5$lYI^!g+Ep&pU$
zKb#)((BtSby$*@>N4g$mI>rz2{PZ~V!@4FSU5{~SM;Ysq!y(o+S=Bq~@#v3wdK}91
zIw{6sy^@Yt&qQQpr`I{@{%A)%B0V2nrq?Zpt7~F?lpZ&njCm2!pDv>w>zP=$q|2yB
zq}M^w9})c!Q9m5g{jr`&M|ypeZpV7)aDMbSlrbJ<dLFE!;`bXU4~OVqPpijb9NH1l
z4`sR^UB-H@>CpEih^R+2AF`uNM^<rkKe}weA}`jL5%K#Llo6lNzI#FaaESR>A*;Mt
z=blLWz6SMJ_eOa*#5{=TPtSug`XQno>+7sgiB^ZF%k(&`)1w|`MAReF^(fQ7H|ig%
z=VRTUZbumr{Sc84K$(uL;?VE=FY*b<Gtl`2)U!gmKV8N+R(j+k=sZRB5N|>M{Sf-o
zksgnBR{G)cBCmlcIP@JC@*;FZJKZ1Ue-*L4;rJCguX3LDT^U`b=cC8dW%SP+%8Pd7
zWsrA4WX0$3)A*g?_#DjtXOYhPpo}~Y9g){TJ^G`bE~Ea>BJw$OM1BXqpF<gu?vFC7
z@AJ^kDxR)K9te@n2mLCK^s9dK{H)qVeu&N+(Pg?H+7Z2m`F;@nFb^w4Kje{+S3*QR
zEBsYIdOV%aqT5lX^Gqls(s?S>)A7y_pM`qltLX8lM?XZ=BO+ggJQO185$Svt{%#ZX
zi0F@Y%!hWA5$SrAk?%s@3K8{)$Xl_pv(gV2hjv74kIq{muZ8>-D?94ZA7w<;BceWu
z_MIlSgMNs}Qz6oMER^Yp{&b}Kp^Sbg(~)jRo(hr9ccGrnYq64t)1#dgV*AK@A&*7p
zxll$#f4Yo%`uCog7wst1k<h*eMP93(_T4Aykq<)|k<O=~A0qv`Q97T7eys8!kA-{~
z9a-7`EI)>MSRp+hE1nDW$aA4gN94Z{(N339kBD|w`H{y$9FG4&J1eB;VdY2X;pl$I
zUm?<YtzYHS=yCLXtavk4dgRMk@mi?I??`{uALGz}I7A*!X6XA@<h2ly*Fydbc{9|b
zAM#)*4~MM$S?Q6_8jd%^?@#G|^!Q(u(GL;xqC6bWMd!)TpRWJ2GUj20!|_^lJ1ZUy
z{m_p586w($6_MvcUW+axzeU&6WmfU1XNB}U*dFp(h;*I{e}|8HMD#~QKa__<x<AU8
zANeap<h$s68Op5e7>|A^(-G~M2kmrT4f(5o6_H0nWR-6?JH{d3g*?}vl`)PLVm>;e
zKPx?*&qCe{^(Z6H#fo=BJ0d+E^Ux9fFb)xUEIJ}DHXLFc-4FHneE};w>d_x%<ge&F
z7xG&NS$r>xadgByh?oauI$}Kfp^W?$o%cc+5&18aF^-OCr{mJ0?-0-*?ZfeDbe;=k
z<iqHE7upfgPM0wr^>oBIIwF6C$co=WJ1ahol^){|hs%$C$a5i|MMt_FWmd>4{#SW8
z<gr+NH;eohtM6dZj(O-Z>gk9)8u}yJ&s)}@6kAWStD-^!T6@Tc1LIt+96QN`+B0WY
zD|C_8GM;%}QZLCRCD-JOxjjTz$YERJ7*4oW_lNy8xonZ^yQX|9@bbCvRvWX!9FttV
zO`HbACMfgW_LrqYd8I*t_iTyyi5alsvS`QZbrV40)Ux_>%mLyNc7OhhQ}4-ujMFR&
zkEx)WY>|`TIuRU_mj$dE`-1%NYa{ttQ<PT^RFfq5w88S)7+7=Mp66tO7$m6t=zm!5
zAqzjLYQnvp&sIsXk4&@YZgrSE5n^5?+-bC6hpTKIAJk*Z$+YfiI#kU68dt+n%cb3g
zvmo%+JJG8-lVD{k`;m-0Q{nab81K5wM?dG`{GEm0+hV^D^=E*5ve{;|xHp_k|8^(V
zw^9P8IV>IVOi%#Y4__0}oWKR^6U+`Bye0*&9Op@>#Eyr#iD&O-zT+JlcN<}vpH$j?
zEs0D~wwF)!=pa0!ZOty&y(0ZihKmBqWWd8|LS?b=cM_lN$Y`_>0I!>G%2i@1;IAvx
z+u=S5eq5K&Br5pwBI-sh$yj6@7&rL|QQ=qczG*ZT_H=Fk#I%?SXD7dp6OWPv&Z(zP
zSlp3?=a=6{joAE&Z2S~r&b>(-rp9fGa6G{WIKMTLt`@Vuz9jk&QlfhbzmVp6o3`}q
zt|zvi&r}Jn9RVwRQa<L(tANsuh6Aaq6=9{~jNM-{B*FA#WZU~!{E(fN^V)IXCE-!r
z>yqCq&D=X6X1MoCH}Mqc?o*T$1(#>5qrUIXC%PTwUNA-ha>p$?zu#LH;Pa}p4;$x$
zk$aYRevk~XSt*}7d`b9czPk6m5oc^4Qq#Y@pS|h4AUUlXvZwQz)z5Y2yICN}VeNa8
zT34lPK6xB~fcWa2z0z<lqH}^u0|PAY7`pnkYrs>zRXfJK(}y>^ADzf}CkA7!7c)59
z1%YXMHER)sk-1#{QO}a@5$>6c;^u;EvM?a{sEyAH((z`RM%(zwP?5XWS#~A^MmsFG
z^(dBvd#>hA^_jEam51W<sd)>a|AQ}CV8ei~QC$x2Ud;erj-sFgTRxK6-G+P{zPM3p
zd|X~zNAD%_e*LYTCC`bV*Hn&#tJC1}=XmYVs4l_~iTv#LMiR1Wk6Sf>0Z_heh5b=#
zAhmn?z4GHGAX*xnw8LEm0-xk6xgM7PnZG-5<M(+L=gIkL`TCky+lgDg&FrOa6CqP%
zGxV<J1H%ZfYVYNya6&!J+*4l#s+}S~Z`r;8_<EbIUQ95DW677s9(&7xDUYptyrd-H
zUg+LN#iL^&`rT&f<+7itSvGklQI38j$MBMrVd(=hFPZt)(6Ew(D(Yz;RjnrxZF}D@
zbesmO()+vn0+nE$vSh#W^m$<HRIz?#st&{*Zt~=aGXaeaCF`2|ResLTk9GHw06%Bt
z70;_o<?#7lM;?Vz(kZk01B&(%v%C9c$fRPTy^>F^J-dlKNo#Yc44)3;dY;C_%vJ_p
z8>JJ*RVwh<S4FupauEnkzBqQ%^QACH=+xX((aV5eIN{@b9|8{7Czkq1Xuu_f)#IH9
zK9iFX73_x&ucY$sZ`m#w*-fcmaK5Z}`ZlSQKfS8*!YG*K;xp;Q%$bnoJxcogPBpMT
zcZlzDwF%hdz0tk+)C4lehgrm4bbxxt+}j3o9Du*5PmILqfz}PvOQ&>|;FYa=?Ev+S
zgc|obaVbBc8rSbRmNQzEWV{@x@ZjE0<OXDJGm5H7f=PhSc8wm=*pPJ9>gr55C-G7?
zvfK!gIoK1oi7kYuc|s~diYwsNyPXp<WL81xp%qJ8m6yT7-ba~Nq%5HMV$)c~Fg++M
zQZJt~c)fBz#QVt1l4go$>08Vn=SjsV=zOo+(M;{_HdWQvqR6U;wcl2+K1TdB4_uwc
z_n5ebGiuL{=YmT+cuvlYoC&(6Pm8AS)q?xq=gkPHw1;SpN$$qqY+&7<+fPRlXUH@y
z@S0ff4EG<ttF)c90A@bfKg;Qa7M!cQdZBM{yzRFug);q{m>GNe9Zr2-MvYb5wBgFV
zFO;HcS!n#e6=bdbsg3P(ONkP@ut{I}C!#9DVf}T9BCOrJQHL6@3})IBO4%Z<!TQ7!
z*_AeqP#CejD%8vc#)kylt(oQwRl=M0+L~I!$jW&Kyrh(&UT{@%FU1WSTRplW76nk(
zr^cRb?)gluJTN}-`zKdYx<Y~Ru4Eqxaj_bwHlvcvdz7O3dHHB)+MgLTT68Wf;n7wp
zojV`49M*V$Yn%n-TT!j8@@wJiw~_&7pC{C<b7Iu<ZiFQ6+A*JESHPj|{1k&%6Apel
zuqq{29M(;<9WNcL%;3AQ&-BCna^|u5@%2J|oW!`c;_}|6CFDbh%f<7r_L2!1sfScF
z8i~o>b(=HB3V`qOtH*ASm<Muxt$iX&wh*$0pSS$A6L^>e&7Yg&1yRdNo5fFfgSGCW
zj<Anj5ZUD=`*ySoczbE@3QpI9WfO%ravhL{-r#RieNov=9k<;4ili`#F{(guf-WCf
zY(FExx6YY7(p9?jy)BtM>K0YfROlcN4&R+;aj%z%d-A0jn9YY%9Tyf&Qn81gdh@~y
zTvkBwk!1z;f$(Z~X~X6dzEESGF0Nm<3C;yeW#?{K1K+k@^>7he3hWz$#%IY-hw7@M
zX12D9%sm&_Lv3%>GEJf?7g_2DQ2lbtVkS4T6WduuANhG#k_yfI`{SHWlRDGO&TIBP
zA+>7pwUg#f0*9&Jf=gYrV2Rp<%P*U);F-wYk43}>0vEsf!Y$|rr3^)RL%&VH&K}E<
zTo?eB=3^p0FZF`k^}@VHv6k?Jr^;NgO$uTP+KhO3w3+fR*9mMi=wzzvbGeTw45H+p
zoIZA^j)SxdO><A&=|Q~uqN1dyq>`L1`tza#`iZ0Zwp0j{f@kx+-MN?Rz;u!9Rk4mP
zu&C_rGu3<E(2~)=s_pqE2u}6c@+M_7eEljlCCzmggr3Wj3h3~LyUGdT)s=Q|_kW(3
zg5Ylxb?ERiZ3f4b&--@m6lC67aDCAv&JmPvNE(-y&^}7XO*+H!FgF>w(rNF8$KHft
znj9Ckw1mubwYWaBYz#;w>`G(WNI+`njp}!2^x<}*bcB=9Dj0KY_lS{my<tq`<0$i;
z+d+T*Xg8DBTVb-Pznzw4B-|;PP@m?!5gt|_QY$tyg)1L#h`9_XLu+mIgN{pvjBzKU
zpFMlf%CzuSa+SCfOC4DB>3rg94&rvKulwFYU!tWJ6=ua#OdfuFryCeN_`Whue<;6p
z8oXCc6gBKJhxCnyB64|Lz<O-xiOpmi+?i*5P+K$<gbQxXad{mJ*^cx1r^ki@lQAV>
z<jbwl61@D`^}8-G!O>)$+gW9J+_^PKZIm3t`^B<DGd|B_Msf8OSIB;57P8G3_tDr%
z<s59N>N+`!tZH0%{_#;Ka=Wbi^!q7gMAArGdaryJ+39ybI{1kkc&w{>d3>fN%)6Lo
zaFu-}Y$aV!I^S&vVXgPmkB^Lkppg=0@2>9y@2Hz!+V)4l{SOxmMDzyd>FCB|@#@Yn
zB{S9J0?%AH^tFEnnP$qEv}Ur{lzM4qqxH#x9_=yI=L>UR&x|=r#U_`?O*}Y-yiN;Y
zQ-?ru=)CI4T(uGs+a34pMfqrOIdOe-<ux^E_rIz)s>uWr?^Sxg*t`*XHoX^dtBiyU
z*QSe1_i(uIJhr7S>L{$9R?ed%atzkaSTsiqg5Y`Ck}WRN)`3-!@4*S$h9G>Xi6fy%
zh2gPb%MFn#ZRYWWJrf^J`oa_nY(Fd1d6?=lIm%z{Ek-^#mTI;JMiDpvC7t&kJtS+!
z)xB9|#tG?*(+q4}8Bk!pV%v0f7dUxj$)ughfpEQX)AP&byTS3;uI|@IkHCYScbdh_
zk|E1u$M;CxW6;)gHu(7bP>|J0%+QSSgWJ;+6-k#WG)^ow&UyS@W2;?&o*$boGghoa
zrnC1ebL662#VWPq)LWaH(K^@Wl8~g7!Z~Yq5)Ru-+9J^pNt0JuNZaR0;Qm2C>zNb<
zqqn(By;HXXfz_Q>yzfI``_9cmy;V`r&^dF%@pDO#xa!@a>HH~>Aa>y5yt?CX>jOvg
zW1T2ywwmuYeuFo7l$^<Gj%L7niJcoGADS_0ZHZp9pBA(2yTmEM&urA81rsjGUrD5{
zpLj}S$xIULlkluXbq}%gx%kN9TqT*EVJ*6;oDUMlTud1V(}ZoGO{W#rc)<ywcFWP8
zk>KL15h|H{0KS_khm-`KghhJU8JQtR!FS7{)vfc7L*w2<XZLK21GTu7#nCxV;QXz?
zWPs8D@dme|+Gulzr;x?NP0a?(>M6T+3m%<9t?kP?xN5{@YU@MZNFgaDBGD@uzDH>n
z`S|M67oqEK$iV>H`-xA)V3B{)%PZSWf$<=JwNix_JUaU7(qTpn^r^jl^{C(|_^)#g
zy)ix&bfP6iPTn~K2B&#clUh<h?@H>xNX~t*Quo``=u2MUR5@V(c9juG7$;8OlVrzG
z3b2$Fk1}UkNS$^O3>Bi-$1R^TapDzfWZRpX#am@b;`pL<vwRPb^b^7L%XL~vz?WSL
zyNw0G+;VmPx?VHbHg+4^c@0k}H8^s)GB5&cGgJfCY)pnJcdk#ovpE&Grmw2yo{$aQ
z_h!g%tW1U&sk?3Ak3-=?fv5Vu^P7Nc+@VvM>RM1?a3#k4oD1XnhS4&9*OxN8;>A*@
z=S`r>BJzFQ_7zbIj&gCjWImZP^Z3~Jb%{i_DYQZ|`6W^IIJHZ%XeK<LTD<<GnK8_6
zyxGj9><b4<N;`Ra_JZYfQOi?O@enO(&G#}e8|n>B1aC33Ao=mD4G~7CfswNCiHBc2
z)X&uq3^?fyW22X|*EKPqd#X-a{0tXHM5C+Iv4<;|24!!Sg~U&%@?1|buR9b_E@GFR
z^x{;>5=NHlPNf9WJ(7BCy7Uu?Sd@25DStNHKO<=g1J*E^LGh=*2n9voQWbC6{cvbX
z_Tde2M<7)!-RJm)4B%W*5_}~l6WYJ<d=H+G0p3Z^s#{kdhOSxR4No3z0v{%ie3CT4
zmc3?<yNeweL8(noyI(mo+uY;k7u}yp8P3w<?@GH$J$U)`LXwja8If=*>op^V*q6?B
zSMd8nBFYlHD@9bmbiBl;HH{A7*TAEh`#Ka}&;4k=)a4-TA9<GwS(y!NVrrbi`njN0
zd-9}Bbp}}7kiQncBNZ0r>%@G&69+1p;nkzoZUmj;#@BOh>jODFL!w<TfU&IKE%!pk
za^^Zpj`;LQQr<p_WX{KXlpxyy$La)Qviq{i-F1~`NJy$vu8LO|aUE+nYHf%zC`xTg
z%U@>!TH1DE&qs%XlK%8G#e>ISm8ilK#qrtTa8}0Ym0A`gB$f&W`d@=dEj8NS4Jm`?
zNgdN9yLfozX53M2<qv1xe?0wt@cy#w*>&hew<jZMZsM`4`Cd%#9mQSyf|aPL(Xabn
zC00;(JLfyTJ!M4F88z!n+tUg2`p#3eo4QDvcZs3CoD#f<mK(!1(G^}B$kjg;+yQ6P
z9z^eYo(!U%o8^_)Ux0>5*Zp@{WkJD#M+t$GvO$zt_9W{3N$3btj53Hn1V$R`S6;vC
z4-E;G5-ZJA;P97^Tko1~W31T0zo@-(8Pn|jgl|=P64bq>wq0fOA5)*MoGVn*H75d=
zb4$<apCzRi#%5M3jRa+q@>=E+1I~mO%*wa*f(f<Bd3G|vFe+e-viGf{uyc&Fe29Jy
zFmF5<*&<&6%{gC1<BHRv$zN%{#M^T~ao*q=h>HXHQtOetetuB!GG|JD7zOHjlbX!a
zbl{YBG^h;zy}H#$7{Uu2rhvzy8xuJLc)-7YrXgGK#36Y)>|}i2He-|kz?oK+mroai
zL-W>Gb-tQE{cAi=9-2MsGYQ!~Fn3nREQnarQ^;1T1%4WT#IZv|v-pZPyyBoR)oMDY
zlMqOSW!*~FXTXIk{!<Mdq`<JfqHk2NA}sw>dQXLII*h&jUf)}4D!958Q4=pKL;LGC
zDNd8m#D09g7^kQ@JZ|MP5<g-}j;7tL$gd<2Jk;*`q5g|o_(<?u>%gt23M<NM>>#UX
zS99mpg`gu~Zol4r8Pt_=T8F3WzyYy0B>uGtnSCH(v3aKp_{bD1b7)$@%l@gt##Kh(
z9%*yG|D_q!vsH+3&7(m0(Y+{%!Fh<t)l{7S&<MU=Bw~_TcEHf#e$F`<_sNCVXu16=
z2+UgB;4#<^&sX1>Tx~|sr1!u?$iJ7glsR;pe^?H#2hT5?VxkCq;)?Bsv$w*z_4cD(
z>(+x}U00+(V-;*mvtBdGb``iz(LBo8#(={0f6OmC)Q)jX?eWxsM5mr?+l)ofZ<_c1
zYk@U1#Fd8T>aT;-;ZZM4{oTRFPCtLr;Jo>)6c!k5A^`3)=Eu$&r3@8YCq_Mwt0wy0
zrp!?tvY=MNA91>h3C(U_TN1bopr)v!%Zk$xiVieu)R~&VKo4K_QY~MQRi1CV(7_ts
zRI-P#1#7~W&Ofdb4@?;oppv_v0}em3m_1_#1u@p*Una;4gEl|!$A-zypv?23CPTvy
zwiLuhCw~tI_R8m%YTvE|jquHUyt+;>s<d;_tL>{`gpq^l;N=3scTrPSP@1V|xkuh9
za1jK&|ID5nX$rT!pGJ?&_k*Cy?Z>_DZh`a>^+%6g-T^zTOAdQ@ZiZINPkDY^cCc(;
ze%i}#E1}}r%-HXRN^m6o(jG%Y0r=W5^7ga6deF`n`?54o41Bd8r@8yNz)_*_-O)LL
zuxr}u71M9|Kv+|aviJ#S2p%9k;c81^xn$&1zG7pzIM(vTJ(vB&!Pw<`Vh$&$&bGb#
zeBuOP+r=<&?off|agzo-CTc=cS8=gNf*t7lRW=JQ^oHpfUsrEXHG;51qu$1N!~lPk
z#8iu0{$Q_UA6jo82I8;HDn~onz}2xMveHKx!q=HpV!c2S`DmX-y&We7U*k%4$CYct
z#hcp+IqnA!+Z$uHeccAeTKP-Niy~mH>7{9lsx|=?{=7gadMkJ_Do5QpvJ6~4_=Ov4
zt3b+HAD-zA+RSqH2#>9wXM;qn`nhNsAy_wXrzOmv0g)bZ%3AG?u<^lI!-RsZu&74b
z+WGT7n3=hH?@8NjVE2O4MlfU>M6b`|Ejt_myZee0oDv;ihHmoH>spGSR@zv7|H}(v
zqOfvOuq`KegtS*TbvS^YlkO%ni+ND6a9;H8=;d%KQeg3nSANhKcPJy@Gzv<_%Iq-}
z-wVF04Er3P?go_;Q!S3?`oYFWG9qs~JRm4IasJy%Q<#{){L&^rX>xs5hxEn?A*8u>
zeO$9!2^r1Rx~hY79MoRR*?7t93%SyMPh&~8J^0Vz@Oio;5af?YzYXWu2Oo<Dn(ylG
zhkj{xKFjIRz_xGEfyZ+DpgaF$#y+i`(0J9s=~TWmcnsYi=AVtLGin_e0v=%k#VJ>*
z9-HS94-Q--x0Y<<T|aFa)TZCEd;44o7M)J98(*gb3g=91UdlMZa;<BR?dFF-%(`Y9
z$;`c=wPgS2?`-EFXoVqDs4ofvh-q7ccqC*$45aR?T>*i4>v#-L*n`oHUjHL0oy=`b
z6Z9S)-a$2Anf!U5Vkh-{8)J%I>`r2rGF83P_a0ez?C74(sxG2D`M#W5nHJbzASLsy
zH^S7FWMjm~K%l;x1<LM-2f|$!v29%foS5=bPA)0|<_90RE>aK=r+CtLZ+o!=#xKn)
z&zrvG=Xm#pZ?$yJ-$>R(jSl2Hd5YMnp0yikeu4C;^l#K!#SPyc*i0YaFA13oKX)V;
z`amw<wR@ZT!$3_Xx7g|1L3n>&{he}9GK|*Sy=3X+<Io@+6_Vr{2kWOLu3yI!2k)Ky
zzsw)K4Gb=<Q~e~a4T5Ip+R6;P7zY|vXUxd|$c%Drz4pbUl6o7eL6)Uhlee-Rkqz8C
z$(x%t<EoaHkvDG2g~shNaG4!q9!)ZVGgIB?z2D;pS6$@_JsyO^K@m&ocfE(8iJw_>
zB{>_$_f$TxzJC^OKab|#BXWFjUrvcil<kMjVS%S>Mu$PiwS+?x9ajO*glw;MlVq4)
zB_SMl+-6YYu4Fhap0|-=lf2>-du$JRF;!ye;<?4-_KN*R%3jyWhiB5|4`W2(8E?`V
z&l?7iaHe&wn6xj5cpN@f{3!@bo3e}h)@FhNn^!>Lu~cwdmlRfAkO61yaz0)d+<zZG
zymimn<b&XGsNCjD<=}f*WxtJUh#p*ZU)4~4T#0Gf&lh5K(~puHaU;Z9hn-mYKby?G
z=m9w}f4jg7<IiNH;}tgrLqpKuoci4O{6@HaWYh9x<9uMk+sh)?*^YsnT;B833(vqn
ze|vA#`fS*DKQv>+_$*l19CY@Z#YHHaviju6S9{=`tVPuvkDbtMpth@G@Z8zDCv0Rn
zQ-er|i&5LTs)&~{SCC_+5WGFRvM{Z6A+YaR=x{v345ati`h?7gfZHo8r1uMkg2C4K
zGi7P1&^vx0>tgwN_%eZ<lCwVzS#B>s_%Kd@+U+%al3i2a?ZW4Kx|RjNt-`rR5~D02
zNJPk6CC`v4{q<U#F#jxyDdyS#Aa)J4-+RRJCUbF8t-PB{QpSqBf8L|e8~&KEE12(|
z!_`lCg|?R*C2An-|7E-HHcv<x2rTqF8443?mQ)Ai#llSX>#wGdJp*}X<4WF?WdhV$
z8V+pCf~1SWdpviahV%s+-k&r(2uYI?l^``7?k&2{Yz#I9`wG6UrO#$EU9v<rJ|p$a
z^;cHD8e4jaawzIAiZ)~;ek)#Hd?~Vu*cM#7-*T&%EYRnFJjZS_WS_oSBYM~j`aDm!
zKPV4?mcqvL?X7!ZNABD*&;7?hwZe`c%JYDs;E+%<?J|__)=Iq2aS_DR(vC33Wr9-l
zvyr>plL6kVtLJukL+oCTcIy|GQ0C;s|HRgavB5+`J>k(OCfsj{_ZP~IiTRZz5
z6=Xl}av!ye_>WhAIPTO#l6OS@&b@U)FwP*}Y<7<s$OWtOzG8TR^nEoAim?Y)*pI66
zUXu)!AqO48gVNyBQby?~D1fno0#A~C@*pze#C@gG!(h_F+{3Xs3Ah&c+)bJv2HmQE
ztWQ;7)X6XSrpiklS>3duXHgI}ZEMdV4=*N>O<MKVsW*tkOCKncSvUCq$vC;+9j`9|
zY~cB@#m^je>Z-=XCi%hgWWVDM-ci6243Xm<GQr@fp_oc)9;iK%KHz)%JY*EHPwkY=
z0rn-iJPzu)Fri-bOY8ara6z%uYmCWua9`LVdT#K&n6GnTk3pjiv9vw2=bl9>*}%8?
zSmXgUQkk&lg!_a2WW4RGf}6s}NbsDLETNqbDXEz%TkMm^5{%1i6If^-koMExyZwhe
zy(Vo>Z+RwJ`cOH2-nj@$W8#mzDrRSuC-?Y}rO#fj*D~8d&`*60<66n?EFv*5Im=mm
z8WY<K`=9=~5sYK|oO1h?#if6UL%W2uTvo=GV}v0-Q`YUl$)Wyd6O6<D{}`wK&v9V9
zvbq%%31^cjc7M#rB|;w$_VY}!dxx$`78y4*A3yqO^Ru5gKI{kc)g953XD$!rP5p=d
zM_9DTd^&f5iSu8yL{3Mi(T&_ZUOkP!Cv<4O;)(eszkwxoSp@rOJUdmc%;FZ|i7AX?
z*xvk^x6|j1dx4}Zh4cSoK3`Bc?xH&%-!5nj{F(p9yk!qv$CD&4W>4k?$D2R<OV5w<
zfc?gHCJRSiJ6`^t5`7*!UtM=DsYw~Q8+PD-&GXicH>3o5Vp)!d+o+=Tn{+wR+Gs1N
zVf%>aZxkK(rM!YTYuZ~C?k*rNz8C8Fr)<zTB+2Ra`uH*iz1<(zLzt8;wU2p~WlC@!
zo9;@hC4Pw}xIQ}c$HZ;dy-IdHOjDf?RV24!GiUs^LLzcJz;{u=wV`%zkRRhLAb7p$
z>l4=#w)f+{a3lX`d)ed~J9!X3*dE?bH3mP%c`fS4b;EIne%X|B;_(x*PROH7aObC=
z_nXR@<-Sszm;W43i}i}PDgmd6?C}M)l-VW1*<iDJ?SmX*Qax(Li8FTy&g1uwTW8L_
zTS?-RgVQ$G-WnQD7P(dQWBp#vojPq(%4!X~j#l%gMV2qurSLi~xwqr=vSitz=Yl>l
zO4+cTkrGZ+@BJ~Ki6nKKW6gHv6{5l(aKz$rC7EG3eME_HAxYeEeTECvlJeZRntbs$
zWTa}^gu3P?@?(6rf3}}ALU8bYmt$i4Y*T-%pRiC3F5Nr{g5x%cdsL+28cViW!_tp|
z*(8V3<_u?+k-mLd{Tm{mkdq$POQ)9C{cP9gf5zQJnC)TR{<UxZVg1}4T0eL_pXjdG
z{HWW6IxvYNN#WeVpX2{=eTB)@c&L`x>`w5$!THC2Ow*KXiy5dVAJ#s8%z$RXtKkg`
zxf_YW<V8HncbkaFmydOu2Kyn!n5Gol`;1J@`!VlqQh%)DwG15phN1a)yiuXy6m*=b
z*g=kUs*NBx9}l)Dr01M)BE{ok4(K-SB)S|Q7joxbBQN4L6drymB6*t@`>&UOOpF%#
zN-2cDAzh`rRTb6@^8Wwxe0W8mQain?{5eTaf7N#R*b8#vL+l|flY*c9!Fll)n0wz*
zoXPAv*)*f=Tru;<efGcRMgQjG@h$>D8oStw#6<jwhF{9*Mbi>V@$2fsWnc4%d-19i
zqCWT@9qP}IdFdok^RGrZ!aI@*279j$t_z{X#^2mtwUYykTi-G-yd+!|3;CS_ZV&Zm
zg9i4e)k4;hOGt%zYv?-R{BZb6q;KBLMP7)_CAP^7GR^vr``(Lqxk**cJ9M1%?Um?1
zu6~p3&0{J(825zSRV>kJ5r0KIMksQ>-q%hveQ$|nU+g6_mq`nM&iF+3*CuRx_oSQ5
z);{y(=-YS1%KqacGpA}&EMwHV;(uGG=hJlSuSgnEKl*W=QkmdjVki5Rs;$}0YwWE@
zPS2>=aPfpMvEuf$m6uB<7p*Ij_<M@UUKi<`RnCn>hupus<;`2Nrf$}V9LFz2K>M`)
zU5l?oV!WQ-p_$)_VsO#A*Mt4|sy@5r;9x&+UVf~XWei-mqE`>6zdP`R`J+FKq5c$6
z+EX{{ESDWa6ck^LWDaPPWqB@jA0`A6vB|A-^fS+rIr}c^m0r3<++B`OtornV+~!|C
zAmG?R_Jo#tfo?A`vh?vX;vNC_R({jsr$#~4k~ap2>N&vg=mnd+vz;WY?f<@Req5&o
z%pH<8%(uby%*PdFg=Nc~sB|0utg`(%)U5q};oAnU*NtcSyI{bKM2AybE}hu^vtNhy
zRX-h<dyjngygEZu{RO%AP5jww@h?QR#Kkkdy`Q*$3^5}bqu|6+>(7tG*dYCr#+ty|
zF(B3b$8)p$=k@AIf74=fHCRJ&RZWld72l!xl4I&D(_G^&{+cP{nKiGZ$ckcQh8k`W
zyG7lYxif#tPJVKEW5XCzE*mnvWkK1EQBmX$L<UnSxuiz<l#k+>s-fr3J0d&&NuA2P
zZZhTSqY*pDu|tpeaRG^h(XjvZk{K~^W8vfFt?x5!I6-mAtC^?ne<dvj&w?!m*U^uD
z?bT4V&yg!%5X_MDlS+FQr^TG0AO74^wTW3H{xtsPV84E>t0KzvQV#pvu4&{;l|!WQ
z5-TFPO}$I+{BE+k*M1uKTq1$jblLY>Jt3-RK61`cYbV-g`NM+V^blwD^(Q8^ae>-f
zb(4fnUhtDB2#N9Ig|W|6?h93phmp<J*=KK#2A}QqleZjb8|qi^&-W_c|5df;w>W?O
z$eeh1<ookc_CNde!K5Z9Y^6A9+0n8yf43t!O6k7oyt|LAzkgKqmi`r@ADy`4yJ-{2
zxz2O$%E4~p&h!gg_kaVgJS<6w;u{BS5&P!NT_XTO91&(b9FrmF%kp;7V-q2%W9tmH
zIqWbtZ~iB1kJmrn_db^lqbR}MjH9Atxr4YdGwWveC)aWBnLplFP88Stg)Qx2&nT0)
z)#Ej4B}hw6x4zzfXX2*>D?@7&$-u0LZ{Hu665H1o6(;pGleX<Ou}&IYWPG(w?U*EX
zXf(`k7tP^?Vh<yO@!`U-{08rr-EqP&g3a*a<f&7@Z0ASq)Az@~fi?HKyrNov&X=_7
z*V1r>2nN?H=`o#4&6!=(`nuMP>t()Ad-7!BLodqO=0(}{1yw`u6IsF|qFdf{+nuN#
z*J;z7aGZqmZEgpjTf`uBO4W<iZ%FPJ*_WL>eWcY?M7+gs9ITD*xNdc460B%3d~e)4
z4MbCf2Lz}Y(3bx}N^*iI%<(_Eb^9XTf9RLm|E#aQj33w0n0dASQmE_TbySTyqUm>G
z%h0-iOliJNYJBS;N5Yk5jIN&dA~(A0qhCEtA(MqgW4MCvk?cJ8wN`EK34dJvJBzfD
zaC>!<6t@{agtsg&d!0WWPOu#*+x&bMJTXXqpl&-0Hl;k6!{iW$(ekefW5!Pakv*AP
z7Fc$WKK8Z7Z!_aGO68<gbq24azvhD{1}kit>XVI3&sO&QoUf5pb^%64^%OJRZqwvF
z%0#kHGt-<akTip>^{vL!WOb93Ub<r?sfw>CA1~8E+$A38$~26DNmFO+cw|2X%2Tvv
z7wCyXnfvuH`@1s0tvL|VN6m&v<Jh&U`)7fRn9_;1fk`00`$xapYAhzV4EF0pN67U{
zuE~rSR}A+W4)#ms^3D2PJ!2@%gi}Scm+qhx&nZTS8njS>_QC>IgR~4Wskv}aWFv`N
z|LNidt_(8qVzt)BV~<E-;bGrBzF){wZu?R94{?J(?|m~i&FN65w_{vJ-z=~hlXyN}
zdk);RTbcGLdk*YpzNu(zmxnY1{jD?XgkZ)p^W$z?dWP25F^$D~TReG|Co)zwj^gJ3
z=FF7;ex+-%Av@*Bw&147)ktcD%2(O1=^rSMx+x0pT?tt;eW$p)z&5gePs*m-W3!3%
zHEHGRXP%JRTP{X<r*@ONXO^i6T>NlaXrKJ@_Y$BM?AFS(Qh<gLW0y9Zm<Nlu#fw*q
z%>${8Y@Ad572%TpD0@d0aR_cJajt4&gBd~Rt7D_nHRjBT6q36V$B>fAy<)Ut4YOj?
zixX|<*r??t53B6rc2TQ&<zVd94yq_7YtG$X9a8>u_WcEU+sVmOJiWFxIpjT)+gC=e
zkto$XcXydR0=_G>-;AFr2(@+%5<ROVA?Tx&!?YYls8T2t-H@#coBKTK@+Yc89QWD&
z%|%L3_`P)HWFCp1^EEz_{dHhUs>YI&8%E`|oMAYA*PIx0ay>J()7wxyjE~ZJK5<Wn
zU@WzK+|v-g=1&yZSzX<C+K6PF3vphsKa8Z9_oZHs%_o;71LQs*ZYE<3>mt&pa=^=c
z^|FS+-$y3sKTSJ+Kn_e6U5wYOREE)MTaIvD)PNTG_Wb721lSg66}>Z7h0okpa)xUq
z;cmzu&(#5<vYR9ycF5DP;r#ezsc$aBYP*tjVDlEHRPFBaEBRBYFCl#&-OJ*r-j^I2
zU&H#SITth&AL}e8=WjL^Z+#g}_?DdTf2w|+Z1{Al$f@)#$(8oL)3lHq`u3ep_H7l1
zzM}q;ci$DE(E7Wq^dSu>uiY%7wMz>Y>pQuwOx1+&IElENZgrq~q^|BeJPY)M^Tr$5
zj)J@s?6Yrvx}fnUP<!s3n`as4R_D)Is=0~zy=HVn>b2=qLW`Vmi)Aup-D+rE@RW@_
zdMq@7W4Q@Y%w6E7pSFiA%}8`+&b>*@&P#1Q&-h3ZO6wMwtmXw5;?H^P^eiZS9(bxk
zP6Za~SxcARALREVxB5kV(Sb4itHX@v=z?|Tof5%H3NDU`uo<Hy1x^>^<Mb^@!i(3V
zg+|@Hp;6vz6mYTSGUL>I%NtxO+n94*ciexcGV^D@xcr<&%bs(PcW-3YSdTFyZ0EFH
zIT;6tzqIg-?RjP7(1^nD8GfDQ&=IfkOI0UA)pB!_hxcXSq@UmEK0^j*o>+c7HbNJA
zbdGG-Q8xf~6Xz=DFnwqiRbQUZqY1so8`72($%EENZn5PaqhWUU%maA~uWNX@J-(_b
zc!QzVDOAwu6v4d6&6cuxyd)K<BcEJ8^9=PRF)B#Tfs;IZ?<61KVohos*4e60K1@<~
zt$dv?Q9+{LCr|Yh?j}Nw$_Z)Rf^fHSSA~3mB4j%>tK8j7AY!`loRD?;&}7Bvw_I-w
zCo`zm+T{x&x#Urc=XNb<FMB50wpsy1v`${!w6~vZdf@69kaA1I^{m`8&%H$qA+|et
zmxN=OFF)M(=ANlYl|8v?$=QE~67kO0&{)k)o@$&|@Y!Tf>d&r=nS3RgY~hjMGrn3)
z-kIb%pH}Q8J-J$o&pZ@?q83WH@PaZdx69?z9{l~)weh}g$=8fw;q4u=AGl3H)@;F;
zWk(jn#8HKjydQNT%vU=~{;UFIbybIN80;7H-#aS1^sYwUg@H}w9HopIX61g9TK6zN
z^`8<BNSjMVxlY-aS9*qW&dv_aTg^|pCJjj5@?1sAJl-hZ?mt0vq~8=4H$EZmdtK8b
zFOP!Vi<iDv(3k-}g;opYR;xo@iHw7&n;yt$*a#=qnZU*Kvy^6^8m!YPm{y7JG6mZ_
z9r4m9Ixs@btg~8pkWab&Y;on_^Xpp6AM2|^Bk2@q-<eg;xE9uN-OzMDlXq-5$JKZh
zs^Q4dtcxEqsb=Es7_fRWkt;v?hAUz<*<Y+(nKt<>nfB7d(nhn9oPYDQSMbYdXw99F
znd>(j#_G71NOdzoGiKt38b@RB_gO1BW1Tt3*sX84zhWtzEG-=udfp78PA{}7*rp3B
za^7(%Z&QTj32#r?xN*S#d!O@fQI9mZ#=PqfBM%wsZLd@a#}TH?tX;ll+m)&3cgH<9
zUYJci_f~rqb8sq&?|FRQYm*0Q>ay*1>c}L@@`Yg@Z(k55x69_aJ!9dd#Fz_j7&35*
zf9}G)OSEB?e^U6OohHzAsq3tw!%{e`Q8Dv<nl;ocz2&b`ZwB0{$*!rpbs^){c*`_}
zIlz{4M8ZF+kLW539QAqdOru}(x_yMg6NYhK_p_|{WabUyPYP>at56v$?3|BrUZ%zg
zIv%~SQ=AxNUurS#_9ono#+4aO7YSF_W*OD~*W^fA#f0*kyb$BwT5(!L5v)YBTyuJK
zA+o;qjPYhO@EZshQ=4J~(p79_&hCS}U1{8SuX=N+zx-js^V@@b^9H^Z>zn0a!>B**
z$1d_X=A+Z*fEOA89Z_A|+nzDZ7xC#uolIix-u(5JUpIqlz5m8?>hyf-%Bxjp4eV!;
zaSl5wl9V@)oLQGkp(Kw8NeS#(e7lvLX_=cDw@v^Gy=xNV1C&5rVPBy0VgopRpj#(k
z%~FWuos!f3bQvgY+A{A_pe^WA{ELrYHiN_~+B@DC=)ffN7rFc!<iP!OL88fuPD1s*
zo0{dPr?EzBa{Lt?U522kpy4*>3C!nh2jhb6vzXcao@`%NC{lNHx3(ADjiVHdcZy{6
zy`<jNoe|j!66D5>8I>bSEJ<?J>4q1!TgbMVSA6&NCJ?n&P4Nw_*<?FcZt(X(9`mtW
zo3j8%IeES;O!J;;6=~(nUob_yhTIno_?jG4L+Y(2^(H*7CXH3QE`_xe63&()0S`T6
z4gWpD+C3qLjLGvVG#3|7WcnL&ID{X*$P~(-<Y1SmOl>oiyFEuCfs&mhdc9=rJIY~W
zzIOT+X<|NQk>_cbWrWu+k6m+P2$|*d6h8MRlS4WDN|t+aN#ArO6Vu{SQV>=!Vz%x>
zBG;91Tg>bUIkhsfGqtXcXvXkAY&U;KEOJBNeq`2?+HsR-c<vhf4khVujoCd@jnSk1
zgunBcFbwTQGIX1!Fi%EBrEluE%ygb0m3X^zJ{7IznO?XtnbLaw=w|4#c8UmapGh36
zNc^1N`a6RIQQCa|^$YDCM0s(vXtcvg^6C8O{IjdC5$-Wj&nr{!5!EUV*`otB#DVAV
zyo$sIa;2>;xiO-Nto-`v{*me?axFu1acgKJ(b@%(DGv9D@p}VtUe~1>XHKqWdaGG5
zIvNv`+{{Fok0g>TVk`=na-#G2zFTTi&jP~KCf-P)azxrxpM3mG88vY$O`oSs1~(Yh
zVY$vkLbFLW%pjULsdGpKy*o?#&V|W;Z!9DuuANM^=6Xmt-?ncv|MZNo)rv1#xS*N%
zcJ_7M<9|(_cT;oay<Zb+D|e<-*-LU^*Ugcp^_67BD^+$K6*~>t__-Ndd~6uC_Tk4C
zUYg0A%`bL!Af|-5I%Qq`y*;{==jqmqD$~zW70TBG_bvNIy)dfhZFiqf1bRa|_F1ka
zkCyIyv)^MMIW40;z~7NYE}t8bbt$-%yqvYJ#rt9nNq-p6Hhb{=JyVZ^LsYwkILf-7
zu)p)3jF7mye2Z2q@!pl%G)?XuDY^bR4Yt$}ZZ_+}?S&2+M=unIYI3Y(C{HLJKl`E7
z;Lm>qqq3XsFjsGE*xk`=NIlIGG(YE=MLD127B-9@LG<U|T9h_YlZ@==c)3T}i%h9i
zqsA%6lha<A18L?tBu908`e%!K<cw{@!E15P$oUN%E~+12ld6NADmpvcNQvp<L$kJY
zkSE!)GIwuxkm*<WFUtkCleJIx-#fVE8QCEDBHGW^MdP}`$bLy#Cq~gjHRTmBhiM=q
z_T|vR2h4}3wE{*Nn^2}lG)1h>U!*pz&JTL`f}N-}o;bFBwl1k~+7x{2v>&-#_et-N
zN+OxKoo!<NlPjcBMqab8x{|b?zacWAzln%Z&!+h|w36DSZS%&ccamME?@iL!x`=9$
z^yY5IE>iNaYu<@#pGkJwW8aF(Cepc4ykLgJIt>v4)eJT#H-=?&_Zq`<^O%-rbH4b6
zR5L}Sa#yHsx1g*S*u{$TT%n{L%g-~6xXFQum2bM28<GV{=O#|l-9!eaa`$g`JWkf-
z3mrN0YOs!H#ciSDT|-(56rb!$c}?;Y0wZ)}KamG+p}FrByNTh9nt_6xUJ}<U(D_il
zk7TV|TeWWLHxeGQa%<$^d3$`~;*`s|UK)3eRCPyvUC%h|r1g-aSe-d1mb1L6@+q@7
zct^y<Y8$G6n{#KlO96G_-m&}bQM|<M`4KC=k4s42&1ioY>0lCeAyF+%B9&}jsCDa;
z?oHBZ)a><0^BM8{zFmbYsFj@KT=9KXLl;S^o4%~fdw?8W!L!iDbMXJiW?u(v*G7VP
zMfK5q_YsgkWzEBUW-D>#Jd|E<?5}a}$%a>r(|sAzH<o=}A_vUC_y&%XnvKkUIj)ca
zt(8>n+U1wFpDdvgW}JIa@koFW#U#0-LKb8|JJN4ULKvwHo}U>fo<8_JSCByb%`&pW
zdUDK@TTSE}m;c&;j~#@a@5Lm$+kK?}>xyzMX*Q4zEp4v-!~tZ*)JN)Oqv4xxf5xqJ
zc2K*d{`Tz64icLyJ%?lBW)01)f;CYVfega?*hN1{hxzpGjsp#TFPQ?Xw~It=aH1ZJ
zDCvCJe}_8d2hsaArV^_SISE0FY{==cGfXed-9_AwpZXAFd7iwOzqoGZ(+A}7l5dA!
ztb9#w8*FRYyRVDn?()1FKW!9rJd&thB*F>qycK)Alg7fK-1STI1jj+|mU}uqqq*U@
zo`L)$`z~@!Pv%X>xnPa4aaqiW=535n&fb{%3kJ;fai=DYD0<5*)Jeb2Q{qO+N%g&C
z*gc@u%uh|T2@xabcf6d`^K=C{Y#KQ|N@Op&8Eh}Jf8r(bN$c80U*$(cWp>(?N7AjN
zUN|;;jB6kH4(kpVl#PZr?ZUd-ta;$som(5+9r$49sqh~CbUwIiQ783oB`=gFhb}Xc
z>LdPAr>gWJTx0SjPS1G$?TnI-I|JqlEn!L@_vSsM+Rhw3@&3*_HBahWv17Z}!D{MI
z!Yot6+Y%(gBS?y8oD=!LJ9FaXUHi$cmqj}*#^w_)@iKvo%v$1Y7*pVu*FnBdR`Xx*
zd?eUy7Wr~DeH?6fYGS%pa3WZ)+SHQTFbOImU-UbX$-u{{5t4aC0B-MBTdA2h5=O``
zlx-LprJ;BIK-riDI~gac+4JVbnlo#pXMK`)`^*eFK~|n$<4YCoNd3tC@Ra&yYC1_z
zT86N@unCM}x{*XV?a2p44wEW7(Xlu6ipZ|anab{EjifF)XJ3y|7vVLK9O<)pH0+Ye
zoo>&~4`(!Hc2%=a0n-yJ)9!={L)k~~FFAd}u;j_vS-l^HpfI%|ValS>ka9kArqr>$
z8pXxS4}PlH%}9+ocQ+x$nyI}@i}4}5hske}w(ZJ=K<d`|84ZU`8!7gLoKB4mie&#`
z&1Y2cdLpfH`LRClF_QhNMf}P4QlfQbX||u(D>BZ&##?lJAIUlY;)AE*I9Ra7r_?KZ
zGB_XoJR?|jI%K?*e?S(C17#=t;B108OkQ!ySkF`p><HY8<sAo6@Au6+GxvZ-tj?-O
zjL~rncqe@K1;+|zsJ7{1H@Oj1&k;YNK=~j_j5Ghx<H(nkVcXH!9ha3!bBbz%O1}@8
z6+uP%2c9GkZ)AohwBI9oqZT#KjeAFYPxA}<q>LQ=QA+8i`a68UxL6b^A2|&af?jY0
zf1d$z19ofkif4geVK0xH#BA_ZHJ>T7LjtPWmH(dx?)x3<FN^~?WhZ+^*&+$q&-3}5
z+b-iVLPfGkiHx$-BQlc69)$>rN|9BzgeV~)zRF7WCXuh7&OdNp=Umsh->)}=)>Bq!
z`fG5#Y%zn7qRa80t6~~4T+?9uPq+j2()NHGZbgN!98p{KxEO{<slIyM`u#U<o>jE{
zn2Uh42>!@T!XJHkM=JkKU4*_1H!=7sbR&m)yH|!4i^%RjBhwKX8j$2#N#VBU0)eka
z7~_%<L=72j)sG$p`Q?A@z89rHx{1KF@l6WO8x1c0`63C|>P-yPmAN1w*(pFj`7wcl
zI2SATE{kZ;!}EzR$pverqCe&Sm>!q-d9S8ACklW4i#5eeXA!>|66(O{s)km_j%cbB
z-a>(9Hv5$*ULlpt!Jfd$KEysY`kPF>hK5`yQW+)=Krn-0a^GQo2v*&#S-B?;n_W8$
zrqm}O3)PS1o|J<({BQg$$$O_%!B0+17NQOKC@(GZfsF|3(vnC%AxHOjbJd?*qH~^e
zqvEnV7PQ2oI3~n`k4DO|#0baX>$_rwazU%O0DFV)XrUIeTaaL_nh8a>22MNLI=n%3
z)u)noI);#n;pF0Lw=Hz4=6FST1S{M~?@c@QOb9k6FAgu69S3P!!n@S(C*g=ydbfXu
zBGlcOzNGd~5ia~7j}5pefM0CweX#@~XkTd!zA#cmAetHV#)%aYcQl8aJOaEisY1oB
zIztXzaAiL2j8P&!@Lt{}tMv~acC4vfk5Uh@Re#ie&Ju;>{?*aT{i;S`k%5~ce4~h2
z<pL=nW6#N&&zGg_{WptO-M$knBMv(S@0DBH<REPLSbqI;Ww<kcu}D~$2uk_hdJNe_
z=(9W>%$=_a2WJG28}o_-qZ41FVe~5kXOYt?!_${UjMKz~_l_U77_>ge8^?>k7$<kS
z#HHa@>D>FuRVk47{nIw(R>sK2f0HdCG8U=$o9o6**P|GDeeZpZdrsS-VeZoq6@)AM
zoJ-5+gS&@jxgYmS!IbZLoxn3n@Nk`MK`O_==nrB1CXG6@og(+no>hmN3z*qFpBhlz
z)nB~)K?+EN`V{-Q$_dRZHwPzXONkE)=_ZPYgE6aHBULZC1o1g(CW{!ZY@Ewjn#t9b
z3b8VJ+FVUGL*tvJ=Kp#VP@h!4iSOP!%}-M*J&-bst{Z4ACdx4YTT`ss=6NBgI8CEu
z?jr|_G0JcGcZiS>w)IRkTLb>tijb5owBhYH32oDQZP0Htj4)8ug2{Ep$Qf@r5DkzW
zf71DuFu&<lD$Q9%9AWzJCgv83HB%UQ$f}9r+Q&k&Zo?BitF=pw^EExHU`kbaFlm8K
zSf207TTDfIDo#P+)$OP~eyU~P>Js8(N!H1AXMron^ULz<#GvO#T{Zf-=bkNj)-Oh=
zLvvTJ$*k!qaD?wWl(l*g;rH=LrG`GVmZ=(W+vq}c31_9(j3Qha7V4lSHxSnCr`_oi
z>WG4?eTH<`BeD3v$W;29k~ocm{X?6;Lj1)#8>{m%7L;(UdEl^~E$V#q_3MCT7OIL!
zmRM@ryCMTA3wt8g(dz~uHv2RV(2@w}`hZEn+m~X}{dPn^m$JD6=d>Zci@N`XGznO^
zgo%#T#_*}jW?6;B1iEu?)A4i~LA}(puv6&Xyz!rVn?~D0nCg7*W;@bE9AcTi8u;Np
zc15Z(=MaM|9+L0VaD}b}e{bx0#NjLl%4R#gqa5Uj`ou%?{*>mTmduR5h8$nf@%HV=
zOto#~+GeVD?;|gWcE&J<(8<Bcq3MxkK6Us4v$q%R^?|9;=kb*C84z8W>*kp-gN^mA
z2XD;G!M5&JWbl+JP_0i5(x|J$5$#NF>9-#V&)vU9l<aRK3N*}}D6>w$2ze8>9x@7e
zFU7kA((y9f+u|DgmqcDvG_1IGwBa&RuYZlK>z^Vo(_=AhUV|vrC@xp#H6@hoY*(@i
z34%-2_*S2mA{gx49cVeL4U+j(g_l@Oz({o>P^0WDOw?C!Ix<>97v#HK$+QH~%5#pV
z#?OP#H?9Ybnc5)Mx#=Z-{R`ooFFL9l^qIK%=Rr-MZ7OEqDiU?eSrsqKyuC4JTZK!i
z`WLov3L-^cS}J}S5416TRW+IYB?>MkS+OPlKnJgHUs@NTfnYoC0d5gdxbY%%X_ZWX
z@pGy6Ju!MPdQRqWtBN`3?zm0xsaiq0+(`xZleSRn!d;OXZwq-m(P?9>Hn5b)rl`cM
z4{90<=4kUPVa~^I+qbKSc=PkY`N`**n4V3$WaJv)-Lz3n|7Pm&4%%?jSrt+AXz?8V
zd$a4vx=rWH&)`z@S`%kkZ5&5MSg>;X(*uC$W`1>WNkU%x*pE&>HIQx5yzcLB3?Dc3
zH%=K^0;!<WwJOLCs&p2jm^dBb{w|$bPNE}-4Hk*!SvbI}Gjz|3ri@`L7uP&d@STuK
zI!Hya*iVeAqi_l!<zf*m=Mqg0Y2gi+vyw4eGp^{JvS0m*B;sp&b9EuyAJO_$cf9tf
zL>=uVUKcigqv}R!V@<>gL7&L_tI{&S|AT*8=CT$D9AWcmwBGY3y~3f~$QE7*SNLB_
zb^;GtaiPvnF7Q#PmPW__GSsQEcv>*H0O|3A4=g<9p#I;JHu_6H33i3+L|Xe1Vil=^
z$x-(yHnbaO$fT=_A2SJR_)KiWm%h*kU6_+WYZB&9>h_b-;Qfb5q>NfbiZdybS6W2!
zgAGz5-kgwY#ye(fEDt74(l-T9>h9%3C*DYMSioOhk4E(k2Vhp~xaSgn1?+{b!sbKV
zVP7+omd&MWu(@b0DbV5uF+uuAtJN)mYxvZs-;)yrRhvPE54B@Nz26@S;$vQ519BDn
zXFnR@wmPRp^~*owSIZ_S_8m|_4sj*zW$a;SG`=aZLZk^<->sacQCLNRd5KPhj6;wZ
zQDyAZrwsF>cgOFo8A9qS5v!GQYbbbR@@Dgr3y6!S&B!Tu0Qsj9N#Lm$U^x3TKTU7Q
z?dyscu=aw9>eeHi2DV`Lrvzt@n<1RxRXsfMZHmZ#aNx_dQYmJVwICC_ZHm948azic
z(}UYnMwRr4sUkOJmHY2RqmW+_=T(u(52%=WM!4_rAM|$myG)0tAP7(^q#5htur!`}
zTX*^lF!n7Njg{C#-K^tAM2;Ic`&|ApTj&j3zQNvQMLrN$a=TdUmM;`9&5};W`oLT1
zJ~!&Uy!N`-q_F?FB|^Wn@3U651>%pok1>9dl^8wiy3X=p3%ri&B-Bt2;I~_AZaxkN
zq&p~@pYIimQjTLLiLIT8Dv^+&BD{;7O1(79ktm$u(aVb)R|g)s7tP}9XF)3Yn~Gzi
zGYFL(EqOHH3E3nASJq-*(BC~<PCgq5!<xN*Jib9-!o2xUC^it*ytr5b2V6kEwK?O(
zmo)-O@ws+F!3r^Wj@j&|O)d7reV4E6o(=xQ>Sk|=%?PgC*#E(^MibdDQN(5CB%qg{
z6KAqRyV3PxdggG#eh?(=-+!S)0-E<cZg7=0tZ_XuU3+=~iZz%eeloZMrzq#4K*K#x
zQSaNlL=1v0_onu9XUVYoaJSTsFJ#Zbjwao#CPVhZk+q8f?!Yn@som=Jm#|V<&COJ?
zNi3RS_V46u!V+q87Y92X@C1QcpS8j<+}M#$@<hpLv`}`ruiP&cCHROo#IN+BS(cC+
zDh0G)wp_O5qi_Nu!b6%wW%R%^zMJ=Ds4cLEUOr4i?*WKP^s$PRKddiWOnDWO!Ibj%
zkAG&Nko70b0k;bSm(vCVtCC@`P)l=XCdmr~U)ekT)1V}BWWfFk)?MOk0c*g*<_Byd
z;UXoW-vz(;UGsxu))f9xiMhMF#}F}N?>iLSGLeSeP=BT65ZdL-nvY*)1erIo0Z(;K
z!i(%^uY)8bNT`z(^onqR9H8j`@3R+NVBBJ=b-W1*rj5TAJVN0aqn>O`b~prRHixJk
zh=f<oEjKdxBf&no$Ny8958N8#Z!0aMCjPQtJ$v=ee(Y-$^VUD^PnfXD$+cCHYq;j^
z{sWiV=kZd8m-lf+Q#5(Q;E?xF4hj?Epnu@`6SW)IL^Ecyg3)_PDbcmPzK)fXA_CE8
zK=zS!p#FDfAdj<I_3Y)y>B4jMD%Qy0n3Aq@@XDUAd8*PUDsvZh`dwHKaomG%al5mz
zN%uhHXWdzhBM?}<G+YR}2Z;4)E;H+4blAdbxEt?}ZY)Q|yEiJ>8<%~QWNW0dg2#9c
z_AmZ8hZxkxi}OeYi07}DWXfK@3@$`|DPhV9rcotdvZV>Y^X9U)>3ef{RYpGl!2b%Q
zONhU4rSXT@U%B(l=R?6D#K&&oZWQ=-ab&P=#DJZ%kWe8Z4${Uqyp`JH;2@G@&lDv?
zHx+f^fG8`mRMEmW_yZHB%w}1@`}7+&dt*wt!p0A`yQ9l{Y<UB}J?uBhm39#+ZTwCv
zEiOiMUeBp_siu+B+4a?}Y+hivv+oky27*!B@ViCk7U1tA5*Kgd4w`Ew&ZC!uAb<4J
z%7wl=Fjcwh^YQaND5$6WEVUH}8`~qto+QP?@I$OhG&%v=yG%IySVG~1pVAG#NKT@0
z={v_&5<9j^-QSk|Xc((Cys=|?H5gav5H|~G+`(fHimbVuvPX+agIN>zU!i+7kC?j2
z3y90pDuE;U2v}Zl%I!L?0Rei~7TVM<?(Kd0+*@=#VMK<~YmAZ%Nwe|Uj=qtg_4J=f
zyI>sLVA||1<4yn<dWvc0!%1-WccfQNbTYUvD9zgGhC>udvbB(wk4U`9Rp}|mgVB`V
zisuX+#j10Qxq_E(<E=`fF6Gzuq3AwwLH|H!l-0M^lVbk{&GpGWuvuP4lbQ5k{sO}A
za(p#zPF5S5*AtR@+iXF$y!qXgByU&|59;z+4*}b$R!?Q)XsC6_&dJ+(0Pu*-h2Jg-
zoSCEt-}R?JTh+4Pc3T?svPZ9<XpVx5hb&xrR|JWCy}72_r;lK}ZDR7IL%*=jaJ*LP
zZUk<SzmM{2BMqXXY)sz>bwv+nsZ3Hvt5AZN`j=Oa*U^FN1x^`eVvzawz-Q|C(~w{D
zhU_Ho0LpZ_RMFx-pmR)m^}R$GB$l)K(B8QZLe28yvYd$^Ny({kcs2#91Pt#>_N7DD
zP-t?g)+6A}SZzu3iUE1KA)j0%Mr4va6jU=OjH&aE$jZK-#rnTn_&Ohr!8P`^D2x6%
zfD&@`(o1SQP_vUH>C53dBp7{2*12yB$)<6#m@G@cQ{Af3_tpk5*4<G)J?{h-Or&3_
zJpPbI+s{N@9uAy@p<FuQ2S6%*a<Dx#88~>P7_EEK!H?vw8zq(rbNl%9N`kT>Qa*u|
z;c*;X7aYuONjXN0tBv|rl6e%{UDa^k8eYPfuSjCV?FTsh*;Q)FL#!x5{w-xbb{%!z
z7<283ZA4A0zPG+7@1nH4BgHzD(x6f!)M;sE49waMJ^3t`fleSRcP=sr?sAnp8+jWA
z8SIt8iWd{WP0c7?#W4-${4<@nWip{#qP>@CIvcb#V}BWX=D?Bgfytcv@jzRlJ6Vw_
zOH4J7oLY5}#=6dQJXk!kj;%O8C;sC~#&^a;UM011Afucdv!zphX!pMhjQGZT^p1nP
z8|1wY@TmsXyk^<GUV{1ArMfd<yW$i3wD2m#RW5}vog%~cmS9gEqZlCa2f1$UO9r-s
zP3hlT9zx`44<|92#}I1qdt-b#2OR3@N{0fTKww%wN9Q{QqLx}<QM<-TEPI9W`Sz*5
zn8y`*PT}ozTx3=-b%cW#xhD%1-Wm@=l7$-?f12A+O%1yYa-fELQNJzcbmSrUvSJyd
zk2##DBm8u=_1N=)&SkbKw_(~~y;ZC_4jLbZyz}8mgZ*<Q)>h`(keNXqJ}RFJ<$}{&
Q>JoWy{u{Yo>ufUo4@^-v{{R30
literal 0
HcmV?d00001
diff --git a/Generators/ForeseeGenerator/python/ForeseeSampler.py b/Generators/ForeseeGenerator/python/ForeseeSampler.py
new file mode 100644
index 00000000..60df119b
--- /dev/null
+++ b/Generators/ForeseeGenerator/python/ForeseeSampler.py
@@ -0,0 +1,433 @@
+import numpy as np
+
+import ParticleGun as PG
+
+class ForeseeNumpySampler(PG.MomSampler):
+ """ Sample from the output of Foresee generation in numpy format with columns E, theta, weight"""
+ def __init__(self, model_path = ".", model_name = "DarkPhoton", com_energy = "14", mother_mass = 0.01, coupling = None, daughter1_pid = 11, daughter2_pid = -11):
+
+ self.path = model_path
+ self.modelname = model_name
+ self.energy = com_energy
+ self._mass = mother_mass
+ self.coupling = coupling
+ self.daughter1_pid = daughter1_pid
+ self.daughter2_pid = daughter2_pid
+
+ self.rng = np.random.default_rng()
+ self.xs = 0
+
+ self.read()
+
+ def read(self):
+ "Read data from numpy file in format E, theta, weight"
+
+ if self.path.endswith(".npy"):
+ filename = self.path
+ elif self.coupling is None:
+ filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self._mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+ else:
+ filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self._mass}GeV_c{self.coupling}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+
+ print(f"Reading data from file: {filename}")
+ self.data = np.load(filename)
+
+ # Create probablity for each mode as weight / sum(weights)
+ self.prob = self.data[2]/np.sum(self.data[2])
+
+ return
+
+ def mass(self):
+ "Mass converted from GeV to MeV"
+ return self._mass * 1000
+
+ def shoot(self):
+ "Choose a random item from the data, based on the probability"
+
+ # TODO: what about reuse of events?
+ energy, theta, weight = self.rng.choice(self.data, axis = 1, p = self.prob)
+
+ self.xs += weight
+
+ # Create mother momentum, randomly sampling phi
+ mother_mom = PG.EThetaMPhiSampler(energy, theta, self.mass(), [0, 2*np.pi])
+
+ return mother_mom.shoot()
+
+class ForeseeSampler(PG.MomSampler):
+ """Create events from foresee directly on the fly
+
+ Requires:
+ * foresee to be downloaded and in python path
+
+ cd <PATH>
+ git clone https://github.com/KlingFelix/FORESEE.git
+ export PYTHONPATH=$PYTHONPATH:<PATH>/FORESEE/src/
+
+ * scikit-hep installed
+
+ pip install scikit-hep --user
+
+ * forsee files dir symlinked to the run dir
+
+ ln -s <PATH>/FORESEE/files .
+
+ """
+
+
+ def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, mother_pid = None):
+ self.modelname = modelname
+ self.model = Model(self.modelname)
+ self.energy = energy
+ self._mass = mass
+ self.couplings = [couplings] if isinstance(couplings, str) else couplings
+ self.mother_pid = mother_pid
+ self.daughter1_pid = daughter1_pid
+ self.daughter2_pid = daughter2_pid
+
+ self.rng = np.random.default_rng()
+ self.xs = 0
+
+ if not os.path.exists("files"):
+ os.symlink(os.path.expandvars("$Calypso_DIR/../calypso/Generators/foresee/files"), "files")
+
+ self.pid_map = {
+ (11, 11) : "e_e",
+ (13, 13) : "mu_mu",
+ (22, 22) : "gamma_gamma",
+ }
+
+ self.mode = self.pid_map.get((self.daughter1_pid, self.daughter2_pid), None)
+ if self.mode is None:
+ sys.exit(f"Undefined decay to {self.daughter1_pid} + {self.daughter2_pid} for {self.modelname}")
+
+ self.foresee = Foresee()
+ self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1", channels=[self.mode], distance=480, length=1.5 , luminosity=150)
+
+
+ if self.modelname == "DarkPhoton":
+ self.data = self.darkphoton()
+ elif self.modelname == "ALP-W":
+ self.data = self.alps()
+ else:
+ sys.exit(f"Unknown model {self.modelname}")
+
+ return
+
+
+ def mass(self):
+ return self._mass * 1000
+
+ def darkphoton(self):
+
+ # Production modes
+ self.model.add_production_2bodydecay(
+ pid0 = "111",
+ pid1 = "22",
+ br = "2.*0.99 * coupling**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10)
+
+ self.model.add_production_2bodydecay(
+ pid0 = "221",
+ pid1 = "22",
+ br = "2.*0.39 * coupling**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10)
+
+ # Handwavey
+ self.model.add_production_mixing(
+ pid = "113",
+ mixing = "coupling * 0.3/5. * 0.77545**2/abs(mass**2-0.77545**2+0.77545*0.147*1j)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ )
+
+ # Question on validity as FASER gets larger
+ self.model.add_production_direct(
+ label = "Brem",
+ energy = self.energy,
+ condition = "p.pt<1",
+ coupling_ref=1,
+ )
+
+ self.model.add_production_direct(
+ label = "DY",
+ energy = self.energy,
+ coupling_ref=1,
+ massrange=[1.5, 10.]
+ )
+
+ return self.decays()
+
+
+ def alps(self):
+
+ self.model.add_production_2bodydecay(
+ pid0 = "5",
+ pid1 = "321",
+ br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
+ generator = "Pythia8",
+ energy = self.energy,
+ nsample = 20, # Vary over out of phi and theta
+ )
+
+
+ self.model.add_production_2bodydecay(
+ pid0 = "-5",
+ pid1 = "321",
+ br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
+ generator = "Pythia8",
+ energy = self.energy,
+ nsample = 20,
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "130",
+ pid1 = "111",
+ br = "4.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10,
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "321",
+ pid1 = "211",
+ br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10,
+ )
+
+ return self.decays()
+
+
+ def decays(self):
+ # Decays
+ self.model.set_ctau_1d(
+ filename=f"files/models/{self.modelname}/ctau.txt",
+ coupling_ref=1
+ )
+
+ # TODO take into account BR
+ self.model.set_br_1d(
+ modes = [self.mode],
+ filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"]
+ )
+
+ # LLP spectrum
+ self.foresee.set_model(model=self.model)
+
+ plt = self.foresee.get_llp_spectrum(self._mass, coupling=1, do_plot=True) # This is just a reference coupling
+ plt.savefig(f"{self.modelname}.png")
+
+ def flatten(l):
+ return [i for sublist in l for i in sublist]
+
+ coups, ctaus, nsigs, energies, weights, thetas = self.foresee.get_events(mass=self._mass, energy=self.energy, couplings=self.couplings)
+
+ return [flatten(thetas), flatten(energies), flatten(weights)]
+
+ def shoot(self):
+ # Create probablity for each mode as weight / sum(weights)
+ prob = self.data[2]/np.sum(self.data[2])
+
+ # Choose a random item from the data, base on the probability
+ # TODO: what about reuse of events?
+ theta_mother, e_mother, w = self.rng.choice(self.data, axis = 1, p = prob)
+
+ self.xs += w
+
+ # Create other momentum
+ mother_mom = PG.EThetaMPhiSampler(e_mother*1000, theta_mother, self.mass(), [0,2*np.pi])
+
+ return mother_mom.shoot()
+
+if __name__ == "__main__":
+
+ # Testing ...
+
+ import os
+ from math import sqrt, log10
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from DIFGenerator import DIFSampler
+
+
+ path = os.path.expandvars("$Calypso_DIR/../calypso/Generators/ForeseeGenerator/data/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy")
+
+
+ modelname = "DarkPhoton"
+ mass = 0.1
+
+ theta = []
+ mom = []
+
+ d0theta = []
+ d0mom = []
+
+ d1theta = []
+ d1mom = []
+
+ # Accounting for rounding
+ epsilon = 5
+
+ # Mother position sampler over Z
+ mother_pos = PG.PosSampler(x=0,y=0,z=[-1500,0],t=0)
+
+ # Create mother momentum sampler reading data from foresee
+ mother_mom_sampler = ForeseeNumpySampler(model_path = path, model_name = modelname, com_energy = "14", mother_mass = 0.1, coupling = 1e-5, daughter1_pid = 11, daughter2_pid = -11)
+
+ # Create decay-in-flight
+ d = DIFSampler(11, -11, None)
+ d.mother_sampler = PG.ParticleSampler(pid=None,mom=mother_mom_sampler,n=1,pos=mother_pos)
+ d.mother_sampler.mass_override = False
+
+ # Loop over a range of events
+ for i in range(1000):
+
+ # Shoot the decay in flight
+ daughters = d.shoot()
+
+ # Get mother and sum of daugthers and check these make sense.
+ mother_mom = d.mother_mom
+ s = daughters[0].mom+daughters[1].mom
+
+ try:
+ assert mother_mom.E() - epsilon <= s.E() <= mother_mom.E() + epsilon
+ assert mother_mom.P() - epsilon <= s.P()<= mother_mom.P() + epsilon
+ assert mother_mom.Px() - epsilon <= s.Px() <= mother_mom.Px() + epsilon
+ assert mother_mom.Py() - epsilon <= s.Py() <= mother_mom.Py() + epsilon
+ assert mother_mom.Pz() - epsilon <= s.Pz() <= mother_mom.Pz() + epsilon
+ assert daughters[0].pos.X() == daughters[1].pos.X() == d.mother_pos.X()
+ assert daughters[0].pos.Y() == daughters[1].pos.Y() == d.mother_pos.Y()
+ assert daughters[0].pos.Z() == daughters[1].pos.Z() == d.mother_pos.Z()
+ except AssertionError:
+ print("Error on run " + str(i))
+
+ print("mother particle:")
+ print(" E = " + str(mother_mom.E()))
+ print(" M = " + str(mother_mom.M()))
+ print(" P = " + str(mother_mom.P()))
+ print(" Px = " + str(mother_mom.Px()))
+ print(" Py = " + str(mother_mom.Py()))
+ print(" Pz = " + str(mother_mom.Pz()))
+ print(" theta = " + str(mother_mom.Theta()))
+ print(" phi = " + str(mother_mom.Phi()))
+ print(" x = " + str(d.mother_pos.X()))
+ print(" y = " + str(d.mother_pos.Y()))
+ print(" z = " + str(d.mother_pos.Z()))
+
+ print("daughter 0 particle:")
+ print(" E = " + str(daughters[0].mom.E()))
+ print(" M = " + str(daughters[0].mom.M()))
+ print(" P = " + str(daughters[0].mom.P()))
+ print(" Px = " + str(daughters[0].mom.Px()))
+ print(" Py = " + str(daughters[0].mom.Py()))
+ print(" Pz = " + str(daughters[0].mom.Pz()))
+ print(" theta = " + str(daughters[0].mom.Theta()))
+ print(" phi = " + str(daughters[0].mom.Phi()))
+ print(" x = " + str(daughters[0].pos.X()))
+ print(" y = " + str(daughters[0].pos.Y()))
+ print(" z = " + str(daughters[0].pos.Z()))
+
+ print("daughter 1 particle:")
+ print(" E = " + str(daughters[1].mom.E()))
+ print(" M = " + str(daughters[1].mom.M()))
+ print(" P = " + str(daughters[1].mom.P()))
+ print(" Px = " + str(daughters[1].mom.Px()))
+ print(" Py = " + str(daughters[1].mom.Py()))
+ print(" Pz = " + str(daughters[1].mom.Pz()))
+ print(" theta = " + str(daughters[1].mom.Theta()))
+ print(" phi = " + str(daughters[1].mom.Phi()))
+ print(" x = " + str(daughters[1].pos.X()))
+ print(" y = " + str(daughters[1].pos.Y()))
+ print(" z = " + str(daughters[1].pos.Z()))
+
+ raise
+
+ # Store mother info to plot
+ theta.append(log10(mother_mom.Theta()))
+ mom.append(log10(mother_mom.P()/1000.))
+
+ # Store mother info to plot
+ d0theta.append(log10(daughters[0].mom.Theta()))
+ d0mom.append(log10(daughters[0].mom.P()/1000.))
+ d1theta.append(log10(daughters[1].mom.Theta()))
+ d1mom.append(log10(daughters[1].mom.P()/1000.))
+
+
+
+ # Plot mother from sampling events
+ prange=[[-6, 0, 120],[ 0, 5, 50]]
+ tmin, tmax, tnum = prange[0]
+ pmin, pmax, pnum = prange[1]
+ t_edges = np.logspace(tmin, tmax, num=tnum+1)
+ p_edges = np.logspace(pmin, pmax, num=pnum+1)
+
+ ticks = np.array([[np.linspace(10**(j),10**(j+1),9)] for j in range(-7,6)]).flatten()
+ ticks = [np.log10(x) for x in ticks]
+ ticklabels = np.array([[r"$10^{"+str(j)+"}$","","","","","","","",""] for j in range(-7,6)]).flatten()
+ matplotlib.rcParams.update({'font.size': 15})
+
+
+ fig = plt.figure(figsize=(8,5.5))
+ ax = plt.subplot(1,1,1)
+ h=ax.hist2d(x=theta,y=mom,
+ bins=[tnum,pnum],range=[[tmin,tmax],[pmin,pmax]],
+ norm=matplotlib.colors.LogNorm(), cmap="hsv",
+ )
+ fig.colorbar(h[3], ax=ax)
+ ax.set_xlabel(r"angle wrt. beam axis $\theta$ [rad]")
+ ax.set_ylabel(r"momentum $p$ [GeV]")
+ ax.set_xticks(ticks)
+ ax.set_xticklabels(ticklabels)
+ ax.set_yticks(ticks)
+ ax.set_yticklabels(ticklabels)
+ ax.set_xlim(tmin, tmax)
+ ax.set_ylim(pmin, pmax)
+ plt.savefig(f"{modelname}_PG_m{mass}.png")
+
+ fig = plt.figure(figsize=(8,5.5))
+ ax = plt.subplot(1,1,1)
+ h=ax.hist2d(x=d0theta,y=d0mom,
+ bins=[tnum,pnum],range=[[tmin,tmax],[pmin,pmax]],
+ norm=matplotlib.colors.LogNorm(), cmap="hsv",
+ )
+ fig.colorbar(h[3], ax=ax)
+ ax.set_xlabel(r"angle wrt. beam axis $\theta$ [rad]")
+ ax.set_ylabel(r"momentum $p$ [GeV]")
+ ax.set_xticks(ticks)
+ ax.set_xticklabels(ticklabels)
+ ax.set_yticks(ticks)
+ ax.set_yticklabels(ticklabels)
+ ax.set_xlim(tmin, tmax)
+ ax.set_ylim(pmin, pmax)
+ plt.savefig(f"{modelname}_PG_d0_m{mass}.png")
+
+ fig = plt.figure(figsize=(8,5.5))
+ ax = plt.subplot(1,1,1)
+ h=ax.hist2d(x=d1theta,y=d1mom,
+ bins=[tnum,pnum],range=[[tmin,tmax],[pmin,pmax]],
+ norm=matplotlib.colors.LogNorm(), cmap="hsv",
+ )
+ fig.colorbar(h[3], ax=ax)
+ ax.set_xlabel(r"angle wrt. beam axis $\theta$ [rad]")
+ ax.set_ylabel(r"momentum $p$ [GeV]")
+ ax.set_xticks(ticks)
+ ax.set_xticklabels(ticklabels)
+ ax.set_yticks(ticks)
+ ax.set_yticklabels(ticklabels)
+ ax.set_xlim(tmin, tmax)
+ ax.set_ylim(pmin, pmax)
+ plt.savefig(f"{modelname}_PG_d1_m{mass}.png")
+
+
+ print (f"x-sect = {mother_mom_sampler.xs} pb")
+
+
+
+
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
new file mode 100644
index 00000000..20f3893e
--- /dev/null
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -0,0 +1,198 @@
+from AthenaPython.PyAthena import StatusCode, McEventCollection, HepMC, CLHEP
+from GeneratorModules.EvgenAnalysisAlg import EvgenAnalysisAlg
+
+import ROOT as R
+import numpy as np
+import os
+
+def fix():
+ "Python Fixes for HepMC"
+ def add(self, other):
+ self.set(self.x() + other.x(), self.y() + other.y(),
+ self.z() + other.z(), self.t() + other.t())
+ return self
+
+ HepMC.FourVector.__iadd__ = add
+ del add
+
+ return
+
+class HistSvc(object):
+ "Class to deal with histograms"
+
+ def __init__(self):
+ self.hists = {}
+
+ def add(self, name, nbinsX = None, loX = None, hiX = None, nbinsY = None, loY = None, hiY = None, title = None, arrayX = None, arrayY = None):
+ hname = os.path.basename(name)
+
+ if title is None: title = hname
+
+ if nbinsY is not None:
+ self.hists[name] = R.TH2F(hname, title, nbinsX, loX, hiX, nbinsY, loY, hiY)
+ elif arrayX is not None and arrayY is not None:
+ self.hists[name] = R.TH2F(hname, title, len(arrayX) - 1, arrayX, len(arrayY) - 1, arrayY)
+ elif arrayX is not None and arrayY is None and nbinsY is not None:
+ self.hists[name] = R.TH2F(hname, title, len(arrayX) - 1, arrayX, nbinsY, loY, hiY)
+ elif arrayX is None and arrayY is not None:
+ self.hists[name] = R.TH2F(hname, title, nbinsX, loX, hiX, len(arrayY) - 1, arrayY)
+ elif arrayX is not None:
+ self.hists[name] = R.TH1F(hname, title, len(arrayX) - 1, arrayX)
+ else:
+ self.hists[name] = R.TH1F(hname, title, nbinsX, loX, hiX)
+
+ def __getitem__(self, name):
+ return self.hists[name]
+
+ def write(self, name):
+
+ f = R.TFile.Open(name, "RECREATE")
+
+ for n, h in self.hists.items():
+ path = os.path.dirname(n)
+ if path and not f.GetDirectory(path):
+ f.mkdir(path)
+
+ f.cd(path)
+ h.Write()
+
+ f.Close()
+
+ return
+
+class EvgenValidation(EvgenAnalysisAlg):
+ "Gen-level validation"
+
+ def __init__(self, name = "EvgenValidation"):
+ super(EvgenValidation, self).__init__(name=name)
+ self.hists = HistSvc()
+ self.ndaughter = 2
+ self.outname = "validation.root"
+
+ def binning(self):
+ "binning for theta vs phi plot"
+ tmin, tmax, tnum = [-6, 0, 12]
+ pmin, pmax, pnum = [ 0, 5, 5]
+ t_edges = np.logspace(tmin, tmax, num=tnum+1)
+ p_edges = np.logspace(pmin, pmax, num=pnum+1)
+ return t_edges, p_edges
+
+ def initialize(self):
+
+ # All daughters
+ self.hists.add("Masses", 100, 0, 0.01)
+ self.hists.add("PIDs", 60, -30, 30)
+
+ # Daughter i
+ tbins, pbins = self.binning()
+ for i in range(self.ndaughter):
+ self.hists.add(f"P_d{i}", 100, 0, 10000)
+ self.hists.add(f"Pt_d{i}", 100, 0, 1)
+ self.hists.add(f"Theta_d{i}", 10, 0, 0.001)
+ self.hists.add(f"Phi_d{i}", 16, -3.2, 3.2)
+ self.hists.add(f"ThetaVsP_d{i}", arrayX = tbins, arrayY = pbins)
+
+ # Mother
+ self.hists.add("P_M", 100, 0, 10000)
+ self.hists.add("Pt_M", 100, 0, 1)
+ self.hists.add("Theta_M", 10, 0, 0.001)
+ self.hists.add("Phi_M", 16, -3.2, 3.2)
+ self.hists.add("Mass_M", 100, 0, 1)
+ self.hists.add("ThetaVsP_M", arrayX = tbins, arrayY = pbins)
+
+ return StatusCode.Success
+
+
+ def fillKin(self, label, p, mass = False, twoD = True):
+
+ self.hists[f"P_{label}"].Fill(p.rho()/1000, 1)
+ self.hists[f"Pt_{label}"].Fill(p.perp()/1000, 1)
+ self.hists[f"Theta_{label}"].Fill(p.theta(), 1)
+ self.hists[f"Phi_{label}"].Fill(p.phi(), 1)
+
+ if mass:
+ self.hists[f"Mass_{label}"].Fill(p.m()/1000, 1)
+
+ if twoD:
+ self.hists[f"ThetaVsP_{label}"].Fill(p.theta(), p.rho()/1000, 1)
+
+ return
+
+ def fillDaughter(self, p):
+ self.hists["Masses"].Fill(p.momentum().m()/1000, 1)
+ self.hists["PIDs"].Fill(p.pdg_id())
+ return
+
+ def execute(self):
+ evt = self.events()[0]
+
+ # Loop over all particles in events (assuming mother not stored)
+ momenta = []
+ mother = HepMC.FourVector(0,0,0,0)
+ for i, p in enumerate(evt.particles):
+ #p.print()
+ self.fillDaughter(p)
+ momenta.append(p.momentum())
+ mother += p.momentum()
+
+ # Fill daughter plots
+ for i in range(self.ndaughter):
+ self.fillKin(f"d{i}", momenta[i])
+
+ # Fill mother plots
+ self.fillKin("M", mother, mass = True)
+
+ return StatusCode.Success
+
+ def finalize(self):
+ self.hists.write(self.outname)
+ return StatusCode.Success
+
+if __name__ == "__main__":
+
+ import argparse, sys
+ parser = argparse.ArgumentParser(description="Run gen-level validation")
+ parser.add_argument("file", help = "full path to imput file")
+ parser.add_argument("--ndaugthers", "-d", default = 2, type = float, help = "Number of daugthers to plot")
+ parser.add_argument("--output", "-o", default = "validation.root", help = "Name of output file")
+ parser.add_argument("--mcEventKey", "-k", default = "BeamTruthEvent", help = "Name of MC collection")
+ parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process")
+ args = parser.parse_args()
+
+ from AthenaCommon.Logging import log
+ from AthenaCommon.Constants import DEBUG
+ log.setLevel(DEBUG)
+
+ from AthenaCommon.Configurable import Configurable
+ Configurable.configurableRun3Behavior = 1
+
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ ConfigFlags.Input.isMC = True
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion
+ ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Default FASER geometry
+ ConfigFlags.Detector.EnableFaserSCT = True
+ ConfigFlags.Input.Files = [args.file]
+ ConfigFlags.lock()
+
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ cfg = MainServicesCfg(ConfigFlags)
+
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ cfg.merge(PoolReadCfg(ConfigFlags))
+
+ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+ from AthenaConfiguration.ComponentFactory import CompFactory
+
+ import McParticleEvent.Pythonizations
+ fix()
+
+ acc = ComponentAccumulator()
+ valid = EvgenValidation()
+ valid.McEventKey = args.mcEventKey
+ valid.ndaughters = args.ndaugthers
+ valid.outname = args.output
+ acc.addEventAlgo(valid)
+ cfg.merge(acc)
+
+ sc = cfg.run(maxEvents = args.nevents)
+ sys.exit(not sc.isSuccess())
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
new file mode 100644
index 00000000..11a10017
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -0,0 +1,256 @@
+import os
+
+import numpy as np
+
+from foresee import Foresee, Model, Utility
+
+class ForeseeGenerator(object):
+ """
+ Generate LLP particles within FASER acceptance from FORESEE
+ """
+
+ def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, mother_pid = None):
+
+ self.modelname = modelname
+ self.energy = energy
+ self.mass = mass
+ self.couplings = [couplings] if isinstance(couplings, (str, int, float)) else couplings
+ self.mother_pid = mother_pid
+ self.daughter1_pid = daughter1_pid
+ self.daughter2_pid = daughter2_pid
+
+ #if not os.path.exists("files"):
+ # os.symlink(os.path.expandvars("$Calypso_DIR/../calypso/Generators/foresee/files"), "files")
+
+ # Set decay mode ...
+
+ self.pid_map = {
+ (11, 11) : "e_e",
+ (13, 13) : "mu_mu",
+ (22, 22) : "gamma_gamma",
+ }
+
+ self.mode = self.pid_map.get((abs(self.daughter1_pid), abs(self.daughter2_pid)), None)
+ if self.mode is None:
+ sys.exit(f"Undefined decay to {self.daughter1_pid} + {self.daughter2_pid} for {self.modelname}")
+
+ # Set detector ...
+ self.foresee = Foresee()
+ self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1", channels=[self.mode], distance=480, length=1.5 , luminosity=150)
+
+ # Set model ...
+ self.model = Model(self.modelname)
+
+ if self.modelname == "DarkPhoton":
+ self.data = self.darkphoton()
+ elif self.modelname == "ALP-W":
+ self.data = self.alps()
+ else:
+ sys.exit(f"Unknown model {self.modelname}")
+
+ return
+
+ def darkphoton(self):
+
+ # Production modes
+ self.model.add_production_2bodydecay(
+ pid0 = "111",
+ pid1 = "22",
+ br = "2.*0.99 * coupling**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10)
+
+ self.model.add_production_2bodydecay(
+ pid0 = "221",
+ pid1 = "22",
+ br = "2.*0.39 * coupling**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10)
+
+ # This is a bit handwavey according to Felix
+ self.model.add_production_mixing(
+ pid = "113",
+ mixing = "coupling * 0.3/5. * 0.77545**2/abs(mass**2-0.77545**2+0.77545*0.147*1j)",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ )
+
+ # Question on validity as FASER gets larger but OK for currrent size according to Felix
+ self.model.add_production_direct(
+ label = "Brem",
+ energy = self.energy,
+ condition = "p.pt<1",
+ coupling_ref=1,
+ )
+
+ self.model.add_production_direct(
+ label = "DY",
+ energy = self.energy,
+ coupling_ref=1,
+ massrange=[1.5, 10.]
+ )
+
+ return self.decays()
+
+
+ def alps(self):
+
+ self.model.add_production_2bodydecay(
+ pid0 = "5",
+ pid1 = "321",
+ br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
+ generator = "Pythia8",
+ energy = self.energy,
+ nsample = 20, # Vary over phi and theta
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "-5",
+ pid1 = "321",
+ br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
+ generator = "Pythia8",
+ energy = self.energy,
+ nsample = 20,
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "130",
+ pid1 = "111",
+ br = "4.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10,
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "321",
+ pid1 = "211",
+ br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 10,
+ )
+
+ return self.decays()
+
+
+ def decays(self):
+ # Set up liftime and BRs
+
+ self.model.set_ctau_1d(
+ filename=f"files/models/{self.modelname}/ctau.txt",
+ coupling_ref=1
+ )
+
+ self.model.set_br_1d(
+ modes = [self.mode],
+ filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"]
+ )
+
+ # Get LLP spectrum
+ self.foresee.set_model(model=self.model)
+ plt = self.foresee.get_llp_spectrum(self.mass, coupling=1, do_plot=True) # This is just a reference coupling
+ #plt.savefig(f"{self.modelname}.png")
+
+ def flatten(l):
+ return [i for sublist in l for i in sublist]
+
+ # Get list of events within detector
+ coups, ctaus, nsigs, energies, weights, thetas = self.foresee.get_events(mass=self.mass, energy=self.energy, couplings=self.couplings)
+
+ # Return energy (converting to MeV), theta and weights
+ return [[e*1000 for e in flatten(energies)], flatten(thetas), flatten(weights)]
+
+ def write(self):
+ # Write LLP results to a file
+
+ energies, thetas, weights = self.data
+
+ dirname = f"files/models/{self.modelname}/events"
+ if not os.path.exists(dirname):
+ os.mkdir(dirname)
+
+ if len(self.couplings) == 1:
+ filename = f"{dirname}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+ else:
+ filename = f"{dirname}/events_{self.energy}TeV_m{self.mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+
+ print(f"save data to file: {filename}")
+ np.save(filename,[energies,thetas, weights])
+
+ return
+
+def setup_foresee(path):
+
+ if path is None:
+ return
+
+ # Add foresee to python path
+ path = os.path.expandvars(os.path.expanduser(path))
+ os.sys.path.append(f"{path}/FORESEE/src")
+
+ # Symlink foresee files dir to current dir
+ if not os.path.exists("files"):
+ os.symlink(os.path.expandvars(f"{path}/FORESEE/files"), "files")
+
+ # Install scikit-hep if needed.
+
+ try:
+ from skhep.math.vectors import LorentzVector, Vector3D
+ except ModuleNotFoundError:
+ os.system("pip install scikit-hep --user")
+ try:
+ from skhep.math.vectors import LorentzVector, Vector3D
+ except ModuleNotFoundError:
+ raise ModuleNotFoundError("Unable to find skhep. Please install the scikit-hep package")
+
+ return
+
+
+def parse_couplings(data):
+
+ try:
+ couplings = [float(d) for d in data]
+ except ValueError:
+ try:
+ couplings = np.logspace(*eval(data[0]))
+ except:
+ sys.exit("Unable to parse couplings")
+
+ return couplings
+
+if __name__ == "__main__":
+
+ import argparse, sys
+
+ parser = argparse.ArgumentParser(description="Run FORSEE generation")
+ parser.add_argument("model", help = "Name of foresee model")
+ parser.add_argument("--mass", "-m", required = True, type = float, help = "Mass of mother [GeV]")
+ parser.add_argument("--couplings", "-c", required = True, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
+ parser.add_argument("--pid1", required = True, type = int, help = "PID of daughter 1")
+ parser.add_argument("--pid2", default = None, help = "PID of daughter 2 (if not set then will be -PID1)")
+ parser.add_argument("--Ecom", default = "14", help = "Center of mass energy [TeV]")
+ parser.add_argument("--path", default = ".", help = "Path to foresee installation")
+ args = parser.parse_args()
+
+ setup_foresee(args.path)
+
+
+ # Create PIDs
+ if args.pid2 is None:
+ args.pid2 = -args.pid1
+
+ couplings = parse_couplings(args.couplings)
+
+ print(f"Generating {args.model} events at Ecom = {args.Ecom}")
+ print(f" mother mass = {args.mass} GeV")
+ print(f" decay = {args.pid1} {args.pid2}")
+ print(f" couplings = {couplings}")
+
+ f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2)
+ f.write()
+
+
+
diff --git a/Generators/HEPMCReader/CMakeLists.txt b/Generators/HEPMCReader/CMakeLists.txt
new file mode 100644
index 00000000..3076e7ab
--- /dev/null
+++ b/Generators/HEPMCReader/CMakeLists.txt
@@ -0,0 +1,9 @@
+################################################################################
+# Package: HEPMCGenie
+################################################################################
+
+# Declare the package name:
+atlas_subdir( HEPMCReader )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
diff --git a/Generators/HEPMCReader/python/HepMCReaderConfig.py b/Generators/HEPMCReader/python/HepMCReaderConfig.py
new file mode 100644
index 00000000..72081a86
--- /dev/null
+++ b/Generators/HEPMCReader/python/HepMCReaderConfig.py
@@ -0,0 +1,24 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# import sys
+from AthenaConfiguration.MainServicesConfig import AthSequencer
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from TruthIO.TruthIOConf import HepMCReadFromFile
+
+
+def HepMCReaderCfg(ConfigFlags, **kwargs) :
+ cfg = ComponentAccumulator(AthSequencer("AthBeginSeq", Sequential = True))
+
+
+ from TruthIO.TruthIOConf import HepMCReadFromFile
+ hepmc = CompFactory.HepMCReadFromFile(name = kwargs.setdefault("name", "FASERHepMCReader"))
+ hepmc.InputFile = ConfigFlags.Input.Files[0]
+ hepmc.McEventKey = kwargs.setdefault("McEventKey", "BeamTruthEvent")
+
+ cfg.addEventAlgo(hepmc, sequenceName = "AthBeginSeq", primary = True) # to run *before* G4
+
+ return cfg
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index 709797ae..d6e254c8 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -82,11 +82,23 @@ if __name__ == '__main__':
#
if ConfigFlags.Input.Files and ConfigFlags.Input.Files != ["_CALYPSO_GENERIC_INPUTFILE_NAME_"] :
print("Input.Files = ",ConfigFlags.Input.Files)
+
+#
+# If so, and only one file that ends in .events read as HepMC
+#
+ if len(ConfigFlags.Input.Files) == 1 and ConfigFlags.Input.Files[0].endswith(".events"):
+ from HEPMCReader.HepMCReaderConfig import HepMCReaderCfg
+ cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
+ cfg.merge(McEventSelectorCfg(ConfigFlags))
+
#
-# If so, set up to read it
+# Else, set up to read it as a pool.root file
#
- from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
- cfg.merge(PoolReadCfg(ConfigFlags))
+ else:
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ cfg.merge(PoolReadCfg(ConfigFlags))
#
# If not, configure the particle gun as requested, or using defaults
#
--
GitLab