From 030fe1d8d1a4e619069d9a115b302fa6ed13886e Mon Sep 17 00:00:00 2001
From: Unknown <unknown>
Date: Wed, 2 Oct 2013 10:00:25 +0000
Subject: [PATCH]
---
.../DetectorConfiguration/Riken_65mm.detector | 31 +-
Inputs/EventGenerator/10He.reaction | 15 +-
NPAnalysis/10He_Riken/Analysis | Bin 66916 -> 66708 bytes
NPAnalysis/10He_Riken/src/Analysis.cc | 15 +-
NPLib/MUST2/Must2Array.cxx | 20 +-
NPSimulation/Simulation.cc | 1 -
.../include/EventGeneratorTransfert.hh | 12 +-
.../EventGeneratorTransfertToResonance.hh | 229 +++---
NPSimulation/src/EventGeneratorTransfert.cc | 7 -
.../src/EventGeneratorTransfertToResonance.cc | 744 +++++++++---------
NPSimulation/vis.mac | 2 +-
11 files changed, 527 insertions(+), 549 deletions(-)
diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector
index efc969276..0c25c3e98 100644
--- a/Inputs/DetectorConfiguration/Riken_65mm.detector
+++ b/Inputs/DetectorConfiguration/Riken_65mm.detector
@@ -84,16 +84,25 @@ SILI= 0
CSI= 1
VIS= all
%%%%%%% Telescope 6 %%%%%%%
-M2Telescope
-X1_Y1= 155.77 50.23 8.18
-X1_Y128= 155.77 -50.23 8.18
-X128_Y1= 133.17 50.23 -89.7
-X128_Y128= 133.17 -50.23 -89.7
-SI= 1
-SILI= 1
-CSI= 0
-VIS= all
-
+%M2Telescope
+%X1_Y1= 155.77 50.23 8.18
+%X1_Y128= 155.77 -50.23 8.18
+%X128_Y1= 133.17 50.23 -89.7
+%X128_Y128= 133.17 -50.23 -89.7
+%SI= 1
+%SILI= 1
+%CSI= 0
+%VIS= all
+%% Alternate Telescope 6 %%
+M2Telescope
+ THETA= -130
+ PHI= 0
+ R= 150
+ BETA= 0 0 0
+ SI= 1
+ SILI= 1
+ CSI= 0
+ VIS= all
%%%%%%% Telescope 7 %%%%%%%
M2Telescope
X1_Y1= 27.07 50.23 -154.49
@@ -104,7 +113,7 @@ SI= 1
SILI= 1
CSI= 0
VIS= all
-
+
%%%%%%%%%%%%%%%%%%%%
AddThinSi
%%%%%%%%% Det 1 %%%%%%%%
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 8e2f05efa..b84050210 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -2,7 +2,7 @@
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Beam energy given in MeV ; Excitation in MeV
-Transfert
+TransfertToResonance
Beam= 11Li
Target= 2H
Light= 3He
@@ -10,13 +10,16 @@ Transfert
ExcitationEnergy= 1.0
BeamEnergy= 550
BeamEnergySpread= 0
- BeamFWHMX= 6.232
- BeamFWHMY= 9.069
- BeamEmmitanceTheta= 0.01208
- BeamEmmitancePhi= 0.01681
+ SigmaThetaX= 0.6921330164
+ SigmaPhiY= 0.963142053
+ SigmaX= 6.232
+ SigmaY= 9.069
+ ResonanceDecayZ= 2
+ ResonanceDecayA= 8
CrossSectionPath= 11Li(d,3He)10He.txt
- ShootLight= 1
+ ShootLight= 0
ShootHeavy= 1
+ ShootDecayProduct= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%0.69u
diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis
index 6c0904b1217ac65c52f5f7efcf964871d826d228..7a15c0ec1ee2c7c99032279403bdbfc030232674 100755
GIT binary patch
delta 24595
zcmai+3tUvy+W*(UY#b03MFj<9R8%ZfG&L$SEHNrGEHx|jkjIQVR@7@o>3IiJ3_=QB
zYM4?~W~NpY=Ak0Pyy2yyGNZg<$IFmbSZda>=KuSx*~?i>@A>>UpT)Dj&$>P9y4UR4
zvyMLPTvFsLjSUw&oUXw_2-l+Se;;1>vJi9@6e3Cp2ZYiMvRs95gF}Uggg(p;4}n~T
z_^S}MYS5yOhlUR>5n^x~#FP%lV_mN9b^ORbT0PpjzHCI+vHhQY<ou?5(LAKNQ7&>6
zdv1hpG~CT_S#V?Fc-;zj8{9a!32?W=@tWut>i-@1f0xaZA@8>R(`|V#WDeX6n<*YJ
z{P222uB~_&@)6trKel`nG7s)?n<<`vn*-+|;Cj+8L;>UixTkE^q7a@UxTkGJG2}wH
zXW>ej%OZFr+)}tQI9|*BLaeYQ#Y(sr;9i77bBkBtR>8dr_Zr+9IDM_v9$T)1d;{(+
zxJtOU;ogC}Sc@tGu8nZ-!M#tw^#NQB+=p<R;68%e4EHhIR=90&+u?ZC!R?0I1NRBs
zr*NOa@v4XW0`5z=y>MT_@lu2A>#79jZ6SA`5BKZ%zXP1NbxWS-YRR9g5$M-3eq&I0
zK~VZF^7F7n$Gh@0{&$z{TP1q_Y}GQy_CXd@I9rCF4&niUVBv=UhM@3wRQ3FkXkniS
zzn}lXJ;b4)Dx?STPI&}o?hDFhX;4Mq4@$K#DE`Wz@bsW+whd~Al`R`idR*<>7q!%h
zTKSSIam0}zZG2DzToM$1Rm%pz8sG>|99SP@&kI5J><!8#R*nrH;0Td(!xKk*5)}6v
ztWNr3ljjGOxIc)WmPcTLH%R+iPz#?3(RuoZ=DTuWMB>25f@<3{$d<gIesIaT5s8jx
z<x3HXp8pA|^Cv+i-X6sN42r)osH_L0h4}gxr;AP1L$nchb`v7|tkVUx;%mAK(H}g}
z&mOTc3chBii)$GZtVf2mLo659ruL$%xGi3YYJ}4chJP0+#J=yHu8ubUEg5U!xj+F&
zx(IPT)N-*++KYA~DFJKeP|aP%2$c6O4D6xOd)kO&4lE;Hr;EMB0^Ub~2Yi7H3p!y5
z2(h%q;yHxZZ+5D7=^$E*FeL1c0@=t+e|ay&M}A4jm%JJ^-hRgEinGIy^%TOLX1Rue
z+j~$$E{OM>2-I+YXCeOAVY%Aa396AYA2mRJszpj*Ko$zT6Wjp>^hCrTj3$2`TxYik
z!+%BF*P;b=mwXe8HR^J|0UkDSG1}rB3`n#KJltN0Q-3;L`E~)rAwNCobfG0xjXIzU
z8ZHFF&!fO>7|14Je*MtIwMLU)+Z&5C(sOXLMJ_`Hvi&XKQ3<|jE5uOrX^I^oxgG9P
zhQ$T2tO6Oj?E<{$<5_58ZSbq$Y_ym*Fcd8~(J-(I2KGicdPW)W3JjWq{M+!pLyOW;
zV{Sh1Uo~kt=Gg+Yh;Fj}XrZoXQnm!+??oTZF=}`OqdWqm8fq2(A;!c~6yP6ZZNy+0
z^c}{4Ht_ivAr3Krl&=z$A!Mjgqut<%D3BJhz*@A(lG9GtXZE1FfIggH#GeNP&ioWG
zxDnGdb*rurMm#L}9xbd};4E680tR^?RRPyxkZnK{jkozlj<#)q45q@MhCiGx2{i}R
zILyDjFhI|1H}u(dfvok7Fz`;y8a*qjF=jd&23Mg4G8&y8*Nt|DKS0spl+%TZDuZ^q
z(PYM`&P5ZR2d}Zi^AX+^vp^f*Kmj#aHxP{7;%VOywjd172h-qxBjQ%0%Z{PKIYT{x
zQM-pNU`(5H;6gMR#*Q+mRTBCbHOjL2ISkq|7=&_E_;46-9%DfF`H$TZj{I~BT!#3S
zMi-nyx{o}<a!s@|VADQ;L7@};9wNkA@MU)R9T-&8j20=wXxj^ebPF{0!NOwCYfi)O
zU`eiF;H?Pnffn8t>NogC%o@)#uu#v-d0o+m$Z)1z;9D?Yu`vryqR$s%jOYTtM3*f?
zjgU;Wz(0_2>6XBlslynWfO(IGRN?W6_z3t`n;+<j`PcBDKp$O)3?h+XUpvA!%)bl9
zTHXdtyyV(IA7728YBWaqm6!#M=wr3Cq6HSA$qyhsvcni?D@I}|AB_s=dxs+q^RL-h
z>mNqMM^Iy3!{1?mIjG99a#a{jT!s~>#ORWzP+&R=WEXPK4MIk-C=hDZ0<(LeJhb=;
z+?^Oc4C}9`wnwd7bP>NJLbg2!Iqzp7WiD#ODzHyq!}O@I!<oVLNMC`a^eE;#!)K!f
zmtmTv*ad!psva=nw{r{O#F*0+^u(eM_h6-}vlArs6k@6|O3$JNHtcY^7GXZKzzo!Q
zCNg9VX<!Zv@L?g*E9E%czMGA8<4eq%DD-)T9lsx!>Oc#6#CD9nMq|*ti-`9cvtlh;
z<atc1$#w?iFr>m5gqtx)UoZx#52M^`tb`?KvG0xWJQ#Ee2Kg7hv^L^#G*K)Hyu~iy
zN~24n&?WytpU;E=dZi3Um(2LH>UeZksJ!Nn<f=zvHukEz<mx9Ka?J2B$2$4;@aQ(R
z$;b?s1AouR&!O36XlBdPE#icT4`lmmq7%MA9PWmq?z<0*PN)dv5+X|Eh-;!fhwK<L
z`!03~5h9xGTxSfnTw$=A8z*ZF?oHrUn0@(_XrTF&7#f9?7<&sT(S@ayVVG2u=(#dV
z<h+s+L#l!j9khzFHD>D?dN2#uQ^LqfN;Jp@O0;w}CF)Q^NrNb{Zf&K+psS-qgYTil
zY^<lm%D$HpL$QGpJ$!)jb&L#344I>pT`^m|^k6|bMTwqirrd*@A0?LJ^ORTyE>L2u
z3baToHwHcBlX$d1*#pZACB}LbWswjaDKRu+DKX^YDbeL_N(`waN)MVlnV!Bvq)=Xh
zyDViT9yCy1DnuIP0wIP`Vj8AX_7h?xWq%<uDF<L3qD;k3j&dL#Fi;M{J&W=(Atq4{
z#$=%!f|VcA(`qPgT2x$)F-Vz)<$>}Fj8V#+*z8hXi7`NV6~+K1R>b*~S7QuN4#yav
zyar={G96=pas<Wz<+T_Cl-FSlP>#eH7-D%`tulmIL&YeJ0m|z!1}HNz1}JY3VtmoD
z+vc5edeQjeuDM&5JR(G~yx{Zs3ci?WdAng)dY5DDc{iZ2iFfVX;u|s(o#pG5hyP6D
z>oo_y9D17?k<8bN!#D?F9??{<6dcbvoGMMViaD-xFl<CqrDBfn95=Z_aVa^OT&j2p
z8I3QR@)ehn(UzhqSMf@6IyqZ0=Rr;;cZ*F~O2G+{lf?+>ia9TGCX!PXbAIGxlii9d
z$*fXSl;REKT(VGHP0k}X{|QB}8gf3lQ8DLHP9e%i3n&FAR8A?mRxu}5P8qpUF(+6~
z1-U|TJ$VheRPkPNB{^Sl1G$=<tM~x9mYl7a6EVkAM@5!Wa7O0TlhYM*V&*iEQx$W9
z<}{Msikrz^a+KmT<Yux^e4c!s-24Yy-~t)TCt5(U!07Vi@Iz2jy;3-#$cZA?Dt3`$
z$(4#D$!>Cm;wW-5xm0mSaw<7raV$BFoU1sVoKDVG?B@8-$)qAnDUui=i=3`FnLLr4
zsyKz5O?E3zCFhW%6b~lnl7-?lavr&vv!&NCaz43HarzMKR&xres8@=Sj8IChRh&sK
zBUdUOO|BqUD9$3UA(tv1N3JC2E1pQMCg&=iM6M-gE6(;%QAb6VQcNS)lhYOFkQ>OU
zif51;$!^8DWG^{N@l0|vSt!mUpC>o}hb=INjN2SqK(U8wS&oZ}dZn08jw078E+og2
zD-|y!yU7)bOUcRPQpHQispNdcW#lw+uHu#CbaFP?(~FB?P9_yuD#9vq7CBw<8uCPP
zs^azJY_eN%B{_#2rFa85mn^^of0$BqWZu!{tZ}#Ia2uJ!kNIOK<ZuJ#_5JRX+1GdJ
zG#cAM@#u)+PFe6&uf)^ter1E?GuOv-p~Ba3X&F@h&?~-u<oI@L^ia9!`YxR-pf0GL
z*@2T>z!Sd1nN%@I{&0QYAugxo0@cD8?=IM}V`P$(o46uplJ^qK&JUQRX>cqyy9IY(
ztWy87yeYG1!ZQdB=qeT6X~&o<i!*x;9|PUVZbe@fZ4bn|AG)IEqHmecGbpd5K&Ro_
zQc`f2okPeU@>FII&yT>9asDiH#0Wc!rass;1XBEhs3irXnI3g0DVSvEbm?z_EOhQW
z7zQhvYOnzbL}kYYXzCw6&&9%ZZew)HyO5HN^(-2&q=4%LtMr-Oge#Rc$JV~5@>Ow^
zX<C0o=`~gT4nt%8b+vUK?~72<(vwkvLW&vDZWLYcV|Kizt=L+jL(|Z0-gi*-qAf)Y
zv~1sb_$QZ^6pW)$$_y@c3~G&_lbwq8740aBPxAI-qVq*RH$Bk_Uwf1ks91<!oKeci
z0}$C(f<|ESYEyJ}{h(8gd!^qLj2+=kHcW|<M{ewT$J}4dQeUBYEVUV1#FB!Uwvojd
z_<AVGYhgEKD?$H4v*|)rMCe|DP8C|&QI5MQvsaDkfIPc~-=TDN)-I*Zx3zoZ=9^-|
z&(JdDFiIY|sf*_;##oxh7~_&$O!F-?2={$nQlJx-u)*w}(z?f40lTAgvG*{Z-BDUM
z){fWo0`|nr-*0K5-=<h^hZd%m+NNHRi$=$UC)uXXl5dah;yI~IO*Kqy4w$MFe#lm~
zP1U;Bq0^??<$9Mgo^7huv4>Qfeu`~;z|=Br>b0m!pj9huQ&Z&KH^+pJV>QvLtL4}+
zT|7T4Q<Du-e+`(b6aH{6V5-)A1D!V2F4wz_@oZDI?rl3>)1Ua{5ioU)HuVuyC17f$
zZR#XBa7;}2d{z^t?vveb>EiiKnd&x7{UcziPMBbrs&(Opsdl;Evu6XQYTXg&RIB!H
zVQRHD^=0G}Ftye;b)j_J5);0GWx&)hWvWk^8f%yuYF!+wI^oSU!5#owcZF>N7Ztl)
zZ!g1Ct!uAMof2iXYMnOq3xo$ut+!3Bk)OlVqf84^hsx7vRSSldEKN5|Z5uFEC;aD`
zfT>z{A9Q*E*yVc1GoC#Fv~GkQujv_l6b%f325o9c3sW0yQ%}j4vSPwJ+NMsIo3pxj
zVw9;a!_@8pQ+2{OehHYWb<3gCrrPCt3mDHfRqGzM<27x<r{{pFUTx|iR3*@=&9<rW
z^5L;D;Uig1wCXasXlxfxUu9~PVQN~yRGskP&jC}lZYOlwRJ&a7yNqX>s&y~f@tS_b
z2mX=*pKZ%7)BrStIMTmLX3&&^^LB*yVP=5lFc%p8G42ON87l(;^AJssO`Y^1){yJy
zut0e$pp%Dgb@$FfDa9Epk(uucZL2owE<0gUvfO`bPse2GyY<rU4;cw?SKPo`Q7-N<
z`mT7Zyy>=n?iSfzyb=C*8kX>u^bc8fTQ6j|>9+2JDq35vgOh7Z3R-a&f=XRLESv16
zcvLk7KiV+-lx#n)<KKNKvlx}HD=E+c8`$3V#MHV9yS-0#4HU2kx*ZvK&VhZ{Ts)BR
zZxD3xZ`vuEX3N5HSNB%?xNCp(51`pFgeG+O$<KOU_lP_>Zs_Qz@Fc^&y?XJqqd23E
zafcu+>Q`LgL|1YYvRnpyhkL-m$rXPMXhMKoKE8YBP0ehH6x6$Td>!tdPA0oT){jr}
zbo{~R!*c<DDAoUlLRX$CsQU~$)$_?E1uk3n0pnffuW}>z6^Qqmj%PQ8j{B^wZ`u=%
z>X;Sqs2cQU4Xicjz(YF^POfYbLq+J-B0}}SjLL)94?T|7J6K+d!VhLtQ{g&AMVV5l
z7{8u`qJRoyTwW1~Ib0=Ob0LChRn$*p>g|^f91DF(K@B4MH)7oSvS!QC_$39Zo@jY&
zcx?S#x$ySn;W;Q*x0wzZ2c4eCI+t|DvnR6FrJRyK-`=N-3@V_;E?~bLG%@+=k5La@
zfDU;RI$c1S9s7C4vkTC=1@gIx2?^B*QT-gKgA*T~ShuTxke^TN(>KW|LkC3}Wi;Al
zT#()Gh)(o=?=QV+C!QLW6eJ@rV{CmdZ@c4?=ob-Z{{D9Jf5S(qgBi64t%}yN^p49t
zm7J-Qk|L>_l;p(h=fh|3rY3a8<l2K7jnH;HXjT6u8lUoep9_js`<G}cw_F{q!Jk{*
z#b`rysr**LdqYsX`it>&G&f7_(In^G*~N7|M;Jz`EWh(+#~1SComYjRQlC%i64KCG
zo|@Ff^Ixp#>up)L{}L@kNBfJ>+8AZ?{u6d?cC>)a?Q}e4^oA44=x@*u%IMaU5*>@>
zrIUsXZe%UX)3A~pwETXz;g5wMAJ|&<>FaXIq)VcwplJ5#-FP6N`!rDfklNO=`HoH<
z@{b2sm+GuLwEDkSk8M3DC|aO;7i))aSG=u{tJ<Zq1KB04U6b59dB{Eg|Ju>?@>Eti
zquyS5emhv63`O9!pKH&LrKob?uAgUXJ@UaRF^<>df+<PEr@~-8J$3jP==Ai|IbCVQ
z)4JZrrFTmF)v?C{S*Ya)6VnNu&O$G-%|`<>T<Z?V5r50{OtJS-CAhcf=LT;fcJcT<
zddaW-x4$OK?b}Z&f>|?t`)P)K7weyq&&8pQtx})QG772s9cGmvlnbF>lZ;4rcz;oP
z73Dmw?4)iE{un2qP*rKel*nYC5ko1q=@^~VgAMz60b-t#cV|am-GCx=TV$b|n8#M=
zbhB$+r4dhko<_Wt^6l*CHidQwQ=XEaXAe!d2Spy7+*q6uN!{(B12;Dw_-6xeCGW|k
zyRY!f$M>4v_m7}eoJj@r_eUM8`5iu*$tAGMXLB)xISgceR-Maoebaa_p%7Yzh6)Ei
z_VZN<??;GH&2MAP<P2f*pV>Rkq|sO{y=h_G33-Yd4q9ISf)snuYWDl0kT}wdFOHbW
zZuT}xIc9h1CX`Ur@JcQ6L|!wP*AG~sxvM&8o%d&u@{Jd3ViJ<M^Rb7xsE<{8l|A)4
zMtYK5r|kSg!~916)Tr9rej;%<Oy)+6<`4N!mCUd6OEUgNO%+1vwRc|;$f6{pk|D3d
zO7xtaOPrC*Pi*XwNl~bt*df22+Qp+*ym2VjPJ08>+G9&)mc}Y`UwVsO?!8b{^DAgJ
zknSXA!=j&nnI#z?(d)Clwtcb6zG{Y?g+=y#G}>R6j&>moO_QHaP4HB60-%l6N?U%I
zxvI2_nbx*XDdw_zwtcGHN+|wfAG2ZRThIU{89V6BRQB!A?blD)SIdy$|Mtae1O0WG
zY!|}NxsCF(sh#D3d&4}nwvpLfSC1cJ)+&Sj$XXkz6uY3%MykSpgW@knG8>i<sf_%L
z-bJ?8Zp1WYWF14E`pz&?<+<7}grU0*Ya?aUJuXiL4IH2hREd)e1C=7$Fi;hI?f+#U
zvtjPz)p+@u-fwNMZQymvz<P$%!*sisJ7%anIr`WN42?c4*WHukX{2F!2ibBe(F06m
zcdt@RveR(*sRE0k_)E(%8|M0g8YqY9?Wqhqth;xdGE6-Y?7}kq+qh_do@sU=OgBlc
zy(dBb?{t@kZ>)JYAE0e2-A3f7ZBvR@q0qgm%G?XZUu<JG%=m<A{uA`hP-dRcw#`ts
zH8Nx>%h0x|Jm=blF!WV<_w*QyytW>`OXiJGwyAVJ<fv^^iXWj+w%t{dp~_77zieYR
z%s3Ucah<dHXa=h(bDx78s>QaJ^PC|_km&WH&ul~e`ERldWOm=l{r8RxNjO~99bX1H
zc2_-lUt32=%C~ay{b^wvu;h6^l=b(28Z!KAS@DnWLhk&!YW|FQR*3agRrn*FLmZz-
z_sl+yuVm)TI!BFc|7hpV1@)*59;bUJ(PNLGe0lYw<3jHGRLVz}IyT6@v!;d&{8X0A
z8r!Mu9^RKDi|QS<j{41wp)~1!`OB=qj$<<UvDt0z+kJ7SPT4IhAM0|(yie85K}5o~
z>pwy0_}zRA%IRH_QN#2Pf5tSs>0PG0^<4g3G9>S+kXt^J*W|_D?O?);J&bZ3`RVRZ
zim#wh-GMu_8rNGH!~ZGRuc)$PG~I{Axujs#pGbh5_UuAk{I}4DWldfe&ts@uAovcH
zIlp%140c{0SdY$AABD%Fiqz@7$J>ZU&M@NiLmW1f{fTu93}BQy5rrkU4a(v{WYILj
zgL@uhO*1k&MjHbeh1nTBX=k(@!KN*}nbCOUQ<AY&weZEtF7y`}S;WlIyx0Y+!XQ*+
zST*n)<W^AIZb|L@+I9sS5Cah-i|qt>Oc0#lb?JPf|MX>sxov~&c@QmX+Bb#lOd5>5
zyg%<mH`w1KEPx*0DP$ZCG=DU*LO<Jg+DB#C6WvGjtm81)%bFEu97Rj|z9ho-$B}{V
zbEU{(ju_5OCsCH&%Y*mG5A$OCT#1}ZGSmt%0Ifo=y~J)nRqPe1*B*Jr>>h(pBZ3}n
zYGG=GLXS3;MLU*kFRM!Nr+jvHY}ajwH?P^jmTbbVFHo-oa`WuI(>~sAXA>6)c@sK&
zN?d0YtLDfcqgZv@>u$tQiYOhUGy5J*;*(?MbZN^0q^I=_bybhf=@8=Cwu1%z+OpQo
zO!EO6xH!XAGQL?4-Os=00E{Y%D#=hWiV#CrS1IzKPz}<7St5_upqRJSeb<!9Tm1WI
zrCnr4LyQpbyHzevyyNmxZcP6Nwf@V-lp)M();9n65uGhAo{N~^V|iEp8s~dk@ciy?
zbG!7f*^KfVXrevYK0wIHu4>?@9&62}wzql@$uV;;889t4$#q*>CaFTPI!VkH`P$qZ
z#{*SSPi8wpMsJpn%qx!Qgja0lor?F~yGgd1KQ^)yUxv;*@0%I@LFphlZGN9tg{c_e
z!?si{p1-<PjDC?715t?dXKSx0zxCx&hZj}dzHqL?!$3z|=vDbxnjS)U{`)>(CicDU
zzx4U?A@73p;&r%A-}-!!kwOfFtU%l=AZs9pLpDI(3V8-{3S<<Xu-y-t0yzsZ6LKx&
zG{}vR`H&w&u5{uZl<(lFh9YFY&$kz{2V^tkwUCi8U@GKF$Y&s{AuAyFLT-R;hHSu{
zvLiPAA-HR$LiUFo4LJ@n2XZRpe8@SF6_77O)<ABAYzV_UC<ov<14R$qD5LP!*l@@c
z$WqABkk3QrKvqM}hdcyX0r?|j4de!Frw%~YLY{}*1KH7qSEI18Oof~PIU15%e!liI
zA96nALdXgaJZ13IKxSj7d;l`_xX*V1au_6jMJD*+e=y`U$SlaYkTW1(fh>f46LJ;g
zdyut|{B(B!at1zr3oM(9A>$$6giM3n4LJ^SAEYN2p5NhF2$|yb`PM_GL)Jk~hdc_o
z6_T$+9*0bVOu~bMbjTFQNs#G~d61(Zmq6YLSqV7}au4JT$dMR3`H<O=OCaY!u4(dl
zd}Z*|K(PUG)sN^C$Rs>WJVS<z!ejcRW^@VUD#%R8G(4c50a*=MNX8T6Gx(BZA!HQ3
z*I5IZ%zPnNA-o>47V;?M0m#enupckbisD~9s1ZCzp~!^%8FCtAEFQe)L#9BkguD^5
z8ge4!UdS1c&5-$!ktlc{BwzFIjYkP1AyXi;AuAx~LvDbqfIJRa1Nl2-17u_xzWo{r
zi{l})A^F8JFUJ#Ew+hD`SD1R6G%nKVEOao;4Og@e-61zGinj0`oOo2ev?$IoMSh6?
z6J-0vF|v7)(=m!ZYdyk;!WZZ0E&DBwlV=vi%Iw8X8MQdpipTx6o16>nUeGvDZ%~Y^
zU+i>TLZ7t|y8iO?;@)ocMyuZlyc&p1li*U{_4)W7t+P5oUM~AtMY#9hsLYDT24It%
zE8`rW%5oWJZT-sUi<WPfMaz2WjHm(e*Cftc<RRJBszGh~GMS88;!M~DN<rpgL|fOj
z1tD#$>*s~C-g}Vcw{qZ;ILA?W^O87gAL4aYsq&XN9eWWHV-3Udb4b2~IQ!&>_}?hM
zSrV5J@wLx~H!K4=!izU8#d^e!knNWyI)=%COXE=IoAEzI-nTR^VJOm|cm9%Q7}ljr
z<V#3XFE=B0E&j(?OQ5Zkzb%ckqOc~lp;h5~pyP_xS~Z9lbXwlLEY3QQ)lT)D%wOh2
zT{bU`4yQH4u#)YSFD*-S%#oY%KSdr{7UzhUf8c)t{4p~1Ij0qgr7A(?cna~`$g$7G
zIieYAC1C*zQSmAe?=;rL)@<;0Tv;1IYfy$W$1QIu?<5<ab0$nfPQA?<+U?9AwLK^u
zWpR!NWdi;`r1H!vb0$ndik@Z)_-T10mb~fm;j%<Wo?L_!x$=#&xL!{nHs5h}E_7Qi
z-Iq7avlcj?QvF%Tg%kG2gs+0{eZ(HFdegn!X>mdJk8dukb#laVw)ES}qh;=Lmt!qT
ziIJs<$CW%nwFZ~hQpAYs`;Q$yAAh*#EHlczEl}=4#B^X28?Vx>0p-b;%44jp@VQma
z>CZbIvGh4E$Pv%SSzLi#;AmO;JbrM3GWoi@GrN1BUbh7DpUCxBe)xQxBU!dzfdRkB
zDVvewLQU3s#ApV^wH}9cU$TvU>dt<#PN}qGS1|1ZE8;?;j>%1nI?IX`cv^<1XKDtu
zt%jijbPE$Y<*pT7BQ6C_Bu2>}S9I<5H;wKd0TYd;=gPiB_FEa}SRij$8E2Iuen-_D
z@mQ-1<y?e3Ez4KNh0O6AP`}a{GE;f2g<Q?$X~bA9+my#y-y*H0W_EhHGh!=fCJnhk
zjw$cz=*X6p`Q;e5xFZczMGV8L+ynYJM_2hFG8l%N<v51O=5nWFraW7I8Tz{K3vu|^
zm+`{Y4=hd<PlZ}5L&efi3;sIw|4^t^<rL>atzxH`@1%ddlc9Lg1^<sYt$CKXp!{nr
z>qARyw5UI<{-3n0k}$C#jGA}CtQ}#ZR{7s;Wi6IJy%5t1`hsxy+q@X#f&Z%r#ylBe
z{nlE{Yd!q))?!y{`hRXs{i~7u|3Rd+vyJ#Na`;Pa#Pe<F|F{kHXVw4Ww$}1?;+3|;
zPq!7{EC1ql4Be*w?{8-{MTuY94c{Ck-jAaHyC{Y(ZEr2^mMYe@x8Ca@wzMC<sDpUA
z1O4xGVDK^gmz(`2e$vtUtdrQ^(fUtEaYp%fcbexwq^~+#i(|wyF(L4Oss68sweVZ3
z%2;bxtk|ObUv=5ZsBgwu>*B?ian{K=aWu}7@nT`Tb*P{CI(~x#;)QP3#(rXBH|v>1
z@o6{fcsFra`JeJ6GGMFve<0D?;1*{Rt%Yt;;HLi#H$xAr|G&GfFS?7x-L19V#Vg(E
z-`$;|e{_F30-`L%n%_@Uro<me5nn0)A1MretDiEXA1l4Q|1A73A7FV0h_?nrY#$&#
z>Zahx00tMQTJNWds#NQ{RIxu53*>)NtIkxUhqSr@KbWZXIptNa6nK_m*~iUCyll%i
zY{|VW^?Pi2B_=t!(U#+pANd|zhG3B)ceZ7+EhpLXaa&@rsd%$7a}eJnp7DFcI$M4K
ziAkY;go&+*`i0NgKYb4^oa_HtN~3gQ7l7U20{>5&&L#`|Ow<3<3cP%vF4!z`@rE&O
zk|GcP@Xb7aKO+QQCKMAta=PfxgzO-C`V-<kEy0E4PoM9h!c9(>4-yr_uZMgO?UP4V
zUGCwz15l`|^T+tT2;5M(QE(IBro+vGTL8BVZZ+ISxb1NJ;EuzcgA2nm`_6E^;fBJE
zf|~$09c~uf0=Q*xtKl}nZHL>p1;3;_4$nEbFnmDi4A&cODBLKx32@WlX2C6hTL!lp
zZX?`wxP5TP;m*N@;W1ZdxZZF>;YPttfSV3C3vL12GPu=n8{xLY?c0`#4B$Bj7q%T4
z!1ab33O5RF0^D@CS#S&Bmcgxt+X%NEZXeunxO3a9rmpTE>d2RiUQdi*S2m*Q&nC<F
zU+)&auLcpZg{zcDUQe{x4TaKL+tXrCWy*eQdxo>`CgMj8dGe06i57pRTQC0!y+Ff?
zYtc9I^|gru@yXOBuE)qhUi!jTQe9>E&tb|D#NUn`n1{Qz9PxUzj9eEP&R-Uvhh~`U
zzAn-7A1;xxu4h<v6yC)xl@G4#5muc5Z@F9sFOR7xd|7T<*CU)S<(9q78zZy2ur4x;
z@8edyEZeT{5vGo;cv%jGcLoOKj8$^V`X1pNpSi0r4rS5$#IS0-oLl{>d=qi_P^;!u
zxo>^%5Qsn5_qKjrEks-Ccq3X4dn3wf#4a{Vj(wx2mAYGqY4Y(m9`mpzS-((KG3>w5
zlzzsJ0T&z7&;CKoCiOG!vAEc<e#XrZm;MEGsG3?(uz_)fs;LPNn^-ZYC}y$#9bPCu
zM&M#gQikF;+Pv77l%eV;Q23;R7WFgdA{fo;XS|?+i*4&?PEC9!yvUr(_#knSIp;?T
zbuI`4LIqyh#l@DU3{{ifza03P+h-JrR`xS~e}t=ch#d4*|GRa$4iSyt+Xc#XU=2wI
z>vA2UlNbU{wHqP?KbE@*%-6R$z8pw>CnR6m#)_#p2eMRdek;+l#I~NL%!PJ?U3&*s
zsb|4vDJ#I0c5AYfx54^y;+L==0{HEsF0~%a?-4a02J;(4)|OUXV0S&+=yY*g$P7~u
zlV3IZBeWC4!B5~Zy+4C?A`_n32Hysr1NK*|o!}Sv`364-o^SAMaG}A)4m@O7XegE=
zLaD*8ftMJ}C!b{o*Me6Xya!xi@VDSq2A=@0G5CxJib^PSjV;_FstxW0#xCe$crS3B
zk-_EQJqC{e*Bi_)UiTV&JGjAM&)rZQFcc4h8x4LOeAM9iV6VZ8!KV!7589i-s2BRU
zo#0P<PaFI(_~!t7#6Bo~2`I!#@ENe~(?7xI4Q`D)>IH+lfkmr8{Qh94!NbA$3F5^T
zoB)ox$ewoMAt+)mD%y!cu-o7YaI(Sgfm03s0-R>>QE<A!zk)LjjzG<_G<(oQZYU;d
z1vm|yZEzMi$Kbz%a}Ayi&NFxsIN#va;6j7>!L!t6jK6kbClqBygu~zpgU^E37#xB9
zNTtDv;A(^WgKG`uuczw_9u2Pdvque*Nl-Kx5oUlJ4W1458e9x+Hh2a2yuq)7g|UwC
z;iSvp9pEU%nEz}69yk+gL^uw18+;C&Y_Np|E7jmwaGJpb!RZE%1ZNt2Cpe3Y@|6Mi
zLov~aFdLk0@YCQNgP#ZI8vG78&){v~e1m!JPNBiSf=j`C$f*oyhr4r`5urP{!r;rn
zYYe^_TxswmaJ9krgKG_b99(bk)8Gade(0tR;2AiLMg%^;^%}en+-&d}@OgtfFr)Cm
zJm>RWm%(GeQ3hv&W5exZ88`-Vp>P`!*fhxoKMPJZ_(gD<!Eb@n4gLt6Y4Bcfmchrs
z6Fr9FKTu>F9Dxl{j=_A~k!x@taGt?w;CzE`0v8(mcW|k}Jdwv!W+;lGs4(~i@EU{v
z1+Fxh#~4)`>;=~v{0F$s;11Z))f?P10{8z0Lva}*G#Z=%_8L4H+-&fp;PVC-gGEGO
z47>z(8C(U9GI(1A=3lI#_!<%12LAv~HrNMFHMj#l?xh*r6P#}F72r&RZvtl-Jjnya
zL__f~INRW-z&Qq&gL4gj8=PnG7I417UxEt_J^?N@*u!0PnW2aUR~Xy}yvE=wz?BA%
z23H$A30!ON!{9oD3&8aO_J|cwGz1i)65MF;PO#VD{orPU&w<Yy9LbLXcqa^#NiA%b
zfL#WU0!LkB-2eX$MeIdIJ24CFHh2Lz+2Cd1RD;)n(+u7QPB-`jIMd(=d@RY*%=`ZU
zC?;wJnBR6}8=MW!F?bF**Wjh#JcHi>=Nr5STxjqKaH-Ae{vUx4Mtl&%XU5!y_Xd|C
z0&e@N32*m*Z~|d#g?&Ev_52#vzsKvU&dTHWubQ#17c+cL*EmS#$8G%k`2Jk)2{Ylr
zgE_MQ{?C1U4HUXS?&E6>=01L_!Q98!8C(qBV=%Yz^#*erzt>=H;~Ny)_y1Zb4j2)*
zk8d=X`}m^<b06<D_zbuitOp&p@#hWZHeTR!rw->f-UW8*`#<;bkyIcaZ${k5M;Xk0
zd`E-1kB>E&+xU2cxs7)l%x!#<!Q93tgFU)H?&DL82=}AU@nEjSm`F9k`NOBR`m<n2
z8}TYwe_|CO(PCo|?*#KX1L`LnA+h2x7*8YA{sW56!V3j|W9*O6QJe+yJx4#c6CrqR
zVx9%oA8ek2v<7UR`}76aJiaX&pYF}W%bo^r*fAh*;LNj7n5UVo0qcuT0*-=v%hQ_%
zgk%quJwEE+-8@Q-=WxyVGjPDr38s(_e{{{==1ELEvrb=(upDe2y0aFnA7Ahop-QlM
z%1<rWJP9fUPayCmVr`<wJj84#wB|uH{|1{!`CI^-M||Cf=N5domP3I@`}|Wsy96sM
z!GVM9+K9IP#~P|44Da0;)7U&&Z!IF|hZ;QC=3ijG!lus-Y5<#uo_!Clu^)o)M6GjR
z^QfN(@UYZ8CT<DXJP7YZ3{0Q)iyX9hjK@3-Z6cnhn1>=Q0h>n#hU2k{c~Vyb*gW;A
z2iQE%&I8t$8=hyt=FwzTUBz`RPWXw$x0`&1wh)oovOU$eo&GG?a$~{fp-K0H7dqs%
zANOx>o=kKN7Ut1%?#GE9^LW6k5^&ov&YqhA))zA`1e-^}Z3Ua>nhn6=HRkDmGr;Bn
zO+SJU;jVxkgXn=(9EaR+wenvqATCjpf1IR%anjnwgVsia%~QLkgUwUBW(B2xIfzey
z&4b!Tx;-sVV4C1=X~2C!{4Cf!@u?DQ9u*hb9plnC;prN%ISxmG`Ht%_6vOVj&jW?N
z*aAy}cyAD&1Dl8E_342-yOH5Yu)dhzB(Qm?)-14jZkp$FBY_(Uj)Be7_ab`g?EGEQ
z6>OeYI00;)^!5PQJi~HkQ2ND8k0<ra6-9Z(`k)9k;916Ld(B|;<h9$AkU>b`q_?lY
z<~eAcd$lw;8El^9mH`f&7KlSMKY$|NezM6v`79`dADF;6hb*jj%L1ka@kVeC@?#a5
z-;SX0BjAgN^La20y7fVgu$bvVF&u0j8=40;56vtDR~UW%Dp+6a^Y_5!DWJcD&9m?#
zlQ9VGF7Plx5)}F(r<3tz2zBJKVDlKUC&1>pq{ZNg#(9O+V0|$?4+m>74q#M=<oP3b
zQ~^Euwmd>|JXk-a=Q*W!fz2Z&cY*m!1>IubgZ0G%!!JRL;@mUb+|)ud2y7lPSqk=F
znEy=h5)|eknAW8{r#W!G;gw+X=(s1q`eOQpLA)AV0R!C#;~01kY#z3knqnW^it?3(
zH$b5;W;i*BR|N4^uz48XkbXL~ze10J%_Amj!R9%rzkqo}B|j}H1MchJGJ|<Ryf}zo
z1{;Uritqs~(~kn1N5?(`HjkA26pUj@Jr~a<y)-p&i|28IL&4_BuGwI|1<hjELTkYK
zB7Xulj|u%9Y#t#wY#{c-#sQ7fz*E2XIPsc?zXjfaLSM|VHi!>_%`+}f1chG+;$DL=
z>Wvd3v%uz=mkYrc?|&?C7ZvtJJ_I&TYW)f9Hr9#um$f|A^GUG2nErEc;J8417cv;1
z0F3htOTixVxZk}{m`6N@;!7a&{KDJ7`eFgof;b;+o&Z?~HV@D3I}{aCXV>~$eD+Y(
zpGUcJh2(~y0TC|XBY=P3*iM`P>x%`3UarYM2D*aHGcPX*3cnF-9vj*a6y84V8Z@b~
tl%}TPQ<!lAV<}i)G@uD=p4fb2)xq7L#Xom*sN+_dofcEox9O_*{|Ag!1up;q
delta 24507
zcmai+4O~=J`p56UTpbV<K?MY1L^N!D0bi-q&}3uF3{4YDT{Nq)tgxtPvmHz+d?C@}
z3;B{{){;??SW#KxONwQNW-hf-RGNb?yJT8a+WddNGxs=$Y4`KLd=Ah3KIiQ@&pG$p
zGk4~mm5c1h7u!psL(=!z9fJ+S*tWdK9V1g-FboD03M1SwY!Fg6$aJJ2qy-yB67+4X
z_#nu!NIw>FV+SqYH8kW{v0)r*Z#s?%o9)q<A3nAv{W{{wb?--9zv}j1qpOFX?e*yb
z<hfBSw=Hqsh_jmz#vzPHn1I0J7KDiilMp5&WFzpH;uG@!ZTLUc;^~kxt@FDr>4Kb%
zaG%8#_iN`kbsKXn#e<L!S?76{d<60_guhu#@i@X02n7f{*c2X5Av}#xM9_{!IA4sg
z#L_K=gnx_@gypE;Sb>w32&)mEMc}d4Cye!$q}YJ45upqL?lv|fyom4;LOH_A2<mu6
zomjF0@->7?gf|ejB2*!;etB%O<lB(%BD`nucF6Y;K0w%kP=l}wVK>5u2zwFsA@JCb
zZ~);Ggo6m5BGe)9IE3&S!V!ee5so78kb|tgs=NJ`pc!GIJ{|u@AY5$Mn%6j5^NBXC
z`Fue9vjO@4DGs$su<aK=;QtWOscnKgA|OYvkk%z`v$tkbn(Nci!U6th!-?~00qOS!
zq`#(h3kMqK0-A7Lu)joBFfky$K^%hFrvi$7IiMK_19FuGq^}H!pBT{W#I~)isBG;>
z@qnXKrwy%j#u2gHkud5)fHoz-AtM9gZ*J`nW4idkk&t>kz@DuE_ILtH=_^Ks46t<-
zb3+nFeG`!OY=B2s2GsbU0Dev!3Q0)47@%z*V9`lVde*^y+An&9CZsM2Xzj27Tb2g&
zL!6i!nqYfDEDud^KOIoRHvxTmZvZb3NdHklT~9_D#-r!$4!SDQXm8vfkCUJ54yYyH
z(gPRi1-m2F$8O^zOzCCKcE>y`K@AEl7-~8=Uv)CN8F$BF68^#NpdBpmr!d2K<h0$<
z+2V(i4Wrx7{t7O{7)GhxbkHZAjE+WHcf&{<rns9i3H42afr&D|yS;G&H=3DEb_e~*
z3J#*ed4_q};%B>Hat$_>#l|a$zqr#b{SslcGrA$;LB!L^%ztAqq(^y4D3|;$T6`B;
z9Bak@($g>+(@e*3a3?of_$p3JkH3Zg>1r79L8ha<mEiz#rlSQYPkQ8e7|@Igc^isA
z1;dcASabPnV3*|)#<zp3a^OMLC1<1IamxTVUAz&#SPcUbtO|eZWEe}mc1M9#!AQva
z&)6OCq-;@dbV2?le|%U6!*Ietx`gG8f{P0@m*3jkFlHe?2RA)(JsRM=N_xO8Gn_>M
zM^D)ueXRt;I~vB7n#Iq;vK$oXv?{oWKJE$^D}&zw8}OJiFbW=Qx#%}=A3U<xW;)O_
z(tx*NP#nr{&-)N9nur#2F@kvM(ha!Idcz~iWuxJtcQK0T38wd;4~J_lyns=@1$_*)
zj6Z}iF$5L(23dO}69zqrcxB*gQHJq%mXG>nh6=<S)w<woa0@D=MXc}$JTmy4-SMe4
zsKU|WpS1LAV8BXTMatmwxJ>8lQ4PXG0ZSI+s!<*YxdIJ_L2h_LR&Xl@*+jUI>mF9P
z5fu&p&|knKFev|byF(zJgUW?Lm4_}+*R>P+^;U&!^>hrfUs0aADh^=Gyrmhu4<0zt
zXm>kqunM#x!&VFuG*lYU;DpPxQN0i@ECrWa@hcJkF0KM)Ko?Zt!n}cKtR8>x{)HaY
z2I<d8xKQh|U(n#>VQ&AZZKMaZ%O(_~G#M^amtkxY`WP+BvJ7|ygEkEYp&n`Q1Q<|?
zE>nFTgp)Lsr#x^y(r0O1a0dDQ>c*(#5|b{cYum#`D#NcqhH(RUkQIL)29;g&NCie)
zp5}q`=+j2Da5&QQGW-b!glh)gh4@|Y@Q1;^`le$Hx%<OHbzQFMhMXuc&#Lfs7%&h9
zs6pjLpZA3al*><|4~L>f%*F<|QSgD?{xMUFG4vDS;aC};f`kz;V1mW-d*b@b*RIpa
zD4+uKKp!i?7pQr?HkbE+izi;^@8j`sRh%|R$Kfg{(>(AZTs{x^Q5*(s2O}F(`5AZ|
z?T`j_iN*C-qRsVlnGh{jE$k35sNTSs$We1QTwH+}sGrs)>ri1MJjgEORq|&PRE-Lu
zmK84Rh5F#}lejxEJ`3}&G0PgYP9w(X!0XW(gmmd*<V-<}*aY_Jd$>H(t#}r2JM!mX
zD&^9J@yp=Bp|~vjS`~herq0vS_jGar@9#o)U-aP|%v7~jhP0lTNwrZL0uP+W%=iN?
zPgXbuE$)m0)gYY*12$qJQDek~+jrl+{xNV8S4}0#%dpBHg+b@^d(drsiP0CQ%~Wq9
zVO#C0_y8`u8kf~HtAICPNV+x%cVdu^)COr}GH$=xO!z!Jwpfc_3WJtn5c(#*wDv|Z
zTvQF0+-y}aPV16zbV)~C;GJQBnkl=ZOTxTWM<cU>!-q|rFysDN58OBS`WX*Dlrd|T
zxc6e;s&P?o^{NUQvCt+;M%dju;2GQx4JShQ?Ae%2+8Re-BM%$?Hfix$y&nU$F}~8`
zmud07{HB&a$oSNXPeFMPCCtV+a2VesKkGz#EO3nARvDge<BOyLpS^xBI?`^OvohrM
zS>Ztc8<#8tIu9`&ZQu_s$V@mW3EUP}buuL!-<J|YE`<_fdN3tQOru0S!-rskZHpn1
zP6dh{LkSx*Dbe#|Dba!~N?zTRFgcqNj-E=19?YRclV(xEAQvU(!#R|=$Z{#sk$IH3
z%JV64$C&5l3C4H<C5FyIN{qoG$}JcQl$e)FC^1p3q{LNTN{OCdPx%R^6-rDmn<;S(
zlv6q}T~YSH7^hr-M;Vm3hN>ws)ZE*7vKXIWDKVD!P+~0AQes9vKnZu(QDWjdLV2ZO
z)KmTm4>Ty>z_db%%lbIw(}v-p?1%Y`vOiXYlml?FQerT)P^Myny~vY+=xp?1+d=4U
z%E1^WltVC<DTiWAP!7Xbr%W@9Xv(W_aZ+B50ZjQ3#sK9t7z32J*Yu?vfiXZi5@Ue!
zT8sh8bT<Y8Pex%3P>#kJpu7%afN~7R0A&Wo0Og-C1}Lw`7@*9=7@+(M#sK9FSiLVg
zam#`;b`J;e?elksKZJ)I;-c5<Ej*lOdg3uWJ+ER6dXiCBw)^(|yS+o=-~(?jC;l^&
zw^tI*Ts&=RL_Tk?WHRHL>Lm9i^9pLJk(@%t=rNiqB@ZS$$z_t$$jRgq$-~L$1f$7a
zAQkCUz^O)4uH-S~baIa5OmZeUOY&H97CBvV7CD=oB6%V?hwPM`O=gRl!X-~7=Yri$
zhE(KGk<SFpzeCY$7P)}jDA`3WBG*ZtLoOlLNX{jfk}D<Wk;}+slJm*s<PypA$d!E1
z(^Mc81yocsL9XP5<Qj60<RWq{IZN_VaveEcatXPfoFaK8xsmLYTuSzk!zHgLH=}-|
z$&iXNDq5JJ`4T;_nT!b-9*|s4=7+4NI>{B}aB_|0N^&&0Qt~#klUydbnw(57k-VLp
z;-;cNDr%@mBj-xqLry2>NUkMklCvZqAZL-&CD)O&$tjYLkaNgR$@OFxIb8BFvOAXw
zLn<1n$R{`Ra_M!PTtIG=>>(GC>m;8cmyl~DH<L@rm6Fes%gAMtTgc_)63G|AZdWB0
z1yaE&$W=|wmCUKgRYT5^%qhuLOU{zasmWDGPM6Fn%2iKJk<6*e)kt<q=9J~~ki#Wo
z>T+WcP+>>~=P*|bx%oGGfODA(_d9q%GUqgxgIp(>bDJxiTqButoGY4KDVcMf%SkSi
z%sJ1MOfHc;d?@aJt`sT?q=IvxD~+5hnRB5lotz_?bD}GgoF$oaqbrM?E}3(rE1R4m
znRBHphwPNhIn(7LhfAJ%Efu*`7*fGW)Rj+e{?+Tn)nKHaoxbRs1;?ASCQfj1?dsyQ
zh4Ht#xa#wq7V+1|b_v4@)OdK*k}lymsUCv|?)}mSiE-CQ#Zckxyeb_kU+h(1KKT7!
zGi;dn`}HwhGodc5$&27+Z{Q?kpFJsaka+d_K10f}EC&_99M4`@wr@<5o$JO$_9Rbl
zT(9pjOVi3|JlZIng3(L;bD~{l&+dy6>(^CCy3I;4UEG-2bL0)sosM7h*`mGvbaz3w
zsCm(sEN3z5D=t)d#>x_Iw@Uc^qNvDBbe{sAj`bCx61K65XnGp!iQ+<)<3C7RTsW5b
z(T3u}sa8o*zxaz#X?HUYRy4)QB2?0`(jZ0CFL~XUD_6=JRL=LXs-d$6n5Lt+kTVCH
z^s(i_;Zp0ev~S9CWg2Ok(l<!GqSrAHPDlINYU$jbbx_jM)8YO~o@PSJDXQX!taMFD
zw9B59hHmq`hNdsty{Mj+?P|gKw36b&i8M-@!Re1tKmUyKx-2@nXy2l^Bu^h^YFTuy
z>B}z0t?;i*h4f1@N|-nn<ycD42wcv}6jfXw=%jPkVgU$KFm?=2l4eS-c>acNQ|@b3
z>qeT#TC<R;xG>K$a!JNYR(dBET$U141-Ge6WkTpSLMJN??;<+inAvNp?0|f$g<qj|
zcGf*oTVQEx#rzwiLYiqAN_a^;e`AdMVWwD>#uO8i9L)1KXb|r`R9vVsE`qVDr<Cqt
zHo)pARqYI>vpP!YZnDxf{f^}_%lBDY<g+Q-)1j5AC6=izV*J>s5T|A8Q89OHjC+YR
zHAOSE#BZv~_zt~n4JxI35jth6Rj+3`(^;k}-4j;2rv9x=EmfwDL{t1;EwfBb7Cml?
z3dv$K;ng<M)D_ayWX;qyep6M(rl0(#D&6PMDO0U_Jv*7sGF9nbwbC_B$A=NWSId>D
z_oFF(Q!6b~v&9cEwSdipsYAu6n_}G0N>iPhsT=&Js*Ew3sY+*RrdsuSez@Q_Rp}a_
zlOy#Re17tqTCGfd4(0ewt+7lk5<AC5g;cT*n0lA!G%m)yQJNa9nfij?RF&~BG{G7G
zN;kwZfs=|=ucwD*s?xPrrhX8vn_8<(J&1U}sdbjA)nfV0Q6Y`Y3sYB%H*St`zbH*j
z*GzrQZ>q}p*m=LHN_Qu8Y5-XEdd4%IH2{=ugq5!8-}sd3_iDW|^>@VkO>MMH^@tf+
zQ6b@$se8rztQhy((o~0L>Q29@D&xyP`b|~3RnRF@t$IE4na(m*>F%@AHO03w)uT-9
z-^$cx%hYI*IzB2Soy~+-e-h)z$GG=PQ^Pe=>;0yxj9;Ado2qmlLZ?i%>h-+Abe5?~
zx4}x+G!D-IiVMA#Eql=d&>Z3z-z<58rWCeV3EqO4e#&PlFgg*prbQX+{Sj^?Q)5$Q
zy_YTIJUZN8-x}z|a}%7sn^DV>jP)qYdsx}3jGAm^Y$_7#CiJvj5<4dRsYk~){tUPk
zZR4$|4$mjlU2&XfcS}EKt70$T2v7Y0OL$BASY+SQ3kClDmL3Ds!%RoRw3_0=wp@jv
zQ7uSik>!dmsI$0`A9fh;5r=N+e8*pLS6hO{*ZLzW>2K@8REgJH{yyE!UqLN&`!Wu|
za-->+P?#>)Q>ywcR*t5BiOh*3ddqcO|9|_&dMpf~2_639W3{fkU%WDLsQV-ya#*+5
zOg#Tsl2OaFJ&+geTT*C8SMq{ky^I)*d%zc-hN($m2zL0Ak2^hyGI`k_lMni@$$vb{
z(#ya2cf(XMcv82n&ot9{ec`GllWK9pwX>!WF>6wi`^Xuu7f%;_u~b(fuksgjP`42}
zIard53mukj3DXVod8v_05u|%mrL(-J(*D)bHyLnjt46qGqi@6Aq`}NYpMSN#!OUXc
z2RPv0FtcG=Wvje0>$U%H)+(8`x*?;o0n4r~=%j`XB`~-lqq-q(WRhp)X=+NPMkZN6
zMIjZiW<!}un&kPCy|AI&pK^*Qx%E$}XQ3`G+zz#Gy~hvq^u`dWm%c#<V_jwGyNffo
zCXckkEahhv(~O$ce4$D?2A#SFm9F-*m@xT@m;nJ5Y_}?SRIHesJYqUCS_5Cj+@MuZ
zYNZ~kRiJc<;^)cTyH7-n92@>N*zqxvZTsSkNX@>YPt7TR9V)5<I@O3qtBz7JKRYsE
zp%v4#AJ40b3zJcnHXff6RoPcYUV}6X|J8BfFZj&XkWtfMriY2<>|yR~H7KC-T^M|<
z?zy1}T{5kvA)^u6&JAYHACmF=3s290WUfCXlcinsSf$nbO3S^ROk_`K-{<py=={sk
z$|imv<$28`j!lVid<XAwP>V}bZnEK%+HKbap)ILXV}c68#Mr4ZS^vX`H&3>1(*H>I
zj7oOdR*^NOOP{*~)c?cA^pn!Yhv98$<Gr^h*y6>K+lLG;Vxu;sVOnZ1edkW?JR0Zt
z^w+xQGQ|1YuZ;ZldtcAZ!1D&xbN<%4|4TDCvxOm--huxk*@+WYvdgWt=7`I~p^W|f
zK|nHp>pNcVBBon&LN@;ZJCt3-=07Y3PaAUQ|KEJJZ9@uul2K<(TE8@GNQT0{9?G>Y
zlT`GEe>s$IX=BB~X;HT8#P`#ZM*a?q)m5tEPeG@yAXU;)EuGTsIw3Nr$Bo!_++T#8
zg>ZeAL8ppPv-47>vsNBTmoJ{5p6OPrwtl#i;g%!sZZ{(Jl1T^D-D2!_zLs2N-R$}z
znk`c|yQS#I%k9r7;2iNXOupRsj3TO5Lsg7e&WAqD9%!a`mPx%#vXWMIk@po}iUm-}
zrnF~H6!M^!LMmpc6kX(F6zk~(Qerj#*T@lhFirU)3*E#ra-dUgSGuuUI{E1x>C!}}
zoXGY`Rt$3%iPW5--EF9{VOryoj4<kcgmaGF)OhTPdfs9tiG?{=xyRs(ThE`rfeY+O
z89e_8S*rOX1^C2K#^*us0%KUhCn%x8bcL9XFB%)VXF|)^V8g~w_WUh~=XPf2*Ufr%
z#xVOc>>Yd3SWKFpG<*@YXt#S&eS=xzo1mf_%u?T3I5LNM@C_CVxtG0-S_-W$y^U_H
zf1v<nqO3VA>ntW}uHhQYGG76G<=P|(*_;Jfv0T<iOT8?f`h8>ENsgsf`7c=IKi(*P
zZax3F^et+`vXAB;f}2n?f3HxS@isL#qd`2i)_Pk|L~%wXV=`bRdd@24ZLO3gboQ8}
za5T@@C$63u<1VA)8oy>yGVckLs>YU79Dza(jhn4{gXv{{$4&>SPGT{vx+^Lx&e*}z
zb=Ik6U$wNanlZ~+hwm!u>}yM<RS9E{G>QQ;<J{G}0^sCE8rVk~C^Prh43vuQnt`%n
zUWb1$kj1d{zrq8>8T)v8+&Z-k{6rd9!<a9>^~Na+eU15rPQ+&c#&)NHUB$Ay+PG^h
zBjv5}!LL}fEZ|WjS6-BgyR7`D;YC^bDk%P7B#U7U<roIV8K3gBpEUAQ)ywtL$Xdqq
z)Qyy7j@FF4>uY7C$erbI=Wt4YuYtwN%rz)h87LL6L!r7^R{Rwde=v~6u<jpZfBciD
zk4RhpsSNy18d%4eIjqCi&5^#waQG=NzAkp&+1*`dIk1fO*`<9lZwT{RQ=C+^pc-YL
ztoMrlW*>`T)mKaVzUJu;>(uJy^U^-~9P-UZt(Rq)N32Siucbk>yE9IFb*IDKNaJog
zM(4?l6Oc?9Cl#YK<7B0|Q2fDpEQTd7kmKYePaV>@lgfEvZFvnhGUj*0Ty~x;v%6;8
zD3O0xFE@WN=UM$P+9mU>V;;+CQn45c>9pI6Gh}`5K=B8=SPZK?0QYmUGx;>g`cxIo
z4E)uXezEiO7BVDw)***wq_6zjtO{A&-^Dd|j|qC`tE%~T=h$p_SAFa1U<<1KLY$kO
z)@CAl!81jq+;cGKqoXVC{Wj?PqgBW6U0?=1f23;m+^#{k*`nZ~D{ON``9rm~OffgF
zYu6riXbV1YdRjj95284+A#Y;P$xp<Od8=&cV$s7hgFg90oPBtFm$e`B-W;~5&Q@cq
z+u0aQi`s}=9vN(#FNz*{y#0lb=r(z6p880HKN@pY!a;e1D90@B`R)K>Cw<HZuN<Jo
z8QYn^J!&q__?V{y<!K$gQYmIWdQH&jgJR>OaWghBW5xj{c@pKR-jIs9P{`iEU0ROm
z9H#Jn{O~Dmwo)|B#spnlxY_z>kPTa{+h}_cmmlMfWG$DYFMt--<jrB<eG2Q*cXD}h
z7TKs%>yMj0@$p$L&3+aQH(4KePgzNV@XiF=uqB{~P{cQla^t?og2RwXRW`qML7S|C
zuCxl88BowJ*xWP;<rHV^ksiL>nq+hw3J*(V;czdO%Ca&DnHW~3K8(@|YdWs1Sy<Dt
za2rxE!yqdIK0gO$$Plw0>wnis&D<^fTl*#y_OMCHSHIE<dlG&WVR)WFTJ@AB5$W(P
zM#gc!`(H;i`q{eMb`mX*^%(W~J`RH;Y}t~G<M5>SGa_v7f&x^ZONE^!;?B}^8g*H{
zT=%i)o*#Y1QIuSqAus-q(PExjv)EKL#hQ$Eek`8-TjJni7_3H{oRj83p+=i5;uS50
zRIC;M`CD|i8Bi~1w$YPKm-JzkFM2%QXV$pARxz>u7GDjWbxC{;yHvIE8aW7ss#acI
z?`tWfqC%zU%D#t7)`{01kLkbx#EjIO*jM%W;}Jpb8PbhEw{CSaQ%tVqCE+NZ)U1Z?
zp;H{m;funHGh~WhNTDo~iq23-%OZ+1er6fhL9t+u^Y-b}cKg=SQrpi;h7>`b@l~$B
z$JvHuadFxat^UK}lrbzT;zQr~F}j+ZJI6D_IPspl+&*bHo*3>Zi0PlX6ZO~AMC)Rk
zjF{8i<iL?V_79YymV52QYXw&hXs!v!@^wI#o3$)eyF}ahE?Y=d?))5Ekf%l*esW1@
z_f)(WbSBR8<NM<71>?ig@b&A07H?kUdufBjo&{I5ElI%uKenst+=7?dCa7QQVjvnu
z`m&m<HcZ&`=dYJv>yyt{y}o3=&CNu%*x+jlFhxCpbkXm6y_r}8hkWMs7C_zt>A~-}
z+JEWwhJ_hMAILJKO@-VJITW%U@&?HBkhemH<JsO!$i9&GLuNv5f}92U8e{?FyO8Vc
z`1Qt@IH`u>SI8rf-EaeJhP)aw3<lf=xgPQ<$ZE(^$Rm&ykj;>FxRH0p2Y_E8Qy_ce
zCN>swETjwaHpqpLd5~q0YazEoz6M#}2EQITf|K)5bi;i-ye(P;*%z`1axCNuNdD-t
z5^^ErQOGjL?;*ECR$!TS46+)s1#%B$X9qlvb!!UbILNV(Qy}>hq<N4FAs0fHxp7j0
zlkJe%ScxBlOvY!gi;#mMqeJla8suQesgPNak3r6XTnAYM`4Z%2$TuKsAh$yvgPetr
z;(YwL7%~p>CCD_$-H;O@4??<gadHkPOCggyUT+0t8e}cxbjah7+adY;pAmf_6
z-gL+$$f=NNkok}!Ay-09faI$^Qy~vP&Vn3+v6By(1GyA(9%NaQ*X=FE$u=mqL9YKd
z`UEl#k1WrVA;a-w!MJ903FLanOvn^G8lMAM30XwOv+DEs&SN2DIKELTgG^?*ked<D
zC*3uW#~_bE{s~_*;Mc;&f*;*z5l)UlkqLPUau#Gb9{CqQCPA);91U3wnFV<Sau#GW
zWIkjVDn1CA3>l9ve8xZ~LFPb~LN0`?fGmUj3UWK-ImmiQy93|5je*6T#Tz9tk^Gr#
zY?$3%#NWyrPK5cQp(NJ4@2J;nif>9HMfh^Nc`GQAWac5oIFY(M);1FVqs$6$f1Zf~
zq>D$oK{DNO=!c5uk*>E)7l(UnSJ7!jlxSXVw{_#0IT5;eoW+Xs%cE^eMSH~Nth9@q
z6&-`(_K3MFBHPt&^LqI+VS7%zShk`cmR|3#h&A8(9GxSLYhGbD_ktp1r4z9n*e*JW
z*k0S9<I6JkpPi;-R0qfCLM%w8EU-U?`GO5yiIm&KcoA#cCGHck=3yA#OUxCKqD<J$
zW1#l3L_8fZ4vL*fcR(~C-BHAJmNkd3w3`QIF6kDes~4#&V{Jj=Ci#E%$^>)pKfT^i
znYs+A6G7d5rDC#45<8Kuzi7b!MDYt!*J3(qCsU`avYXY4%#}!8F2=8lwK>Fnt71(~
zTfTH<NY?~P^p%RKsafnqx?`dN{||^?RwdZ-M9k_0TauW&I!a`%wws+X_xP;XhV)&a
zkG0uF@#<L9M2sv+)U3AK+A&8+5@sn*bFw}&9n>T)t&R<?hprrcw1<ml)^;@u5HVGZ
zSQBgBg+gLwBXZZ+CHX2lf)sa)<!fSXw~ANsf4bPeCc%~~4wXi&z^hd*Wa}+#Ya%<=
zz^Vs9(^&f89^y}>F=mCb%-oFGZmXD48f!Zup1}Va79c87LN(IJ0c&z5+ys5BZLK(4
z8fzBG7?drX&)VDFvCr#$7slArdx!zg#+bigof#($&Uw}jgLbZuGWVc_XfglUSX(&D
zH#s#rWKDHQW8e&z#Tt9Z5?Ip$Pw=Prcu7P2hVr)KQ)Pkdk339>PLaAc*0vD;ajAek
zIForr!CHGeK4%bY=bA*ZY%TJFWnNLU)^0j6Ir~aV!i?w;-=HJ||B-*4V?qmJZ)V<#
zME?lC(Ko`K&~sYvi6?&9OjE6;`O5<HDf`pGOk=>J$tq2)xMy8Ab1ATcG@sM(-^KIm
zVr@m@{dKWs0b*odisp58^GT5J5}AnUHcxa~A8VTf<Dx>Aa{5K=B%f39uElepELr9r
z#2ps%*T=To4jtzb`wwh-7t@ZpGYrn>U~ItuZtMsVzQJxr;NIfv*gT{eEK)bHrkgel
zvQ>-uIIj>JHjIcomts5}Y_1PBRt1|luf>0Hc0-gK=bziD`NnQ8FpZ1yyxcV3H;uPU
z>c5u%Pn%|O8{?@q)NE~I?rUSz$n!VbQuA?J^9zUZUE7gwIE*b0o*#BlUl_vw>q5+z
zLyfACk@G{1U*-9xQ0hO&|M(Q+bg21DJ7Yn+k%!tDAGPE8xps_xF$^X?|6EjCoL^{<
z^BunP4?E(#;kl?0(68#ml&^F$--$4GcN)1o!dM)^^Q{ropOF85i!c{;HbmseZz7E^
zBYFO0XOxh-F<LghtGOb|cqS?c=|7YI*F~Ebqm9aF^P^~Ew><wm=7J5;U&fmMj5EBk
z=9)O8B+lFsXH>_TAM`gCbobaGevda-_cvB2nEy^NUQ00dCKx;9`N;(8H#_-%yVG3R
z!#LtJ&p3?}PB)>r2P1aK|6lenw<a3rdYFq7jRlE3e<hKz4T+0GAzJ#H$NL#2{o=Ow
zGv4gS^M-zmUe=#xyxpH|{-ysT_<w1DS({=kO9_1`#n_m_^PMS-J)L5%O*KTSxhvIp
zFV)0U@~>3UZ&UZG*=0S0+BTiEJ8}$Dh)sp=moZ`E{$RXi$(@$0vt*+s$Klduy7QKt
zjWWrPS+Xmp67rRn9B#>6OJY>Xc)lUQbbR#!(k<W7pn~smQ1YD%jBoiGDAo(|9it1s
zdmmUj-}honv#>i*hLQ7>@8znHJ_cSC^}Tk1FXZI`n~iAvKo&PQqX7T#;{d$DV;J~h
zpfPWU-NExYkP(LHx;dhA#RqnW*Xw-%O);>{dEnRWco%Qf<^dvWv+3p*2++%8^e(($
zg)kH0L4+p}Rv>Ib*ov?V;V{Ct2<H%j@B}>qArWC9!f1pE2s05LM0gTm1;QqTtq8ji
z4kLVva1J484_@Mnz{(^MVIab2gb4^U5gtT%5@7|xCWNgByATc|e2Z`nAqdZ#A`lW0
z1|p0`n1C=7;X#Ba5mq2<LfDG13*j)rw+QDDg7%I@0r=@?BEmp~(FhX|W+FU@@Fc<t
zgiQ!r5q2RQM)(%t96}I&=o*2Lh%j(p)mtz254M$v6Xgk^?8j#Opdh4==<srU2z~x5
z-gww5hP<3$vPUY!l$U#&?4~(l{>wc>*m3i&MTd%4k%DgtUlfO+?+Zuw9gP7XF1?%(
znvROnug8lKBSh~l@oo5q+PsD0hAoL9d~vySAzqOX4{u2@eS29vyQOEFoKB`AXSvw5
zC9zFTJWihzr*O($F;ZR-5w9e+;SX@rUl3R0ls~%7d_l~_soWal1+f@cZf6X&&cap^
zDfYY)9>Rea{USaEiH27a+RVa_aA&<Ff+`Z*@KKWMC2?g%?;wZ?5HI{g{H>z5x%5>G
zXYpc1PxH*jSlEhBDjs#yJsfkv1_wJ9AMg)lFkW%uK@U*|<CP>H^b=(;?`1gXEy`ei
z48ft^st%SH3*Ib&*QgJc7aKk@!EZk1MThHBz1tVe)jAHkk1`nVee<9LDTC$PP=-M_
z`j}T27{2r|el&rD4)rlQ9S6PYV_uzjW92gQDjoy>T~^>HMmXqY%3ygJ`z~-Fb0Ld+
zE!^v4yn=+IV5oTNwf-|yy*48f@9_ESwP99B2CI5)Mi*lUIM;GS5MIc;5zOCbb6&C`
z_id2;RW{~8$u7t`(dG4ofx|6hS<8HAr&_JIVcuB=)@xY@&a#WUUXPdG+o{6qpyoGp
zioa$yFdN0nFG3FZ#_r&vjmGswLVkzmi|=TR1b4$@Z(sh7MkY>THNFKL2lg4?(ZE;d
zZo{c5?ni<ojUNXmYrF*9SL0{FDH{I+JXm8sPE6Cd20UEj12#MYNY@l!BEcAqPl7Ww
zJ`c_Us}`G>+Osw80?yI67ucoc9|oSIv3&fRt10-EW}e2A!TB1`0MFC-esF=tkAW9z
zybxTZ@e1%#joo|~z61(141L_u;9J4VHQohYq480$(D*cXC0O<8@8D96+u=^QUgLOh
znV;Q8e<(Kl6~;(#Ias;)R&b@p4}hyRE&|tRTn4Vy_#JSa#)rZ6m)YIXI1WYQWkpBh
zXRt@(P&BPsV<)&p<1{wi;Ws!7?9liQaJa^ggQFF@;qv8BIF$nY5;$37J|Iicct1Ez
z<FCQ#8eafsY8;AHMV7`1;B1TW0#iq$KNLAyf)QYs#$&;`8czl1Ydi;Bpz-72B8``T
zOEg{wF7>fn4!W12DAN+~QD3>n`@oeN9|c!y{5`lv<DbB_8k?A~>NJi9*GtCrPY<L*
z(WoUD1NLZq8@O5Hd%!IkKMppunQ$@Kq48R9xW-$-(PY#w4fqfWr<UL=aI(ffgHtr_
zh<k9F#y!C48grAFOpR{>XK6eYoDJq9NNK=5P~>O{9s}oUycnFX@pIq;jrlONNaLg6
z5{=J;OEu;uNo65;flC>{r={guf^pzVjdQ@&8s~y*G-gX`HC_g;(|9AeUgOunjUm<!
z5gY?Mpzvr3j)0prJ^^mg_&2Z->c8$ovFveZ%%`{E8eaj9);JC9bZd$mp-9&F4seRb
z4}sG(UII?n_&IQ<#&3hOG_C_@YwQ8Lb2P;zC|nvxVBwalaZhl*#)H5G8fSouG@b@7
z(fDC-sm4n}asMyV6wf0;xyDuCN{v4RS8M!FaE-=i!L=HD!F3u(;A2|7#yz=1OQWW^
z3JE+K-w1Bjm=C;LG@c7K+WE)8(_n|j+|Mdp<2S(38t-;P;nWnLfs-{p2~N?NYvnYJ
zqrvGKUjfe4_$qLg#$&<R8oQ@Lk)tW*f?XOHf^#)q2hP{H5?rA1esGb-{{ok2%nhDO
z{p>cvu-Y&4D~v0_<r@DPT&eLL;A)K@0oQ2!6u4I7)!;gfUjf%&X59Zjgrf1XqN8yV
z?9n(BpFWy39sq99m|q<j_!^$m0q2n%utVc{;Bbvsfuj}w1sS$N;ZzFn0dTU$C&4Ki
zhvMT8AE00#1B1)(-r%$_G4#E7H^0L4t>n7N4fgomrDCq$!wl0g(aOfv`b>O;QY*hq
zoZO)?SL?IDSWU=w@rSN1jk#E#qcIliZX;JyaJ8PNF<0yP8gsQiPh+mu3pB0)FVy$|
zxJY9z)|YC0lB@X=O~KWADOmLkSL<aObG2TsF<0xAT7It9w`n{IT&*z|>)SP+%pX(K
zXbP^@_h`)3dacG>tsl^stMxjKxmrJ>@jd8iJT`09^YvQ1`xzv-!TOjS6x7~$5v)Fb
z8bOi9TLHWu%%=tDV#5{`ZF~(@pL%UUT@4SIF9iGII~y0k{F$JSI~qaw_A1}M4bx<!
zKNR{lN#$UDW2nPmeSfz|d{ozWH(LzWH>6nxcKs}Re-LlmD>6P95ab*t9{Heu4}C|q
z_wn&Q%<JC|>?Al{eE7k&L0P9on;r3^Qt*M9CUXm&E%@xO9xQOvph~d5r%w%7-}fm9
z&jyZ)J9i|E*0(I%53RoK%y(dY<DQFPeM7I=c%o1+OzhneA2jK__<l#cTi+CK3mz}1
zX9wJq=4~*4m!-D-sR!%(m7M~+tY;0}4C^PbzDv)2c*>~n@U{}H?_r0Xoq`SfdGS=u
zxC#0`X4!abq3>t360Gm`7lNl3`X*c5!TN?xiQt?cMDLyPiFw-QVO3D;TkwR%!Wl!w
zT|46=^v#Rzfx2&ySiQ4<Cw+&YYcb885ufi&;8Iy0e429zT7BETo4aFD*7lkz0IP%P
z*Ms#f+Umjjeq+~R^A~+rz5=kmaZ*qMCO5gK8wyLrq#K){is#NrcrFsa>0s=Dc6kG|
z9I(Ed*5APTZd!{2^1l_pzp?%L4r<xX);lBJ16BtOdcvv5XW$E9eJ7?pV0{a>s2&)?
z+RjVd%26L=x4EI<Z!(9YTwX~}fYreY%LCXGz@dq)H?_M4?9|H32CIYR<%0FiuH5*l
z!q7J~JE3K8BEv6WeOJ8ro~k&X2l|uILJkdf;XPn|2e$cOK4at!k9$iNwa)K0UV=jX
zu7-8jv`%2uf%T0G?MbLu+Y#+vu)YIZ6IkD$?8;t<w+?z>IM~0#8uogc0)@Vf->2Yo
z>**l7<hy_ZY`u|Q+aIhOSRJe&KY%|1Cu@`Czw~&FKjQ)uV81Uck)#U;UV%QuOr^%e
z1hBrH=Tfk~ndN$Lj@IYzfYm{RJ^|~yb9PAf?@{MwpE;q>_k>FWtAhnh1oJ1Z)RFH5
z>)U%h1J?Ho-3V^+-@uFmV0AG6_h4&dM7P`_Ev}D>_Ej(xtZ(f&8?2smbN|o>!1`8?
zUxE1pd*!j;!Rlaz-LJ$y>_^5&)NJsd!M?2@;eRSNGl6w5!<%4zbIO=MVdFt<hqZBF
zeM`4z!0KTBvH*S;oC5=$h~t?01S~hfWB=2DF@0O_s5%X-4i-2!fL{;bda%Ap-Cz2t
z+&+&k0qa{o9szTIQn$KpxgVCkcjOabb<lt{0sK+`zXjGdsWrL}P`P~u+y>USbln8T
zR*&w>yW)Kdg}z_twJHAlH@Eq_0lXTU(j}l=dMppjAEUAw8n_*-4)Ql(ef!Q1srba8
zZS9x^t~llP@0^$qg}#m7POv&y;E@162iEtV{4F3pe4xthGw>>KzIHdA1<nf#_HQg$
zCg$#s<baR|t2~D!`ktx5gRrjBuJm4DeSgi>V0AG639x_TzV=4?V64Tpodq|8^?i3e
zgK<OFH(QLtS10;Tf_L+)6zgCG`2oC=Y|V}Aq+?)x8`}{>G1JIhWPP2u49q=MIVCdv
znW1R^dhK@c8xp946-EtH<a26&u)h1`$bk6iV0{D6GXe3v(%?~T<`|O#g}xW!X0SSF
jz$LK0ck=YA9v^=i*P->Uh9dm<WN}Sebk!8kHF5t33T_Ic
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index a09db085a..b51e713f1 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -68,16 +68,17 @@ int main(int argc,char** argv)
myDetector -> BuildPhysicalEvent() ;
E = M2 -> GetEnergyDeposit();
-
- XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
- YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
+ XTarget = Init->GetICPositionX(0);
+ YTarget = Init->GetICPositionY(0);
+// XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
+// YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
TVector3 HitDirection = M2 -> GetPositionOfInteraction() - TVector3(XTarget,YTarget,0);
- BeamTheta = RandomEngine.Gaus( Init->GetICIncidentAngleTheta(0)*deg , 2*deg ) ;
- BeamPhi = RandomEngine.Gaus( Init->GetICIncidentAnglePhi(0)*deg , 2*deg ) ;
+// BeamTheta = RandomEngine.Gaus( Init->GetICIncidentAngleTheta(0)*deg , 2*deg ) ;
+// BeamPhi = RandomEngine.Gaus( Init->GetICIncidentAnglePhi(0)*deg , 2*deg ) ;
-// BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
+ BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
// Angle between beam and particle
@@ -98,7 +99,7 @@ int main(int argc,char** argv)
20*micrometer , // Target Thickness at 0 degree
ThetaMM2Surface );
- // E = E + ThinSi ;
+// E = E + ThinSi ;
E= He3StripAl.EvaluateInitialEnergy( E , // Energy of the detected particle
0.4*micrometer , // Target Thickness at 0 degree
diff --git a/NPLib/MUST2/Must2Array.cxx b/NPLib/MUST2/Must2Array.cxx
index b627d4524..0c740c710 100644
--- a/NPLib/MUST2/Must2Array.cxx
+++ b/NPLib/MUST2/Must2Array.cxx
@@ -239,9 +239,9 @@ void MUST2Array::ReadConfiguration(string Path)
//with angle method
else if ( check_Theta && check_Phi && check_R && check_beta )
{
- AddTelescope( R ,
- Theta ,
+ AddTelescope( Theta ,
Phi ,
+ R ,
beta_u ,
beta_v ,
beta_w );
@@ -504,14 +504,14 @@ void MUST2Array::AddTelescope( double theta ,
for( int j = 0 ; j < 128 ; j++ )
{
- X = C.X() + StripPitch * ( U.X()*i + V.X()*j ) ;
- Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j ) ;
- Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j ) ;
-
- lineX.push_back(X) ;
- lineY.push_back(Y) ;
- lineZ.push_back(Z) ;
-
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j ) ;
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j ) ;
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j ) ;
+
+ lineX.push_back(X) ;
+ lineY.push_back(Y) ;
+ lineZ.push_back(Z) ;
+
}
OneTelescopeStripPositionX.push_back(lineX) ;
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 0d7bb48ef..1c4d652cb 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -23,7 +23,6 @@
// STL
#include <vector>
-
int main(int argc, char** argv)
{
diff --git a/NPSimulation/include/EventGeneratorTransfert.hh b/NPSimulation/include/EventGeneratorTransfert.hh
index f66562a6f..09c3229bc 100644
--- a/NPSimulation/include/EventGeneratorTransfert.hh
+++ b/NPSimulation/include/EventGeneratorTransfert.hh
@@ -26,8 +26,10 @@
// C++ headers
#include <string>
-// NPTool headers
+// NPSimulation
#include "VEventGenerator.hh"
+
+// NPLib
#include "TInitialConditions.h"
// NPLib header
@@ -52,8 +54,8 @@ class EventGeneratorTransfert : public VEventGenerator
double BeamEnergy , // Beam Energy
double ExcitationEnergy , // Excitation of Heavy Nuclei
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
+ double SigmaX ,
+ double SigmaY ,
double SigmaThetaX ,
double SigmaPhiY ,
bool ShootLight ,
@@ -118,8 +120,8 @@ class EventGeneratorTransfert : public VEventGenerator
double BeamEnergy , // Beam Energy
double ExcitationEnergy , // Excitation of Heavy Nuclei
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
+ double SigmaX ,
+ double SigmaY ,
double SigmaThetaX ,
double SigmaPhiY ,
bool ShootLight ,
diff --git a/NPSimulation/include/EventGeneratorTransfertToResonance.hh b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
index 312a3239c..348ee5ed1 100644
--- a/NPSimulation/include/EventGeneratorTransfertToResonance.hh
+++ b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
@@ -16,17 +16,22 @@
* Decription: *
* This event Generator is used to simulated two body TransfertReaction. *
* A Phase Space calculation is then performed to decay the Heavy product. *
+ * The TGenPhaseSpace from ROOT is used to calculate a phase space decay *
+ * with flat distribution *
*---------------------------------------------------------------------------*
* Comment: *
- * This class is not yet operationnal. *
+ * *
* *
*****************************************************************************/
// C++ header
#include <string>
-// NPTool header
+// NPSimulation
#include "VEventGenerator.hh"
+// NPLib
+#include "TInitialConditions.h"
+
// NPLib header
#include "NPReaction.h"
@@ -37,118 +42,114 @@ using namespace NPL ;
class EventGeneratorTransfertToResonance : public VEventGenerator
{
-public: // Constructors and Destructors
-
- // Default constructor used to allocate memory
- EventGeneratorTransfertToResonance();
-
- // This is the constructor to be used
- EventGeneratorTransfertToResonance(string name1 , //Beam nuclei
- string name2 , //Target nuclei
- string name3 , //Product of reaction
- string name4 , //Product of reaction
- double BeamEnergy , //Beam Energy
- double ExcitationEnergy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path); //Path of the differentiel Cross Section
-
- // Default Destructor
- virtual ~EventGeneratorTransfertToResonance();
-
-public: // Inherit from VEventGenerator class
- void ReadConfiguration(string) ;
- void GenerateEvent(G4Event*, G4ParticleGun*) ;
-
- void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
- m_TargetX = X;
- m_TargetY = Y;
- m_TargetZ = Z;
- }
-
- void SetTargetThickness(double TargetThickness) {
- m_TargetThickness = TargetThickness ;
- }
-
- void SetTargetRadius(double TargetRadius) {
- m_TargetRadius = TargetRadius ;
- }
-
-private: // Particle Shoot Option
- bool m_ShootLight ;
- bool m_ShootHeavy ;
- bool m_ShootDecayProduct ;
-
-private: // TTree to store initial value of beam and reaction
- double m_InitialLightEnergy ;
- double m_InitialLightTheta ;
- double m_InitialBeamEnergy ;
- double m_InitialBeamTheta ;
- double m_InitialBeamX ;
- double m_InitialBeamY ;
- double m_InitialThetaCM ;
-
-private: // Beam Parameter
- double m_BeamEnergySpread ;
- double m_BeamFWHMX ;
- double m_BeamFWHMY ;
- double m_BeamEmmitanceTheta ;
- double m_BeamEmmitancePhi ;
-
-private: // Target Parameter
- double m_TargetThickness ;
- double m_TargetRadius ;
- double m_TargetX ;
- double m_TargetY ;
- double m_TargetZ ;
-
-private: // Reaction
- Reaction* m_Reaction ;
-
-private: // Resonance decay channel
- int m_ResonanceDecayZ ;
- int m_ResonanceDecayA ;
-
- //Other
- void Print() const ;
- void InitializeRootOutput() ;
- void ResonanceDecay(G4int parentZ ,
- G4int parentA ,
- G4double EnergyHeavy ,
- G4double ThetaHeavy ,
- G4double PhiHeavy ,
- G4double x0 ,
- G4double y0 ,
- G4double z0 ,
- G4Event* anEvent ,
- G4ParticleGun* particleGun);
-
-
- void SetEverything(string name1 , //Beam nuclei
- string name2 , //Target nuclei
- string name3 , //Product of reaction
- string name4 , //Product of reaction
- double BeamEnergy , //Beam Energy
- double ExcitationEnergy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path); //Path of the differentiel Cross Section
+ public: // Constructors and Destructors
+
+ // Default constructor used to allocate memory
+ EventGeneratorTransfertToResonance();
+
+ // This is the constructor to be used
+ EventGeneratorTransfertToResonance(string name1 , //Beam nuclei
+ string name2 , //Target nuclei
+ string name3 , //Product of reaction
+ string name4 , //Product of reaction
+ double BeamEnergy , //Beam Energy
+ double ExcitationEnergy , //Excitation of Heavy Nuclei
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path); //Path of the differentiel Cross Section
+
+ // Default Destructor
+ virtual ~EventGeneratorTransfertToResonance();
+
+ public: // Inherit from VEventGenerator class
+ void ReadConfiguration(string) ;
+ void GenerateEvent(G4Event*, G4ParticleGun*) ;
+
+ void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
+ m_TargetX = X;
+ m_TargetY = Y;
+ m_TargetZ = Z;
+ }
+
+ void SetTargetThickness(double TargetThickness) {
+ m_TargetThickness = TargetThickness ;
+ }
+
+ void SetTargetRadius(double TargetRadius) {
+ m_TargetRadius = TargetRadius ;
+ }
+
+ private: // Particle Shoot Option
+ bool m_ShootLight ;
+ bool m_ShootHeavy ;
+ bool m_ShootDecayProduct ;
+
+ private: // TTree to store initial value of beam and reaction
+ TInitialConditions* m_InitConditions;
+
+ private: // Beam Parameter
+ double m_BeamEnergy ;
+ double m_BeamEnergySpread ;
+ double m_SigmaX ;
+ double m_SigmaY ;
+ double m_SigmaThetaX ;
+ double m_SigmaPhiY ;
+
+ private: // Target Parameter
+ double m_TargetThickness ;
+ double m_TargetRadius ;
+ double m_TargetX ;
+ double m_TargetY ;
+ double m_TargetZ ;
+
+ private: // Reaction
+ Reaction* m_Reaction ;
+
+ private: // Resonance decay channel
+ int m_ResonanceDecayZ ;
+ int m_ResonanceDecayA ;
+
+ //Other
+ void Print() const ;
+ void InitializeRootOutput() ;
+ void ResonanceDecay( G4double EnergyHeavy ,
+ G4double ThetaHeavy ,
+ G4double PhiHeavy ,
+ G4double x0 ,
+ G4double y0 ,
+ G4double z0 ,
+ G4Event* anEvent ,
+ G4ParticleGun* particleGun);
+
+ // This method return a random Vector of dimension N and magnitude R
+ // The return distribution populate uniformely the surface of the N-Sphere of radius R
+ vector<double> PhaseSpaceUniformGenerator( int N , double R);
+
+ void SetEverything(string name1 , //Beam nuclei
+ string name2 , //Target nuclei
+ string name3 , //Product of reaction
+ string name4 , //Product of reaction
+ double BeamEnergy , //Beam Energy
+ double ExcitationEnergy , //Excitation of Heavy Nuclei
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path); //Path of the differentiel Cross Section
};
#endif
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 12a9205e5..569c489fc 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -464,9 +464,6 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
particleGun->GeneratePrimaryVertex(anEvent);
}
if (m_ShootHeavy) { // Case of recoil particle
- //
- // The case of recoil particle has not been tested with the new version of the event generator
- //
// Particle type
particleGun->SetParticleDefinition(HeavyName);
// Particle energy
@@ -484,10 +481,6 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
// write angles in ROOT file
m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
- // tests
-// G4cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXX" << G4endl;
-// G4cout << "cinematique dans ref world : " << G4endl;
-// G4cout << "\t" << momentum_kine_world << G4endl;
//Set the gun to shoot
particleGun->SetParticleMomentumDirection(momentum_kine_world);
//Shoot the light particle
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index 0a1fb9e26..df57cbdde 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -14,9 +14,11 @@
* Decription: *
* This event Generator is used to simulated two body TransfertReaction. *
* A Phase Space calculation is then performed to decay the Heavy product. *
+ * The TGenPhaseSpace from ROOT is used to calculate a phase space decay *
+ * with flat distribution *
*---------------------------------------------------------------------------*
* Comment: *
- * This class is not yet operationnal. *
+ * *
* *
*****************************************************************************/
@@ -30,6 +32,7 @@
// G4 headers
#include "G4ParticleTable.hh"
#include "G4ParticleGun.hh"
+#include "G4RotationMatrix.hh"
// G4 headers including CLHEP headers
// for generating random numbers
@@ -39,6 +42,9 @@
#include "EventGeneratorTransfertToResonance.hh"
#include "RootOutput.h"
+//Root Headers
+#include "TGenPhaseSpace.h"
+
using namespace std;
using namespace CLHEP;
@@ -48,21 +54,21 @@ using namespace CLHEP;
EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance()
{
//------------- Default Constructor -------------
-
+ m_InitConditions = new TInitialConditions() ;
m_Reaction = new Reaction() ;
- m_BeamFWHMX = 0 ;
- m_BeamFWHMY = 0 ;
- m_BeamEmmitanceTheta = 0 ;
- m_BeamEmmitancePhi = 0 ;
- m_ResonanceDecayZ = 0 ;
- m_ResonanceDecayA = 0 ;
+ m_SigmaX = 0 ;
+ m_SigmaY = 0 ;
+ m_SigmaThetaX = 0 ;
+ m_SigmaPhiY = 0 ;
+ m_ResonanceDecayZ = 0 ;
+ m_ResonanceDecayA = 0 ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance()
{
//------------- Default Destructor ------------
-
+ delete m_InitConditions;
delete m_Reaction ;
}
@@ -74,10 +80,10 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string
double BeamEnergy ,
double ExcitationEnergy ,
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
int ResonanceDecayZ ,
int ResonanceDecayA ,
bool ShootLight ,
@@ -94,10 +100,10 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string
BeamEnergy , //Beam Energy
ExcitationEnergy , //Excitation of Heavy Nuclei
BeamEnergySpread ,
- BeamFWHMX ,
- BeamFWHMY ,
- BeamEmmitanceTheta ,
- BeamEmmitancePhi ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY ,
ResonanceDecayZ ,
ResonanceDecayA ,
ShootLight ,
@@ -113,13 +119,7 @@ void EventGeneratorTransfertToResonance::InitializeRootOutput()
{
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("ThetaCM" , &m_InitialThetaCM , "ThetaCM/D");
- pTree->Branch("InitialLightEnergy", &m_InitialLightEnergy , "InitialLightEnergy/D");
- pTree->Branch("InitialLightTheta" , &m_InitialLightTheta , "InitialLightTheta/D");
- pTree->Branch("InitialBeamEnergy" , &m_InitialBeamEnergy , "InitialBeamEnergy/D");
- pTree->Branch("InitialBeamTheta" , &m_InitialBeamTheta , "InitialBeamTheta/D");
- pTree->Branch("InitialBeamX" , &m_InitialBeamX , "InitialBeamX/D");
- pTree->Branch("InitialBeamY" , &m_InitialBeamY , "InitialBeamY/D");
+ pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
}
@@ -140,7 +140,7 @@ void EventGeneratorTransfertToResonance::ReadConfiguration(string Path)
////////Reaction Setting needs///////
string Beam, Target, Heavy, Light, CrossSectionPath ;
- G4double BeamEnergy = 0 , ExcitationEnergy = 0 , BeamEnergySpread = 0 , BeamFWHMX = 0 , BeamFWHMY = 0 , BeamEmmitanceTheta = 0 , BeamEmmitancePhi=0, ResonanceDecayZ = 0 , ResonanceDecayA = 0 ;
+ G4double BeamEnergy = 0 , ExcitationEnergy = 0 , BeamEnergySpread = 0 , SigmaX = 0 , SigmaY = 0 , SigmaThetaX = 0 , SigmaPhiY=0, ResonanceDecayZ = 0 , ResonanceDecayA = 0 ;
bool ShootLight = false ;
bool ShootHeavy = false ;
bool ShootDecayProduct = false ;
@@ -237,32 +237,32 @@ while(ReadingStatus){
G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "BeamFWHMX=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
check_FWHMX = true ;
ReactionFile >> DataBuffer;
- BeamFWHMX = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM X " << BeamFWHMX << " mm" << G4endl;
+ SigmaX = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "BeamFWHMY=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
check_FWHMY = true ;
ReactionFile >> DataBuffer;
- BeamFWHMY = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM Y " << BeamFWHMX << " mm" << G4endl;
+ SigmaY = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 19, "BeamEmmitanceTheta=") == 0) {
+ else if (DataBuffer.compare(0, 19, "SigmaThetaX=") == 0) {
check_EmmitanceTheta = true ;
ReactionFile >> DataBuffer;
- BeamEmmitanceTheta = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Theta " << BeamEmmitanceTheta / deg << " deg" << G4endl;
+ SigmaThetaX = atof(DataBuffer.c_str()) * rad;
+ G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
}
- else if (DataBuffer.compare(0, 17, "BeamEmmitancePhi=") == 0) {
+ else if (DataBuffer.compare(0, 17, "SigmaPhiY=") == 0) {
check_EmmitancePhi = true ;
ReactionFile >> DataBuffer;
- BeamEmmitancePhi = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Phi " << BeamEmmitancePhi / deg << " deg" << G4endl;
+ SigmaPhiY = atof(DataBuffer.c_str()) * rad;
+ G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
}
else if (DataBuffer.compare(0, 17, "ResonanceDecayZ=") == 0) {
@@ -335,10 +335,10 @@ while(ReadingStatus){
BeamEnergy ,
ExcitationEnergy ,
BeamEnergySpread ,
- BeamFWHMX ,
- BeamFWHMY ,
- BeamEmmitanceTheta ,
- BeamEmmitancePhi ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY ,
ResonanceDecayZ ,
ResonanceDecayA ,
ShootLight ,
@@ -353,10 +353,12 @@ while(ReadingStatus){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
{
+ // Clear contents of Precedent event (now stored in ROOTOutput)
+ m_InitConditions->Clear();
+
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
-
G4int LightZ = m_Reaction->GetNucleus3()->GetZ() ;
G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
@@ -366,380 +368,346 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
- //////////////////////////////////////////////////
- /////Choose ThetaCM following Cross Section //////
- //////////////////////////////////////////////////
-
- //Calling RandGeneral fonction, which taking a double array as argument
- CLHEP::RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize() ) ;
- G4double ThetaCM = CrossSectionShoot.shoot() * (180 * deg) ;
- //G4double ThetaCM = 0;
- //Set the Theta angle of reaction
- m_InitialThetaCM = ThetaCM;
-
- //////////////////////////////////////////////////
- ////Momentum and direction set with kinematics////
- //////////////////////////////////////////////////
-
- //Variable where to store results.
- G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
- //Compute Kinematic using previously defined ThetaCM
- m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
- //to write thetaCM in the root file
- m_InitialLightTheta = ThetaLight / deg ;
- m_InitialLightEnergy = EnergyLight / MeV ;
-
- // Vertex position and beam angle
- // 11Li Beam@Riken
- G4double x0 = 1000 * cm ;
- G4double y0 = 1000 * cm ;
- //shoot inside the target.
+ G4ParticleDefinition* HeavyName
+ = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
+ // Vertex position and beam angle inte world frame
+ G4double x0 = 1000 * cm ;
+ G4double y0 = 1000 * cm ;
+ G4double Beam_thetaX = 0 ;
+ G4double Beam_phiY = 0 ;
+
+ //shoot inside the target with correlated angle
if (m_TargetRadius != 0) {
while (sqrt(x0*x0 + y0*y0) > m_TargetRadius) {
- x0 = RandGauss::shoot(0, m_BeamFWHMX / 2.35) ;
- y0 = RandGauss::shoot(0, m_BeamFWHMY / 2.35) ;
+ RandomGaussian2D(0, 0, m_SigmaX, m_SigmaThetaX, x0, Beam_thetaX);
+ RandomGaussian2D(0, 0, m_SigmaY, m_SigmaPhiY , y0, Beam_phiY );
}
}
-
else {
- x0 = 0 ;
- y0 = 0 ;
+ RandomGaussian2D(0,0,0,m_SigmaThetaX,x0,Beam_thetaX);
+ RandomGaussian2D(0,0,0,m_SigmaPhiY ,y0,Beam_phiY );
}
-
- // FIXME a changer pour prendre en compte correctement l'emmitance
- G4double Beam_theta = RandGauss::shoot(0, m_BeamEmmitanceTheta) ;
- // must shoot inside the target.
- G4double z0 = (-m_TargetThickness / 2 + CLHEP::RandFlat::shoot() * m_TargetThickness) ;
-
- // Move to the target
- x0 += m_TargetX ;
- y0 += m_TargetY ;
- z0 += m_TargetZ ;
+ // write emittance angles to ROOT file
+ m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
+ m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
- // Store initial value
- m_InitialBeamX = x0 ;
- m_InitialBeamY = y0 ;
- m_InitialBeamTheta = Beam_theta ;
+ // Calculate Angle in spherical coordinate, passing by the direction vector dir
+ G4double Xdir = cos(pi/2. - Beam_thetaX);
+ G4double Ydir = cos(pi/2. - Beam_phiY );
+ G4double Zdir = sin(pi/2. - Beam_thetaX) + sin(pi/2. - Beam_phiY);
+
+ G4double Beam_theta = acos(Zdir / sqrt(Xdir*Xdir + Ydir*Ydir + Zdir*Zdir)) * rad;
+ G4double Beam_phi = atan2(Ydir, Xdir) * rad;
+ if (Beam_phi < 0) Beam_phi += 2*pi;
+
+ // write angles to ROOT file
+ m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
+ m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+ G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ cos(Beam_theta) * sin(Beam_phi),
+ -sin(Beam_theta));
+ G4ThreeVector col2(-sin(Beam_phi),
+ cos(Beam_phi),
+ 0);
+ G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ sin(Beam_theta) * sin(Beam_phi),
+ cos(Beam_theta));
+ G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+ // Beam incident energy
+ G4double NominalBeamEnergy = m_BeamEnergy;
+ G4double IncidentBeamEnergy = RandGauss::shoot(NominalBeamEnergy, m_BeamEnergySpread / 2.35);
+ m_Reaction->SetBeamEnergy(IncidentBeamEnergy);
+ m_InitConditions->SetICIncidentEnergy(IncidentBeamEnergy / MeV);
+ // Angles
+ RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
+ G4double ThetaCM = CrossSectionShoot.shoot() * (180*deg);
+ G4double phi = RandFlat::shoot() * 2*pi;
+ // write angles to ROOT file
+ m_InitConditions->SetICEmittedAngleThetaCM(ThetaCM / deg);
+ m_InitConditions->SetICEmittedAnglePhiIncidentFrame(phi / deg);
//////////////////////////////////////////////////
- /////Now define everything for light particle/////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
//////////////////////////////////////////////////
- particleGun->SetParticleDefinition(LightName);
-
- G4double theta = ThetaLight + Beam_theta ;
- G4double phi = CLHEP::RandFlat::shoot() * 2 * pi ;
- G4double particle_energy = EnergyLight ;
-
- // Direction of particle, energy and laboratory angle
- G4double momentum_x = sin(theta) * cos(phi) ;
- G4double momentum_y = sin(theta) * sin(phi) ;
- G4double momentum_z = cos(theta) ;
+ // Variable where to store results
+ G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
+ // Set the Theta angle of reaction
+ m_Reaction->SetThetaCM(ThetaCM);
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi),
+ sin(ThetaLight) * sin(phi),
+ cos(ThetaLight));
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi),
+ sin(ThetaHeavy) * sin(phi),
+ cos(ThetaHeavy));
+
+ // write angles/energy to ROOT file
+ m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(ThetaLight / deg);
+ m_InitConditions->SetICEmittedEnergy(EnergyLight/MeV);
+ //must shoot inside the target.
+ G4double z0 = (-m_TargetThickness / 2 + RandFlat::shoot() * m_TargetThickness);
- //Set the gun to shoot
- particleGun->SetParticleMomentumDirection(G4ThreeVector(momentum_x, momentum_y, momentum_z)) ;
- particleGun->SetParticleEnergy(particle_energy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Move to the target
+ x0 += m_TargetX ;
+ y0 += m_TargetY ;
+ z0 += m_TargetZ ;
- //Shoot the light particle
- if (m_ShootLight) particleGun->GeneratePrimaryVertex(anEvent) ;
+ // write vertex position to ROOT file
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
//////////////////////////////////////////////////
- /////Now define everything for heavy particle/////
+ ///////// Set up everything for shooting /////////
//////////////////////////////////////////////////
-
- theta = ThetaHeavy + Beam_theta ;
- particle_energy = EnergyHeavy ;
-
- if (m_ShootHeavy) ResonanceDecay( HeavyZ ,
- HeavyA ,
- EnergyHeavy ,
- theta ,
- phi ,
- x0 ,
- y0 ,
- z0 ,
- anEvent ,
- particleGun );
+ if (m_ShootLight) { // Case of light particle
+ // Particle type
+ particleGun->SetParticleDefinition(LightName);
+ // Particle energy
+ particleGun->SetParticleEnergy(EnergyLight);
+ // Particle vertex position
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineLight_beam;
+ // get theta and phi in the world frame
+ G4double theta_world = momentum_kine_world.theta();
+ G4double phi_world = momentum_kine_world.phi();
+ if (phi_world < 1e-6) phi_world += 2*pi;
+ // write angles in ROOT file
+ m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
+ m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
+ //Set the gun to shoot
+ particleGun->SetParticleMomentumDirection(momentum_kine_world);
+ //Shoot the light particle
+ particleGun->GeneratePrimaryVertex(anEvent);
+ }
+
+ if (m_ShootHeavy) { // Case of recoil particle
+ // Particle type
+ particleGun->SetParticleDefinition(HeavyName);
+ // Particle energy
+ particleGun->SetParticleEnergy(EnergyHeavy);
+ // Particle vertex position
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineHeavy_beam;
+ // get theta and phi in the world frame
+ G4double theta_world = momentum_kine_world.theta();
+ G4double phi_world = momentum_kine_world.phi();
+ if (phi_world < 1e-6) phi_world += 2*pi;
+
+ EventGeneratorTransfertToResonance::ResonanceDecay( EnergyHeavy ,
+ theta_world ,
+ phi_world ,
+ x0 ,
+ y0 ,
+ z0 ,
+ anEvent ,
+ particleGun );
+
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::ResonanceDecay(G4int parentZ ,
- G4int parentA ,
- G4double EnergyHeavy ,
- G4double ThetaHeavy ,
- G4double PhiHeavy ,
- G4double x0 ,
- G4double y0 ,
- G4double z0 ,
- G4Event* anEvent ,
- G4ParticleGun* particleGun)
+void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy ,
+ G4double ThetaHeavy ,
+ G4double PhiHeavy ,
+ G4double x0 ,
+ G4double y0 ,
+ G4double z0 ,
+ G4Event* anEvent ,
+ G4ParticleGun* particleGun )
{
- //FIXME
+ G4double parentZ = m_Reaction->GetNucleus4()->GetZ() ;
+ G4double parentA = m_Reaction->GetNucleus4()->GetA() ;
+
G4int NumberOfNeutrons = (parentA - parentZ) - (m_ResonanceDecayA - m_ResonanceDecayZ) ;
G4int NumberOfProtons = parentZ - m_ResonanceDecayZ ;
G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;
- /*G4int NumberOfNeutrons = 1 ;
- G4int NumberOfProtons = 0 ;
- G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;*/
-
if (NumberOfNeutrons < 0 || NumberOfProtons < 0) {
G4cout << "Error input for Resonance decay" << G4endl;
return;
}
-
+
else {
- //Obtain Mass of daughter Nuclei
- G4ParticleDefinition* parent
- = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation()) ;
- G4ParticleDefinition* daughter
- = G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
- G4ParticleDefinition* neutron
- = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
- G4ParticleDefinition* proton
- = G4ParticleTable::GetParticleTable()->FindParticle("proton");
-
- G4double M = parent -> GetPDGMass() ;
- G4double md = daughter -> GetPDGMass() ;
-
- G4double mn = neutron -> GetPDGMass() ;
- G4double mp = proton -> GetPDGMass() ;
-
- // FIXME Hexaneutron
- /* mn = 6*(neutron-> GetPDGMass())-100*keV ;
- G4double E = (M*M - mn*mn + md*md) / (2*M) ;
- G4double p = sqrt(( M*M - (md + mn)*(md + mn) ) * ( M*M - (md - mn)*(md - mn) )) / (2*M) ;*/
-
- // G4double cosThetaCM = CLHEP::RandFlat::shoot()*2-1 ;
- // G4double ThetaCM = acos(cosThetaCM) ;
- // G4double PhiCM = CLHEP::RandFlat::shoot() * 2 * pi ;
-
- // FIXME Hexaneutron
- /*
- G4double DaughterMomentumX = sin(ThetaCM)*cos(PhiCM) ;
- G4double DaughterMomentumY = sin(ThetaCM)*sin(PhiCM) ;
- G4double DaughterMomentumZ = cos(ThetaCM) ;*/
-
- G4double Q = m_Reaction->GetExcitation() + M - md - mn * NumberOfNeutrons - mp * NumberOfProtons ;
-
- vector<G4double> DecayProductsMomentumCM ;
- vector<G4double> DecayProductsMomentumXCM ;
- vector<G4double> DecayProductsMomentumYCM ;
- vector<G4double> DecayProductsMomentumZCM ;
- vector<G4double> DecayProductsThetaCM ;
- vector<G4double> DecayProductsPhiCM ;
-
- G4double AvaibleEnergy = Q ;
- //Calculate Energy of Heavy
- G4double MomentumCMHeavy = CLHEP::RandFlat::shoot() * AvaibleEnergy ;
- AvaibleEnergy = AvaibleEnergy - MomentumCMHeavy ;
-
- G4double MomentumCM = 0 ;
- // Calculate DecayProducts Properties
- for (G4int i = 0 ; i < NumberOfDecayProducts ; i++) {
- if (i != NumberOfDecayProducts - 1) MomentumCM = CLHEP::RandFlat::shoot() * AvaibleEnergy ;
- else MomentumCM = AvaibleEnergy ;
-
- AvaibleEnergy = AvaibleEnergy - MomentumCM ;
- DecayProductsMomentumCM.push_back(MomentumCM) ;
-
- G4double cosThetaCM = CLHEP::RandFlat::shoot() * 2 - 1 ;
- G4double ThetaCM = acos(cosThetaCM) ;
- DecayProductsThetaCM.push_back(ThetaCM) ;
-
- G4double PhiCM = CLHEP::RandFlat::shoot() * 2 * pi ;
- DecayProductsPhiCM.push_back(PhiCM) ;
-
- DecayProductsMomentumXCM.push_back(sin(ThetaCM)*cos(PhiCM)) ;
- DecayProductsMomentumYCM.push_back(sin(ThetaCM)*sin(PhiCM)) ;
- DecayProductsMomentumZCM.push_back(cos(ThetaCM)) ;
- }
-
- // Applied conservation of Momentum to daughter nuclei
- G4double DaughterMomentumX = 0 ;
- G4double DaughterMomentumY = 0 ;
- G4double DaughterMomentumZ = 0 ;
-
- for (G4int i = 0 ; i < NumberOfDecayProducts ; i++) {
- DaughterMomentumX = (DaughterMomentumX - DecayProductsMomentumCM[i] * DecayProductsMomentumXCM[i]) ;
- DaughterMomentumY = (DaughterMomentumY - DecayProductsMomentumCM[i] * DecayProductsMomentumYCM[i]) ;
- DaughterMomentumZ = (DaughterMomentumZ - DecayProductsMomentumCM[i] * DecayProductsMomentumZCM[i]) ;
- }
-
- // Daughter MomentumCM
- G4double p =
- sqrt(DaughterMomentumX * DaughterMomentumX + DaughterMomentumY * DaughterMomentumY + DaughterMomentumZ * DaughterMomentumZ) ;
-
- // Daughter EnergyCm
- G4double E = sqrt(p * p + md * md) ;
- DaughterMomentumX = DaughterMomentumX / p ;
- DaughterMomentumY = DaughterMomentumY / p ;
- DaughterMomentumZ = DaughterMomentumZ / p ;
-
- // CM to Lab //
-
- // Initial Lab to CM Momentum
- G4double InitialE = sqrt(EnergyHeavy * EnergyHeavy + M * M) ;
- G4double InitialMomentum = EnergyHeavy ;
- G4double InitialMomentumX = sin(ThetaHeavy) * cos(PhiHeavy) ;
- G4double InitialMomentumY = sin(ThetaHeavy) * sin(PhiHeavy) ;
- G4double InitialMomentumZ = cos(ThetaHeavy) ;
-
- // Beta and Gamma CM to Lab
- /* G4double betaX = (DaughterMomentumX*p - InitialMomentumX*InitialMomentum)/E ;
- G4double betaY = (DaughterMomentumY*p - InitialMomentumY*InitialMomentum)/E ;
- G4double betaZ = (DaughterMomentumZ*p - InitialMomentumZ*InitialMomentum)/E ;
- G4double beta = sqrt (betaX*betaX + betaY*betaY + betaZ*betaZ ) ;
- G4double beta2 = beta*beta ;
- G4double gamma = 1 / ( sqrt(1 - beta2 ) ) ;*/
-
- G4double betaX = -(InitialMomentumX * InitialMomentum) / InitialE ;
- G4double betaY = -(InitialMomentumY * InitialMomentum) / InitialE ;
- G4double betaZ = -(InitialMomentumZ * InitialMomentum) / InitialE ;
- G4double beta = sqrt(betaX * betaX + betaY * betaY + betaZ * betaZ) ;
- G4double beta2 = beta * beta ;
- G4double gamma = 1 / (sqrt(1 - beta2)) ;
- // Calculate everything for heavy particule
-
- /* G4double NewEnergy =
- gamma*E
- - betaX*gamma*DaughterMomentumX*p
- - betaY*gamma*DaughterMomentumY*p
- - betaZ*gamma*DaughterMomentumZ*p;*/
-
- G4double NewMomentumX =
- - betaX * gamma * E
- + DaughterMomentumX * p + (gamma - 1) * (betaX * betaX) / (beta2) * DaughterMomentumX * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DaughterMomentumY * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DaughterMomentumX * p
- + DaughterMomentumY * p + (gamma - 1) * (betaY * betaY) / (beta2) * DaughterMomentumY * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DaughterMomentumX * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DaughterMomentumY * p
- + DaughterMomentumZ * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
-
- //Set the gun to shoot
- particleGun->SetParticleDefinition(daughter) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
-
- // Shoot the Daugter
- particleGun->GeneratePrimaryVertex(anEvent) ;
-
- if (m_ShootDecayProduct) {
- G4int jj = 0 ;
- for (; jj < NumberOfNeutrons ; jj++) {
- p = sqrt(DecayProductsMomentumXCM[jj] * DecayProductsMomentumXCM[jj]
- + DecayProductsMomentumYCM[jj] * DecayProductsMomentumYCM[jj]
- + DecayProductsMomentumZCM[jj] * DecayProductsMomentumZCM[jj]) ;
-
- E = sqrt(p * p + mn * mn) ;
-
- NewMomentumX =
- - betaX * gamma * E
- + DecayProductsMomentumXCM[jj] * p + (gamma - 1) * (betaX * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + DecayProductsMomentumYCM[jj] * p + (gamma - 1) * (betaY * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DecayProductsMomentumYCM[jj] * p ;
-
- NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + DecayProductsMomentumZCM[jj] * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
- //Set the gun to shoot
- particleGun->SetParticleDefinition(neutron) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
- particleGun->GeneratePrimaryVertex(anEvent) ;
- }
-
- for (; jj < NumberOfProtons ; jj++) {
- p = sqrt(DecayProductsMomentumXCM[jj] * DecayProductsMomentumXCM[jj]
- + DecayProductsMomentumYCM[jj] * DecayProductsMomentumYCM[jj]
- + DecayProductsMomentumZCM[jj] * DecayProductsMomentumZCM[jj]) ;
-
- E = sqrt(p * p + mp * mp) ;
-
- NewMomentumX =
- - betaX * gamma * E
- + DecayProductsMomentumXCM[jj] * p + (gamma - 1) * (betaX * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + DecayProductsMomentumYCM[jj] * p + (gamma - 1) * (betaY * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DecayProductsMomentumYCM[jj] * p ;
-
- NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + DecayProductsMomentumZCM[jj] * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
- //Set the gun to shoot
- particleGun->SetParticleDefinition(neutron) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
- particleGun->GeneratePrimaryVertex(anEvent) ;
- }
- }
- }
+ //Obtain Mass of daughter Nuclei
+ G4ParticleDefinition* parent
+ = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation()) ;
+
+ G4ParticleDefinition* daughter
+ = G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
+
+ G4ParticleDefinition* neutron
+ = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+
+ G4ParticleDefinition* proton
+ = G4ParticleTable::GetParticleTable()->FindParticle("proton");
+
+ G4double M = parent -> GetPDGMass() ;
+ G4double md = daughter -> GetPDGMass() ;
+ G4double mn = neutron -> GetPDGMass() ;
+ G4double mp = proton -> GetPDGMass() ;
+
+ G4double Q = M - md - mn * NumberOfNeutrons - mp * NumberOfProtons ;
+
+ vector<G4double> DecayProductsMomentumCM ;
+ vector<G4double> DecayProductsMomentumXCM ;
+ vector<G4double> DecayProductsMomentumYCM ;
+ vector<G4double> DecayProductsMomentumZCM ;
+ vector<G4double> DecayProductsThetaCM ;
+ vector<G4double> DecayProductsPhiCM ;
+
+ // Initial Lab Momentum
+ G4double InitialE = sqrt(EnergyHeavy * EnergyHeavy + M * M) ;
+ G4double InitialMomentumX = EnergyHeavy * sin(ThetaHeavy) * cos(PhiHeavy) ;
+ G4double InitialMomentumY = EnergyHeavy * sin(ThetaHeavy) * sin(PhiHeavy) ;
+ G4double InitialMomentumZ = EnergyHeavy * cos(ThetaHeavy) ;
+
+ TLorentzVector Initial = TLorentzVector(InitialMomentumX/GeV, InitialMomentumY/GeV, InitialMomentumZ/GeV,InitialE/GeV);
+
+ // Array of masses express in GeV/c2
+ double* masses = new double[NumberOfDecayProducts+1];
+
+ // Filling Array
+ masses[0] = md/GeV ;
+
+ int ll = 1 ;
+ for(int i = 0 ; i < NumberOfNeutrons ; i++)
+ {masses[ll] = mn/GeV ; ll++;}
+
+ for(int i = 0 ; i < NumberOfProtons ; i++)
+ {masses[ll] = mp/GeV ; ll++;}
+
+ // Instentiate a Phase Space Generator, with flat distrution
+
+
+ TGenPhaseSpace TPhaseSpace ;
+
+ if( !TPhaseSpace.SetDecay(Initial, NumberOfDecayProducts+1, masses,"Fermi") ) cout << "Warning: Phase Space Decay forbiden by kinematic";
+ TPhaseSpace.Generate() ;
+
+ TLorentzVector* daugterLV = TPhaseSpace.GetDecay(0);
+ G4ThreeVector Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ double Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(daughter) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+
+ if (m_ShootDecayProduct)
+ {
+ G4int jj = 1 ;
+ for ( int u = 0; u < NumberOfNeutrons ; u++)
+ {
+ TLorentzVector* neutronLV = TPhaseSpace.GetDecay(jj);
+
+ Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(neutron) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+ jj++;
+ }
+
+ for ( int u = 0; u < NumberOfProtons ; u++)
+ {
+ TLorentzVector* protonLV = TPhaseSpace.GetDecay(jj);
+
+ Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(proton) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+ jj++;
+ }
+
+ delete masses ;
+ }
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::SetEverything(string name1 ,
- string name2 ,
- string name3 ,
- string name4 ,
- double BeamEnergy ,
- double Excitation ,
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path)
+vector<double> EventGeneratorTransfertToResonance::PhaseSpaceUniformGenerator( int N , double R)
+ {
+ vector<double> V ;
+ V.reserve(N) ;
+ double Norme ;
+ double Buffer ;
+
+ for(int i = 0 ; i< N ; i++)
+ {
+ V.push_back( Buffer = CLHEP::RandFlat::shoot() );
+ Norme += Buffer*Buffer ;
+ }
+
+ Norme = sqrt(Norme) ;
+
+ for(int i = 0 ; i< N ; i++)
+ V[i] = V[i] / Norme ;
+
+ return V;
+ }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorTransfertToResonance::SetEverything( string name1 ,
+ string name2 ,
+ string name3 ,
+ string name4 ,
+ double BeamEnergy ,
+ double Excitation ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path )
{
//------------- Constructor with nuclei names and beam energy ------------
@@ -750,12 +718,12 @@ void EventGeneratorTransfertToResonance::SetEverything(string name1
BeamEnergy,
Excitation,
Path) ;
-
+ m_BeamEnergy = BeamEnergy;
m_BeamEnergySpread = BeamEnergySpread ;
- m_BeamFWHMX = BeamFWHMX ;
- m_BeamFWHMY = BeamFWHMY ;
- m_BeamEmmitanceTheta = BeamEmmitanceTheta ;
- m_BeamEmmitancePhi = BeamEmmitancePhi ;
+ m_SigmaX = SigmaX ;
+ m_SigmaY = SigmaY ;
+ m_SigmaThetaX = SigmaThetaX ;
+ m_SigmaPhiY = SigmaPhiY ;
m_ResonanceDecayZ = ResonanceDecayZ ;
m_ResonanceDecayA = ResonanceDecayA ;
m_ShootLight = ShootLight ;
@@ -764,3 +732,5 @@ void EventGeneratorTransfertToResonance::SetEverything(string name1
}
+
+
diff --git a/NPSimulation/vis.mac b/NPSimulation/vis.mac
index 7130261ec..e8345fb35 100644
--- a/NPSimulation/vis.mac
+++ b/NPSimulation/vis.mac
@@ -9,7 +9,7 @@
# choose a graphic system
#/vis/open OGLIX
#/vis/open OGLSX
-#/vis/open VRML2FILE
+/vis/open VRML2FILE
/vis/scene/create
/vis/drawVolume
/vis/viewer/set/viewpointThetaPhi 0 0 deg
--
GitLab