From 6474bd2f891e563a1e065639b89d527d2a56aa03 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 5 Oct 2023 17:00:03 +0200 Subject: [PATCH 1/7] fix: (Layer-based) Geometry construction handles layerless volumes (#2514) In calorimeter developments with ODD we noticed that if DD4hep produces a volume without layers, which currently is the case with the calorimeter implementation that we have there, the geometry construction will fail to produce a geometry. This changes the layer-based geometry construction to ignore volumes without layers for now. --- Core/src/Geometry/CylinderVolumeBuilder.cpp | 9 +++++++++ Core/src/Geometry/TrackingGeometryBuilder.cpp | 9 ++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/Core/src/Geometry/CylinderVolumeBuilder.cpp b/Core/src/Geometry/CylinderVolumeBuilder.cpp index ceb72c7affb..a4520d3b828 100644 --- a/Core/src/Geometry/CylinderVolumeBuilder.cpp +++ b/Core/src/Geometry/CylinderVolumeBuilder.cpp @@ -161,6 +161,15 @@ Acts::CylinderVolumeBuilder::trackingVolume( wConfig.cVolumeConfig = analyzeContent(gctx, centralLayers, centralVolumes); wConfig.pVolumeConfig = analyzeContent(gctx, positiveLayers, {}); // TODO + bool hasLayers = wConfig.nVolumeConfig.present || + wConfig.cVolumeConfig.present || + wConfig.pVolumeConfig.present; + + if (!hasLayers) { + ACTS_INFO("No layers present, returning nullptr"); + return nullptr; + } + std::string layerConfiguration = "|"; if (wConfig.nVolumeConfig) { // negative layers are present diff --git a/Core/src/Geometry/TrackingGeometryBuilder.cpp b/Core/src/Geometry/TrackingGeometryBuilder.cpp index 911909ad1b6..90a10b0ba04 100644 --- a/Core/src/Geometry/TrackingGeometryBuilder.cpp +++ b/Core/src/Geometry/TrackingGeometryBuilder.cpp @@ -47,7 +47,14 @@ Acts::TrackingGeometryBuilder::trackingGeometry( for (auto& volumeBuilder : m_cfg.trackingVolumeBuilders) { // assign a new highest volume (and potentially wrap around the given // highest volume so far) - highestVolume = volumeBuilder(gctx, highestVolume, nullptr); + auto volume = volumeBuilder(gctx, highestVolume, nullptr); + if (!volume) { + ACTS_INFO( + "Received nullptr volume from builder, keeping previous highest " + "volume"); + } else { + highestVolume = std::move(volume); + } } // create the TrackingGeometry & decorate it with the material From 3201e75b2a11af565c7fab5e07c2892146f160b3 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 5 Oct 2023 20:38:52 +0200 Subject: [PATCH 2/7] ci: Include ttbar track summary in physmon (#2490) --- CI/physmon/phys_perf_mon.sh | 17 +++++++++++++++++ .../reference/tracksummary_ckf_ttbar_hist.root | Bin 0 -> 40994 bytes 2 files changed, 17 insertions(+) create mode 100644 CI/physmon/reference/tracksummary_ckf_ttbar_hist.root diff --git a/CI/physmon/phys_perf_mon.sh b/CI/physmon/phys_perf_mon.sh index 8b4e023be7c..fc6cc6f3d3f 100755 --- a/CI/physmon/phys_perf_mon.sh +++ b/CI/physmon/phys_perf_mon.sh @@ -1,6 +1,7 @@ #!/bin/bash set -e +set -x mode=${1:all} @@ -296,6 +297,22 @@ if [[ "$mode" == "all" || "$mode" == "fullchains" ]]; then --config CI/physmon/vertexing_config.yml ec=$(($ec | $?)) + Examples/Scripts/generic_plotter.py \ + $outdir/tracksummary_ckf_ttbar.root \ + tracksummary \ + $outdir/tracksummary_ckf_ttbar_hist.root \ + --config CI/physmon/tracksummary_ckf_config.yml + ec=$(($ec | $?)) + + # remove ntuple file because it's large + rm $outdir/tracksummary_ckf_ttbar.root + + run_histcmp \ + $outdir/tracksummary_ckf_ttbar_hist.root \ + $refdir/tracksummary_ckf_ttbar_hist.root \ + "Track Summary CKF ttbar" \ + tracksummary_ckf_ttbar + # remove ntuple file because it's large rm $outdir/performance_amvf_ttbar.root diff --git a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root new file mode 100644 index 0000000000000000000000000000000000000000..d9592b163a99e5c0c4df551d2685e97bf58b0ae4 GIT binary patch literal 40994 zcmdpe2{={V_wc#blzAQ#H_1Gw%=3^jL?|*}Q|4=)$B;xQii{OXNJ=GR8Iwpw$(TZ9 z%shSP97FWJ{eHdg|NFk@d%o?tgLC)UYwfkxUen&|xH&t!1Hf!G00341kjVl788^Z2 zf#3@lx`qwDp#cDI3;?)n0K`7e^VSy0%{&gfWC478`s&AA}zQ!q*C0Peo9X*hv#@j`46qJEEQ3F~Av( zeTNYM8y4V=z?hloTbOy8nI)L@yO^29Tbp6aqs*+$`uktaYHITE@v6Pje8tUa4X9xR z1O$~rl7J5mL>jVa{tzBM6d{LoUp0UvAR!_;MPw{{62VG9fn=z!K6Bad^688AtlF^z zScF(upJj^Gc*p8GpExX83OiOZeZ#wkYk`F>4gUBgn2v~X+>B?u6bnm%uglY%s+(Kx zQP?xtf(Wz7u!7boyBQW6O~fxd1A_%z`M2 zFt>)$sx~A}4Ypg8VdAX4}9WviBihUKB zqv5@UL{}p3n2uu)XkEw3+Rut|CrGhynF@QurhNt>uN#965K5f73Q|%FrY#dpTP}#U zNV6T<;yY=1xI5UP-C(*hA^@f>x?Y9pN&(ZA%w~tKEwlW*O-&lSTp&$dVVd3yaX5y{ zmYI+&OG6uvL-oMf3_#%E;hq=5EA+E zEQ%-d23@@Us#9SNl_=|wd67o@E+bi3NZ5gk5jQ}lHsW;iT(rO(ZnTb{jebZrZW~WM z&n(O8)F0ea%XR(qKBAYjB%~s|fEa#Q;ETryo^_*J*ow?4bXeK_!YN@Iq5=Z@go@>J z_D$HZQdjS@)()ONFhRu#+&`dSiA5-;Q4-wdNQ-ixMeaL$56DulE;M&(6a@y@gE?^{ z_u=j9%R$IXgmU8Ik#oF4WZ@Cw1820gKZE2v43m=%CMP>YPOPFGa^gCvI6I(WTFwCg z-4-or%pn;aw18xk%zT@cW-n(OU#s$30~a=hb0`=loE1$R^deh|qCiXpv)Mq@%;8H_V#GS~!wGU^XtqQN2h)+CdxKMG`Ag zEl*-?cA-WNTw$^!1IkE1`ISN%NF;NZNPIAngdieezuh4cx#Q2}tncjL;fRLTG#dio z+oIQdSoThCuIU}wBWwA&r14+A=6t1z2Dk>y%q%V-Ik8GU(4_>cnIcac%hG7KRE8k1 zsX(;9(LD~nNuWZ-iY*);@{r{|%LBR&!(wq7T-^E~vv+6S#}r{c%rw}3^p|=|Dn$MQ zx7ijE8;@ln;8GlRzqkWb3{T8x_vJj!q!>cRi|ElsJs=zR{#6J&XfK_lkV3JGsWIZq zW>HFicLJh|+U{{O;wX|OGT$T84_iOT{Lmj@!Q(8(R* z5<2PHn!BN`RO~=i2{K{>>WiZW_vL{49&y;|i)U$TC#37)2wngXXn0VDpB-_}|Lo|b z>u&D8Dc8aXfMHv%En&Gfhvgc*BiFA9`BZrz!B%sDH3aLXh6w#tLr_2(0teC%vF4i^ z;ygWW(y46?k+WArU^X>G&Mz8*`=W*O;YO4A&ya?&`jt6NTN`f`~Z-JxB=O9zzRoC12nZny6$dhb4RqBs*|-dg`%^Eg@X^u!OrO%+6w%@ zxXS`uhX5StE>ML-Z(t(JV1e(1dJ4`C4rohvJ7*`r z;PXpxp#oT7*Bp2NEem@&H#c)1=;q4+)?PSR5HN5&k3F#D%`s@`swp%Ka0A*t9f@{z?PTORRS*{WC6fu4`6Vz zc1|`BXmbDwz8k@yId_A`u|5IPlEw*qLpfWc+-=b)%b&9c0673u(BSL04G$tI2t^hE zLiYfN=c4E10#DrrK!P1WbqBs5*$ZFS&KHex@^G|3yFnPB9X1^pl%3OFOd#4q-p+3u zmYtJ3G}pO3Xh2c|UG0Mz&A0;_BtD6rF4}b$uU8!;>Y>rW-uK}#cJUg>_?6c%_rR{_ z_ki2xwF`g+{$>qNSY!t%IQVYqV)_XjWp06i{7ZL~1t={jO9yidq?rEFw&Zc=`Zn{i zAgsstpt3#l_wdXqcaQ)P-&{$6UET9?a3mb?*8g)8{)Ni>L7ID@oh;Gj7NAT0Q`S&& z@(Me%+Vy3dDchH}Ck(Tx_ircM2v3-H_k{8FRLn89koo{x0#$3E!DPW00=|7Wmt%j& zW%oUhwz+Hqb6F4O^5GqL;NY7C{u}KCLH@OEv0>-(_7rhI9B~%*0N&y$&Hy|`njMfJ zPxtA`xx34G+hG7>=M8CaFmSEmJo$BoIJv*hkPw0^3jm3G0PR`5-T)Hp;08}qaR(*< z3bXYID?3j+45%qF+$VO{@Xj?3LD2gW;Bq>G;zogAxtxT+FHVA-;H`hOoI$8Tzemdo zK!P1iV6>oU2qePhxCGV^9l*uq>_!8>Z>u)YOrTNb5REo%0}{I9cMS()>+A-alV6nT z9?F;iNU(z;i~{60?4Uqstzm9~c7UOi?%?16UHfOiy?;OdvoLT27&zo&?*J#*1g`At z+Xdgc%e3Lv|Eld4Mof`)99qpVz(EGMt;K$!LTq`JQLdyccA1?7x;srx+x{&AdUji05ep3Wt z{GoyQ1-vpuYN{QAquiZQR(2Q{2QVi5bGH8_Li{OfX8bn;kWZxXoikrk1bavu*wzl4;Js5(5>IRHrJ>vvO@q5tw8A$JTg~1HK=3tO7 z3t1jVw6e3Cqq#e~BA5nZgRucUFYp2n3hgb95)k7T6yg^XLkS9r9}yN85CCtJ17H6| zTzBH@y#{e}wP-toQ7_sJS~0tQ(L=jY!`pPpFeFg?v&R<$?}D`j`0@9RU?*&>zOoZe!r+Y-ML{2W4!wJz{Y9|BMPQoc-YfXNq{g2aheg zh^U6`BHhl+K|T{eF_9+ddw~?<4&i*ZHkN+mz1-IoKfQ=7a3DE4oyP5yOSmIw93)5Z}0RCNE z1p$QFN(y}gEf3oKUvT!{!-n~~3-&)`Ekxnwzol^LZz!w*M$AAGo|5Vog)jbv!oI&! zczfw!_5j)wHXB0R0#5$bodUx*5Fkup0_fYnAnv~h3lnz}>_4O{IGKMB*1*on-F6eK z7xWeQaH>!Jc3P_nnG;Qvu}R^VnhWN9GP3k5v?F@sFB66|35cfL8~D^Wu} z5(w4-ba?i7N~FL)!2zOs;~pfo>0ScUooa^&AcMC8Ghi#Ug2?_4R_k9>{{O8li@nj_ zPfu%iif50Z+1#=q-(vJ04*wHx%4H8-f5P3Zn@~IS1=)tXnOh**0=%71`D2(MGf8;% z0Nr9H2@T9l@*Q9xGqJ!U9RQ<<>x9xZyDf>5hTUewm0W9#h&8T`$lmq7xAhS;O)3*O-;9H&t3>bRb z4g4=Q9tLg!1K-RS?+L1Rd8*`F&eMN}z2(mSE$p~GQnp8hx&ugA0C0l!7N{V=Kd?a- z8_)WLCECFO!~JtB4hpk51Qj`XQ1CA(4k1NO6!$59El;!?n4|UjCuxqo)&rDdFw%#F z)d)a>@0)v~dv;L3Fds^rksHGKdXyC!W9eq+;tu=JQUCySfv?{-96cx%B>SsX{v)ub zVX&wFZ-YGpgFW*<0(fhS@YgpmK$G?Dt5bD}a;*fOD`UO6_!EBHFvP{HHS{@pseAG|BMh(dn97hOl;+$9RZ{)0NB85e zuuHZp8w&=l@DHHA^xWL;{U=H3=-=bw1mj{0<8o{_F1Xg9X7CbX*~A1~m!QsomgnED z)PI7E=N@E!I?210<(}LZp0$<>#0GNcRJ5~=n>mOFd;!vxf9xuQ6i?d&d6Ol(lAQrC zDZF4(DDNf(fwi)S6MS|H!^bPVbqB;8=xi876YK4tQ5@D_q=qsFGuj>)v>UWdB?G=7 zMe^IWm4uZ2n@y^p>hv^#1Uu-$xb5+22(9HktgX>*P>PIb>pN;IM}`ahm%tD@-UJSN zZ4X^w3$}|aegG2evw6;HPlSxO1&+aJ7Y=7pkR9_uS-@BU003pA|Kt4eHg}G=pnkh& z1BELf;Ol=yZhHqU2ta}zd|;yNaX?7G%?0@2q>`H(C_4u?#bwj`0QYpD@EPQsoh4}g z&~9$d;1TRU*kvIH@evnS0!Ic%b%w?TffxNAyCeV!cDMv%r@dR8NUWhVojV?{Ftq38 zyqRLKb}$FG9GoHH0n-guXg7Z6ZO0hejM_o`j}WHZM0gJ{&_6~P?7?|E`j>+PMv5Io zK`$Pz0RWcn>t_HUeZV^pz;}WRXFFm<5dIo5LVg_Ng#X1cfnyabXc%Dgh{v8iZiv|o zg1m}$j;Noj+~PkTTic3cA$V}xj+DUA_MdL|;_1N?QY70^cAGK@k@+vpRcumt_r-0v zAi)3x1C$NI`hx?);JNEB0Kl?+<{v)KO97wfjknr9&)eTh*s8*#ra9fn165ME!Do7@ zHY+KVeyyaShAJsY;7W?Rtx5_9ly<~+C51Te-bxA_xRT=fFO?Kd=2Yd}a-Lz5P$h-U zuay+G+m#e_Ka}zo<@P}(6e{#s7jZH1_;|Rua=7}otTk&FZzAdHhzJCvTg2&tZngmu zgvy|K__O2CAXQosfEo^1Ve&BXHFv2FRYM&8b79W#1fKK=D%Gw^VH*NMZY6md;hC|3 zcRwl*J-6qO4ElttAG{4si_KijlT+`E6#QHF;nAa6v z69PECO|OTzb5NA;v@c6q1ZJDrxm!Q7G&4_UC1=fWe#=_Ne=vheGMbh`%VL0)KGh;y zG{UOjW$Jlkbo7gHW|3=F`}0!rM&7*_-ViM(^tcjt)T7e)`!%17{MW2mt|xUe#EDLx z6Ka3<-DoarGGX)_jj63dIi`bU)uDc4e^8IhQa3}hkIc+-rFx#xBZqjWl_=wQOgCZ) zAFNm^S9V-j4T)i;HK-3`60!PB_)6}*vDWahqtEY|iHHKqiSGx*w3b$kpQPv( zdN)5zaXXarsls5QkE_E?{e|cao%EDf(bCGRkBNM-3f3AAY0VAf3#}DME=GD2{@_P* zA$*%}n17j%b>qJEDNU!)yke~X+mrglX?+YUJoWVQUGfE4FK6xQ z%yNTxoaX_30PImmh~D%F&v8{@N|eA4tC3CK@o2lK1Bx0Z7e*+a;z7ph)4$Ts*xl*+TwIv zrfvw9^y@$tH#N%rU{W8^S0Pa~^PpR8)+5NZ*f zyw)j&>hS(nz`+~`=&FD_muKKb0L5W|80`*`Nq5!eD1g-=;15Vwxc^DGKU{7i5x>)) zaEmEheG&h#?gt#_ckQ2{39egNBya{+*>8+gI$WvnrXJjD>QeB;i3z2aXOkux%_mIP zA)TY}i^}W@liyFT$)fqlXC?XS6{SKpHf0Vqk*7ShLgOWhx+Qt2;z*`>iqS-JDsuOw zwZOjFIC|WQANyB=TE#XV^evgz5A~|@u?-W<2(+uY<;9)hW6=r}Hr5FkF%xAdyz}8% zSxJ&!SY%pj&r{VppCP^24z8IyD({+k-f-6{&JCSvcqi66v9qIw%Uz<{|PQN;`_%vK|{Y{*|1>>d8Z$i^u42$lBj6SR+ zr^WjRsQf7+uhBl$Nc%v-92iO3Ta(4{Y=enpJ$ymGrhLssb{WsVak8|!JFd>em0Tl; z3~TB=iFXdRNd{?AZ1r9Gc2fSJ&}6l`xndc_ivM6>S1P-HO?ufNg}3X0XJ)90t542` z<$1-NJvmwW(53SEs(O?J>!dJxFsVf>S*h>C@F0}~ZQ@1AYYEE+-cOHzBIx!Xr1l;Z zl^ngiD(N7FYW66*dPs(?6Qj4BYMCc;#DOz165V&`^w4|m{9JS0at;hbXybmKeNkhY zYM+vnD%wfMh88>$*d1=k&?uVvhoqc;NOZh`=fEoRBy|UvF#dCS2Zc@&37p&%y0a>Y z=9@xi*p@5jd5|z%;~CwXQ_SAVXWzz-XQx`w;vT(jtHeCxWL=q$MrYlPjwZkGhBiP0 zlse1kTR-G*Le9G1PgNLke5jS3nUXqQ zz2fZ**%Hl!jN4=JN5lPJEA9)OOqF&k_xgrQJ$%{dsAhg+5z8s=GSwmu%F{yzowaE> zjpx3XSVX#!>K#xj4Q@PYxvX~VhvuRP@1=sC#ekcd8MiDS>dGsYzw74VbIV>rXEPk} zTT1)jDrFQD5_G^eJ(x*n{gmlOrT2PB^u{r5soAEsJ6CMxcq0hz)M9vg2QJ*9{`h`* zR{iGUaW?w)?jhTRPlS}bhl_mm?0we#hto*rD0Sa9`PeJTHb}a5AFB*P^j#@w>chTV z#Tfbd;ZqICN4U~`$m7g0&mWO^v*Kjbn}2=CM`^>=pJO^Ta9q?Fi`hR^;lb3py>}5a ziQ0x4u+(EMI;A*Hom}(7Mi6__A+k(u;FZ!;J;$?j^029IBHp(&j7e80i}_Ld&wC|Z zna7tpzU;TGsNNLB7IhJi#|NCw7fwmjB>H{IiavXkDRUd76M>)9A$Fvl>D( zR%+imri|krQ$NCeac=4dr*}R3l2VtOY|@HTzGC3XiQ(`cBB%321RV;)8@_z1=W$!P z_qL-r>ui$PocYk)G3*z|L4k``S~Dbj-!weNswJvcMXhY=9}NjyM(hXiMdA*x*Qub@ z15aOq?m7>wUrY`N{TlQsUV zm7G>B(BQQI?)1a{Naf9uiEQK8Gv0#Ovu(*Y9 z_XjGUL|nAN#bK2u=osiWMOYysGs@U9r)1P~nuDl3WiMCpw+Gi&1`hfY=D#5s;qau)$J<=?SGN02!SfdSbhqx_Co0tV95cn;w1z|8sM)hAPL0c2B)E z>4;*Z6zQYLN0q`ADFx{j^Ay9e4=SqRp()f=i5?!ktQ`}+h{f=r&Pw1z7mvQ{Ys|vx z$#)jh!amXq73ISd3(u0sSh4r5lhtAiS2xHoAyCE-#LBNC{hm*)YY2k3`1<*p< z?ekSmO{dQWsqfHH&#zggsXO~6$kBGep}C>2#W+Cw%Xx{7_wl4F5$Vn|UlT{fRwC+t zJm-v4@cs7fK!yEBUM8(nVbL=>;^f+KC}!%o;IiQaqDrI6rU42qoQ@mh5B2ty#C+pg zxsjZf<-Syyp8IkpF!AjV4DLg~2^(W-uK0dv#Qe#*H3*I|&?M_^uS!@P!Kc|Z=@L)%0pX2T$k zZ-D2VA&!Gj4VZ%~vMJz-?0B>7RnxD=%XLAEI9lr^k2%o1S>8UhS-g6Emwz834EpzO zkbfUzMX}}IC-}(Xa1rqk5nbGMzAks0A4fZ=Y{F|1`hUNk#(D2Akwf+Vj>P{Cks9w93~cy z0OIZ*&-3&K$%Ewjb%m&d8fc4Z;0p1C4*YI6@VEt}o>*ZLL#H1B`AqP)O%fAr`F}gA zi2wnxZmA}HST#{>swRsaA~$XNptl)Do%%JOa0SXIFhTi*>)ZJR3l{C@ z?NCa7Z$5#3^T`ZPf5|7fnNz*wmb)6}59JfQe$6MiRJv@sm_3~1qrYA({sph8C5*_*M`U}+Or{T zV`-=#`qcX>6eydELfg&_%Z)??p?-ve<=4sxd3O1wUNF226HL+%ZhOmuFsU@E0%^?y z3$_LN{3wv_=ar-QJZlI@t#evW7@DDS0(2AVhP z*sV5I!G&K{!AVFJ#9EPTse(_!nWfwE$?}VQVng{ZIy$Uqc2|0y;3QRA5@c z53*2(bi?Ml#&tJ^PI*IW^B@vX*>rKM>raE^fz*q$wKe1MXzKZE`VTVR>eUkCR`N zIUj0xmnK0isLS3yKvnJJeoU+$6@fB4UVBu`37U=2+{BHL6U(36$5=$8htnrjRkK24 z%3NfI{U$#5F1Wo+5>vIittWe$?W1bI4ErI?r0+#kTt{xHCIntRng23S^IM==(=us3 zGXG{WW~hOR@(O3XkJoFyr2N=U>*hzvJS1&2*_D0Ejs5Kg8KIU&4fDint5c26pU<79 zXU&dYjF99F(E$5$Xd-TIa=i%eo7c)3=3Lr`#^%;;73oT6bR`C@Xy^|%SV|Y(&GNFy{+@I>LrvLtE!8hQzWYkKzmKAk|Iu(- zi;*H_+;Sgp?Sp*mjSpFPUKE@Ah$%8f4k1P zS=?G@9*kLoMNZ7E3X^t==g{gKH*Ml zDm!_YJ-p|+(YZKa^5dMKgS%NrO}V$6biw7tGZ)uGefpMdUfCV_W?OTMOqpPsS2^9o zpUiY^Rjo>MeQm2$lX?TojfC*@H0^{(Z^3Jou1OsP_lh~vLK<= zODOC*P&u^Px^Tm9^2uanLz0!6KwUX!onio0;&J6rj(2lb*@i5xx>Rv}BGgUgoP6&x z5^f!)`}4v4pB!{L?6lpz-jm{mxan$6s6V$F8nFb2=vY zMDKi*t?x$9$(8pj>uSoJz0qOOZjVMksC1Or*`HXiBK|IsN^%X0Qi4?a$-tZZqe(B< zYrgjoV_2OMX}_*VRt`0EeKlYjxn-Kc#N&IffKykNyEY`+my^@r@S*9+8wy8;x|@>= zo*o?I`0}!G>G})_Gxg!mNy`hLqWV;)UY(#`l>C02&1L%5yJrIUC&c>G@CvzX7)MSs z^OW~T;Y_EI-pIfn^Xq?xMI-6*KJ>8=?pY~p=B0i@+;c`|b6ke>`uBz%hdg>ox=vP& zDs+1Ea4CrRRU54JFC!io_&*@-@g^=D4Lp2yVO(;xU)tWU?Ul2c50S&7*%;wtA)@zE zF1>gkaAKYYwucou1i#J*<@T9+isP!eaL+V+qYcm5h%7ElZm5>z+h;yQy&KPoTE$d^ zXK?1;eJwagTsd|cH>3#rbW+U*qemRKXEqb5{Ordp&2^+gR`!`ca!5c;w7TuGC&<`B zqEp4D8GU-yAkPw)`g=b|sG9V{?69!JAhl90-UvNK1$7nOtRlxV4`Xey1EQyP3=B~Ju$h``b2EFcY*2hp%O)Fzz4a~HXO)TB3l{cTW}03 z`PHFB`Z?-cym;4l%qOVR+@72yVy48t!|~W7FBi{k&SedA zuW^ENl8Uc%e$IT<`Xnhi+@~4yaU5xq;Ol3*yKf1{6@y@b)pX&s>f7^1e44AnA}+ z>fh^sh2qhSecDX4Im@HV4M{5BG?tmlR?~kldmBkCDZU;MQNu|t10`kC!tYWtntXed zj4x!f*8eIg_@^%hWs)4YV3k)_BX?HM7rWt2)&&}Tx`R{_Dl2;S*VoHq9ql-Ud1(s5 z$VLV^PFx~&b#d>=ziHI-bl-_VP{;7!e3Ts{;u5$ZG5iVB8XyvS5yh+Gage*3(|ttC zfT7V-xj$F+Tfs|S6)k<<6%WSTDV=F5oCKph4^^vL{m>m% zE1URqHOu)fzfUxim|u!&?_{6KjOvS-2ld?bax;xj$Z3($GFdIt>ACi8=|KIYNa-+oCkz88P=;FwHf!l84dP4_6S zwQTrD_@YydZ456llfu$nv~A2;3z!>JM97i!3>fAxOJoOm5uY!7E^%R*ofeRP7)6Gv1ajeZvgTFSo`T$w)tQiTLzM*;nyL zu{36>r)TnYVR4%;I|Ce%JQuh6U}1D*tzCJ|tsz2nwC-C>bgS-5rYi?teU>V+>RNeN=0i@> zFwSlq)oO2hNW;~T(O@cvB286tRNf}RCVUB4&jb=FB*|{@WfF&JDF|HkP@Z9@5Rul9 zyceH4$I||pVZAY|>qpiDv6KrLQOwS@%rSu{Wl7Zpri(q0?FmHzbDw4iT&37w#Z^{k zL{N0jy?IUaf?ux3R5+4}U-n)`-o$fePjnAOTGC)-=GRn)j}D;<<|a`^M41l=c+$h3 z#_wNgA?Ro`PX{L+PQlGwH`S$)!Cv5YVV%q7bLXmD$E&kkd`^}|U^$wz zkvp%hKBSb8yO1R4Uoo6j@H#>@V(}>HGI17mt>}*M?8+?YB<~0h0a{N*Nl$KDeDIN* znT3aMWNayQ81bf@DGv#%-{s@8$EB5^HsD?JAQN|wYLmMaRJep$QrJ&nblvuyF?TCh zo@hGF!r+j%sF*&Pv9614TUJaMYS3?TMdFFQ60K?D%CpJIit=PXS~+V5JGpPzCCIee z_@$JqxFt`TTi$-T)`U}4NpK`tiOZv>HFa9&X3O;u`Y$bwZ`qJtUkTgfT7TqHddc4| z@48pSI$=YT_MF^oq}=G3uj`nClfU=;>Med}rY@ro_pj08Xdg4}9O}B(X&~Jrvwzr> z+U$d4z)bafk?-&47dFPuPb__Z#Xj%2G0Ya}>^)ZHXnV@Z0j5sm}e0 zbjBk=-3Q{BUVTJIc^v$J@4`rr(Yt*sNE%}i;q!coa;Ys#GHCMK3`@3_^cyp^3N;Y| zydM9%K8ig#WnQ+&-#YpZXMGj?XrF(^nKq@nyObthe*g9E>o0xYuNl%FB#xVVsx5`c zkGzg`CrEvIxka7n?7k5TheqR<9VZoI0=|72J&JdeK_EGT=G~&WOsM0a7-||rb2|RQ znYT+%8i}p=GEZgAlJiEso$N@+RW4-_mZmkG^{@6ldQj+Fv7U=+O72_LNiE=5>I!E1ll5FdVpq2~wo+kxTY0~yo1+s0H^ zN?k#RNK=OM>?-yeLa8(}DvIM2176*$x=`C10wq_DQoUjiF7XRga&L9Wuq9MoUHpH@ zMD@Ri94?Cq%^81v?z0*P*@sz0o-Fhjg_FqxHO%MtA%BF$1zI73Ra}C&O>s{i>R)h9 zA*wT<(w)DO&7pMjFt?cl9)f`Eu*}dXBI}K(U6HMX4i~N*&=<)Y#qN}zLnt+Re`!Ja zov@%FdJu^Buo!C6j{l$)2Z%|y_jk9I(A#_VewX#L?vP%x_P2T7hsU(uExX*s>+rtl z)zxC33Ga5!aXJP6C*JmY9iLlIf4nKzPEn)~EWNKimxb`L&AeV)&l^n=-=v^@_ZKE_ zcU`U5q0#G{6a93XJ1&YZ?L%_I8n%XaN^)35MN3%Y=0Y|3nY!osL@d2Ta(7lR*z ztYdg2pQjIfR*ji|eV*MPJ{2w1%^W!SaSk)I*Q22Ba5r-)yKXPY$2vpZZ*{ zao*=n?PGgHn^c_3LfXm5>{(aYzUbb#buX$@4frOxCY_vd`3-@}Z3 zx$#)C;ea$Exx)iI!-hQmD`WK`!TZTr*h`|{4#nQcoTHMX!<`cGdg5Yixr#Mv6+?BM z>sIE$BASrkl&A};k6%a}!XJ<<>2-7)lXE#NCyeXnb)$>u=6#~T&>Ckm>X+0!R6HkB z^!;tpy}l^o?6)xwoNZ}o9Ute@4dHy9v_G4L@xBCh`*qN zzQ#Gn`mvrU>j?bb+QiLSRFobhBSp@{_`&G|vF>fsHEPNmw?m(CaV}-pOr##{)33Ce zOq?9d$8J62tvT!Ap?tPvh`@lz_y-xk&`bcqV5{6rg^?d=44HBsA&IVL6S~n*Bi!&>%6*R@Gh#@$>C-3SRZ6bNY<%Ucj1pKzOIHZ(>BUqKcZu zX4p=^&9hKt(k~^u%G``Q_>}^k-K9&HMYvo#s?QfDRy15!3gV5(ihOJx--$(c=hVZA z@@NCfzJ8zcFPR3r@ns;MtER~$QJEx4ylNJ$gSg234*p%a#ZI{CnTjvRB|AlmmauD2iP#rB1GGqy{*9+N9Wb<~8eyDSNPqD_PFl9KVtm$ePhKobXKPm`<+Rm#gX@ z7X?m@hh`kJ?`rf-y#7h$8b|K=8?TA3WSyyc_j)4I;+S7L33=?={5dv0$3RD8G3kJT znIclHvj=FqTd>LN<|px(rB){DTN9@F@>*{QWh1-NVUx&+&#aKD^!y^~3a)p{fapAswJ_zwsCo@1F8s+|Py7^#&ahQsfl z?|qrvCdxY~9D0db)F{}Im<^Zcb>vfCXG(U0rQ+AAgJlu-w-t@VDDoUg4mdRx1{g!2+xr=MB3UIB_7AqI7PhX5?VMy^Qe}FI6 z!QPQWXGA|qbC%`fa*y?Ruhg2*;~+=>;ojG;i&D>Cf9UENJ(6EpGCSJ7U?6hZ)`FIV z?Rsghr#gGC?FmPt$H49D}|*AU&l^kx_B#?6;ZuB4%dh z%Hi4BrPjr75^L6G=78ZcmPN7PvEIQ3!)E+uXUf#h_xVOSkSOPJ=^Qy+e|4q*+)|r z<(uN___B4XU6fdH!+~da6N{$@g*z2#fVF` z5LMuW`%9U^vo+O!CoE?{f!WeBps;M|7%-;a78d6YSx{IGT+HSPHWczT%hIhQ$kV#0 zaZTseQ)#>^@Tt-&4m<0-gNBNA7%x@HWqx;f+r9Gs6e>NLM_YCD`Ai%Gukv$B@l+QTg&GP%0^rlbCe%gvm?bN zAAUYS80$4)mYp^B@soH2ru&`Kn(6+hy4@t*+!7lBQ+32~>qZKg?!Kk&hGwVessyK# zeUsKPBbXO5VWbUOuboFlB#vpMpAcqed#arE#OlbIu_GcvEuXLFN*amAn0nOEst@3c z7FDZ~rWID}5VO?=;`em+iz1nm%_zbSWhMAt7`S^lMS^UYf;oB0W>8Bjb1f*6!M6m* zI(zn{e|}wIl7eK<%@EluuJnB{9a@$Vt3_#}t4 z?aT?RPkI?g7uJN}MpBPVo{n5isx{_v#pkb^}P{|1w4#A*-4C)yjLs%aO} z?QNzmzr7I?s9zhz=f6a=T;Q@w{LZ*00u&ZSShL8(ngu@LX8JS#g9U09NUzW$0eM)j zB$l7klLK?io)64F)Kk}pE0JkXhrXz7NH-%tjObJ;@!=U~<-=EN>9v}y%pS>PAs4RS z=$H4!r?K$4pP)w*8lm~djuxA6|8$H0a5FATmp9EHSC&BFs`vGlbd?sMym{9xKzS*~z(+_lZ99T8o zzs&BbnNZIs!?(D2e*eb3tY>3I{!?sjkpuk0h5CFm-_y$$>+q(Qh!%NBogD8vSHGh3FC~-)GsSrQiC)c` zQ*uoHhI$pGy#`fHZ42%J4YFlCUE1&I#ozjjnNWy$^j-gUM#y>JtO0LIfl>E)Ck_Ry zo~69X!gES(+%!a`y{!smhHGbw^8<63pB4r_COwdMqLp4~_3Ha%mXYg-PM*hU0_RPR z_6vp8kUeG!^Wu4&?Rz<^$5cP`F8SlU0$)|h*Hf_wZb8xKW*~Rh# zibu6+U&uaL!~8@i_kw8NJ~nNh&`>s#d=_j)T24m~7h$t9zI!oI>YV%5@y;@Sd*5|S zz4MJxvhGJrx(9bK~9|OjRFyV0JXZ*1bll z>jX=4JA;wkBg2=c`riv7_N&bH(l>@EKG8PeT+(3troEZvxTnK{)FKlQn^Sp1;m)iSl?B7@AH_@spqD?S4gX)Y^~#RWbWGV1uh zMBHV{oC};bHHZ)*);3x5SNdW1<(b?2mzim9%5%xzNV(as=r+@;`S8TAjLu~1-0Ddb z%#^OOVPMuF;rL#8V^X&#Lz?;m!Qk8 zQxH;Faq!`PYjAgE!F`Nf5>Lu@jqZEpoMm9jO(b=n-x#ylkefNELH?3K{!8=?a({OQ znM#MOzUtGLTG!ZV=G(Jmm+5#(Bm~{UF5Doiled)35#&-77?QLMU8K-{Ez?d|BYziX z|MwTotx7y6Z+FauYK|9@^RS9f2F%vpV2!Uo?OnB0c`e|lx_+yQ{HbcNrS&mJH0dY>MwekI`U*WFegm=`qjiiBu# zoINQuut(HtEkqaGY+#zdz!RjhdZ@`w>wYYr=SI8Xl>;#qxq(4_*IG%Jgjm9C&ti97 zNMpL9?D_qMk@%19Z^_SO+AJK?)b3ZzPkvAp8N095e0z=keLT}=xz9|85S-Vf9WDk8 zh8|>pC1WjV`7HjtBA-a3`SD9^LtTg6BoE(ZWW8chUMT)8pUtO8?y2+8IE7e1Oq$@@ zk3(}|7%Sd`8CC4V1x39#(f|G41{b1M}>WUp1%=tY4mn#Ob)o!;74g9s3s>FGAq{S*D5hwxzbChgPKuv+NXqU z5HOQ@AgRsh@9mdQWW>voTs|s9LweBWLlMmto|1=6bgkC-Z?%PjSURvGmTxS7!s=@} znL~$-8@L>LKUyS1jW$jrqhe~!vg4@heRMtWNprX*RZ#nzmsaMtNG2Hp_jY0r37K)( zpx1?KLn6Z)4&yQ|qgBME&(~|$yrR1)7QAwJIx2KekeCLpeI7C*d^2C#A?1771N&08 z*>UvM$sWw@+H_v)Bg+}`=^bCvzV)5I^4K}_?t!Y_02!`xX=vSBM4qI3ne-b-Xe`rv@Yj|F2 zG&XneWu}Pb=qs1Cy8<@6<721qte!aZ_EpS&>icUZ#V!M+;Q}M`sEh93JSw>oI$-<1vlT^J?Ol_)hgOSZcuok zjN?R|g{Fr|%4ce&G%~-*hEXS_MpN~Gwg6KSYtkaaRX^UlKX9smC$Hmgfg8Ld17Echh8NE329kM1 zS4U_oGn0jYXIR^_9To1PY>i%8+enBw=FCoP1f46@#I69(32tC|%NHtq3=Vg2+g;H! zDd8$Hy*_Sp{0wXPxo4B|24h#66;BmhxuRz;YS-10OoyFMRq@ zY~@ku{QY@vRy*~!$oj+!>N_>KAu<9(zRI*=LGN?O{&Sjf-)E$pbMBEQ6tdW+6QN&Vo@Y#~UM)5{k<;Au*nVLC;ADqT_Yu`C7 zTHr{hhI2Fb&)^#0slme;q3X@Tm|h;XFgzuCHu#8-w^V4SWHz-DQ%rsstyFf7utmV* zuY?Sh-N8Z}=^6Ey_c$QFQmsB2Nrj z;N^She0-xD!@?r0=9;H;=%dA3wV)Tr@1m9k&5FZL~;jBr(C|J+yj-2(bdpI z?}@cQ(*R{czVzy~v`Mky^w|l&`h|nTm zmLggw`B+d9Yr#;ekKHWKF>M+fN-iOGJ0l6pSk_a|bXs1o+uJPVlpMeLJ+kQAx3%JH z6EXsx1WGbFZM};IJShNnAV}987<11b|v78v!&zCo~^ke#x!qO4U$h-Ifw@oQIWA zt*L;OFyd=cY=l7qtEm^f{IP6}Wq$w$W*2HX<#S?o#0?gD6FGy6i zw`4k`84Y79t1uH{3+c=?0fZtorx2uwD~$t@uZX<~a)(1}F1KHrc(i=^ev^RcWXQ{- zdG?~MkuhT?qBCwUYy$T;GvKF0CI@(&5gAfKUT)IX*W$Y)#}_}Fr5^Aqm$dhNcGlCY zWvQC)PSmAMBWzuoKRmdcl{c9l%qdkVZ0WBTYp^lZ8RAQNr()yOQeSK!6WvN#$g_DH z-=`XRTG>N6i8XZ}V@kSoWhzX0wWvN#ighb+#Zo;~@xo3*#TnFnU0)t9C8$JKb}R-< zXP+?Yx`p_r-kxTpOqxJqE2?z-hxbXk0-@6IBMh0^^6rXPPdxf`>|V+$mulem(~mt+ zpW1=%eu~UzYCkIvsK!kdbh+qzU+33IBaCZV@iT8P@2vGWV#CWR-pguiPs^85*)Gv; zo})0Aom&?E`Kn83cnD9kV+gK8)510XBJ(oouN4yFnZ)X|l(lRjmq5<;{>oV==uudk zS^S@z-O&-Zw>g`j#i8g*)?m~3;NyJ179(rqDJnnkc$01JG1smcb~Oi2gS=M>tNC@e zlTT5d;S<%4sXJB3SM;jANA=#x_pfS9(1BG;_dM4)X=m-FXK!4y?p;}93JoL_57j3= zhMk8Xl(G^JzGvkr$T+X%?I-zay{QeU-Gh7m7>_RpA8zcj2}g;pqF5QC)_7kr(sSma zPe(NazViQS@2jJ#+P?j1q#LBWkyJtnIY_r4NP~nT(jd}ZDo9F5gCL#K-Q69BLrBA+ zQ{cD3d%3=Qzwch}AHOl)c;hh!W9-fNu+Bbf&Arx~YkuZuf%f!F9;=Y!!#)KYXBTg7 zfWgXZX0Q6v=byy|WX>aa5Ej1`BxWf(-NMn=%*$tbT1Pkf&hmI2^?N1biL|B;C>(@3+p6StbMe<#Fx(T@nxkCMqMPk84`h3FS zoZfQyB$O&*9iLLtm75l`u~pdEjEz$*ioeGA7Lh-AU97BTvvoE*Ajd?yJaQ6$DYP8F zw^j3@BZitWyh)6mAuEI1u{rKh-kOz>n=t1BbQQyANlki?-G*PnKCn4*F{DU6a2mHm zGedk*Gv-D7^n8ue6Ao#ftN;axWi=Ts8F|`Byhzva$lB5_MhD5Br?zKWCh87Dg+?yB zPVaUDhz!MmH#d-QW7VF1Xl&oE#kHjgjLRO5R#5RQWqK{ADod~wTYIRjN*%?fUb*l!s=N=*pMYO~|I5IF zgd&&ef+?b3wVwcI9AcEWJ{^Vy{tMAIaKk1Jr-3*Bs;|OHgE#5+3aVfxTp#|MP~@Rw z=VL*$T?dFMnV^EVzsZ&x-#(Z&`V}0WGvp&;YYHSCZcOg79b^A7g*LiZv4B{rSScMQ zi4K38XqUovh(Vp#A!&L+pW>0?+b_VV&R>se=k=&wS65p1-i#{Nk5TI1`Q9qsr(Onu)9>d~`c(Ici*{&(a|h@F*DLA&;Df&2_9J06K#EKUolr&1zb8L=#Cm`$EZ2Kt5}RDFU(e#7tnMj>bM2V_6kSM@}qg* zzi4@b3%T3qq)L)3aw#~PJeCtRdHU>*#(qac2p=02t?Sh4leP!X&?=@F(FArA3@2oy zz6qVjD)O?2^^i?ug(wQzOu_BOMbYN-t~jhMD#UImggfKdjGa_v#BGH@6~ZjzXRY^i zXVYeb!~4J?S{xmtWEBEbq`mnXD4D3|2f~hZl1skY_T5iWI*#{wglweSBNV|!9(3f1 zHc(5A8t(@{DGA`WF_!DPBJr#8S({D6LuL#hQIk;RO}60x?WR(BOSQiaM-@$TBC@Nf=@%S5Mgvx420gE4LERJUX4=hmg( zDDksZr)oE*rtqM8bV>a!2?bzmriya|x!+XmI^Q1u8J{8n0xS9;zTAX&H~Bp!I@ zr$N6xSRpEcCLC*3Z4aS?@v>-F&N(2WcrLDZF6_3i@9n>a8ST>d7Ck`Ij;cr?C2VT#`WfQ!e5D_+4$#iq0O%U;;&Lxbio@b3UX9lW5yR3%Y@F?TK2eg&lD;$e z^;IEdn@7UzMYty`wDw~j1*M3{{`yd`R%;JXuTXDy+UjQR&s0OUgR2iN6LS2r`POXx z@FE}M+`a0+|7hk6Di&-}emYAsxo|-}J6ksMad9d{EQ+#=vNc{$Lc5MJzCY}spBc*;75i3MkiY~bJ{=Y3evm&Z;Pg6y6vg#%yQ%SLT_iUa zMg;^m6KimJa?w!v$b*jHv~vmWR-J<{VIWUFw*poqc0CuV{o%K~DY&sBjO&uS8{QExA)&*;sQt*r z{8--2+vwD{`$543ZD;jr#-tpd4_bjga{kfe_RnnK^yL;vewbsiU zy{Z>E3?~n9Qi&Mu%ZR(mH}QN>Esv7R#c@pvd9Kt-6v?Lnr+`uJNDS}+%8($b7ne*(657=4Bmx5uJnG)uyQdmE zWVDtD&0KK26}Pnwt(5I^oEpOPT*MwWo%(8{2M=3v+O6R2N}MmjoG8uWOE;Z|ih_;w zUp2gd(d}!n3a*rz7dOsvOs!0QrHS%JL7*=y(!gS8%=;cfee}E8SVJ1&{n3Dubr;cA3^m#2pmyU6Fp#G zDjEKuz|2FD08!E*BkD^}cMo(m)Y6yEuipj*jHESvNHwst?@1583?SVJgzc_wM z7Y(u#pJeow^V9=U~J8n%UU>Oc^)fom8E^``C9uT3m(;=|@| zd%cL7G*Cz(&9%CH&!H$j2>F0ktQS3=D)zKDrg$t^xWqDUOQ-x|b3|A0PM(+%?<
8oor zMvpgE8|Jzs^Oeirp+R?AphwT5-?K3=V#g<^*QPH#P&D$1l50XaSCrKbv(r{%|5*KS z&G%EJ)uNA9=cdelsk*ST1xx<$`JI4TsE-4lF*Di2uQnk8y0(ql8t=vP)j?BH&I5=( zCr`_JsqjC)c$c4zGDFa(i#L`IDt`6TMu%M4B=gmNAzSi5N6^_ zyoU{eE6nOOD~QIuVfyD_&*of?Y-rrES{r>z6tRzSKxyr!qmkjx=J)sX4Wcep`3Er{ z=ex$_+F5$`?A;MlHp$kXXX4MEmR9F2IS9xV92e%FohIPvAL$IbWN0IBY+D4Mzz->W zDFN)|2c1pz4+qMX_Qqb&ucw@!A|)I+7}eicCBKuoKT>_UMZ=q3(9hZoD(t82Q$VC1 zyCp;AOs&CTvlG^{v?~@=yo6O=T#Wv7f?Q3rVnu|C5Y?+-%+!JAH4Q-N4l=?*BxN)b z-gAWZRL^LngM4@8DpD7@2COyvMTB-Afjv!bkH$u_s+)^q5kvl}6* zVOfg|jwFb|;T#Bf>moW#^7Z=2PFL(ofp&@>y=@V{DQ~H+5B+0EKHJd}Nnpey4*X zd%b%xj&v7x3?18(UF6W7vXbvFW9HGIeBo_R!Wkr4_wew@)4gh0gm&v@9rHl$hNw$| z#dZ!&vbXje=LF`Ajftw`Jf0xIBbNH=XmeX~&JUoIB_cBVH{)&BZYM7qV(qLlbAEJ* zYozRO`-+iB0l$FyrO*nv#}oCj+6@(m?POj2zI9CuoN=Dd%LoF zjEGR1kLEiG!&Q}epc!+tdHWrzLfG0SNbH*-BBD}!L8%4_6m0@a?&)~E2#7tBx=*hw z<<2R3js~_m!FRuM5$`HD3`)0ffP{0$9X9ZGHRg>OJ^jONcmKox!TeDye&e09VU`&I!>*BR)G^)kuH2!9DeDm- z)_OA*_@65H!xParnZ=P`kSY&ZYE3%X?1-hfF83#wCgyL&Y|62-=03910P}{mY)m+a zd~_4Q3J+v>btZ(d`CQ(44;|pRC3sH=p3J=0LcXjq0614{iv7dTGOywAi;V#sobvE)iVJ+o289SC}!mmQ*&t{p!W(Gw|U6;+hPcFH?P zO4;xr6$0WrE+=xCgdZJkc-{mFM?GE0}y4i9$A%<+mW zgXF?gSD>*uL*i^o({&N@v!X^aMlyyONAgdMyfT~{2DzWAav-ds>vnjYz>sMnqEYZG z!%OmXKCxnhK@)z{0Y608tl{^3gact4_FDlKrCsl~@WI~N^2Yv>P@=o}x7vz;I9h7wrMh-?aw}`YVdnFWeG_ z7V5Fh-b~SNtzC&!8KZOlaV-2F3S+HYVkGaOAxALXcF`m7$zzkEEAsQ$PG)RA;C_VU zC^_%`cg9rvjfA&^ka53Y7WY0p9>RQRN_tC(839yA@NmTo1N;6GjL#vG3OXjc(j=

u)LZ^$Y%s1v?Vp`fK<2L$;c72DtdG?P;CQrM^q11JVrr3rOHibJ!i_V= zd1;!gs2;V4x zXAy8e`5@;VuE-o(Ls)a)tbn^Y7Ze=G3Ly* zTW3{8C?69(k7k+3XAg5-1EW4-1!ch6ASH>MjKF-!WNT`=IL13B+~5&-9j9`+yA8>1 ztm=}Ihfnr!ZXM?CsmFre61*oNU<7o3HUi(re>H*|_j2>b2!526Dh@pYJ=CiX_8ps7= zlmzFxJ&+2U-40SO)MXaPiyYZds@VA7kHKYR?s=5)Ma-?}e&C=ak!2VG7#CK^ULP`_EjOp+!$3wJWfn5`cXQLIpVNFb5wY3v+ODLIzo-;`g@3c7_ zqF0}X5wD^4VP2eWcgHW;0D2II)sYawvo}^DP7{hE1H3hB=~X@c8gY0%q^XmXHJEDoC91ktM|Tw9Tokx zi?mI}PK~U4z4j-n>nlgr9eK8!Vo;kxAeMcbCylv{hymqK9m zw=T4btc*qdc!Y6xA=-G=M4o8nMtN_ftBvMF-q&Zr_PAQvqq(Hyj3g)OA@C)79rQsC z&8*fR953-7?;@5DV2rwXughDRBe|p@ralSl*$N6LYGPYp!Tz84}VS8zPYzTQz%mu2pkiP}m?HL|Y3%&F7=lOY(QXL>F)z z1Ib2!jj&(a$WME6{;zLrK_)Kd za>~l{7Embo9SrT+sE(h&DY2ksmWl((c<{U^{9rCkH6pB%^**9#Mu^s-v$!CYnGRRX z!Y}9pmBDk@GUZ7L-Y71s&2$)Bsi_O?AIrw1h{oOACI`ITIH2^=0Ak9l37Yx~MI9k^TC zylFVORO_cW<(rs7nP%*WVtQo~{e2R=$KW8ace1&fn-{-6*IGa-A;-}Z57Rba@}aYI z@gm`q$+-Kg%auJeVfq?BvVD(zRM&@&Q~k&z{=>Fw;v(oB{=(zSPmybGhtN}_g;v^3 zvo2cHe@|n=5Ga;c5I1Upe@x(X@3BHC8U^-{JYQXV+uFOt71_$m43Hj}MjgFkl!}jn z>7b^kc>npmZ)i3=G-K`|J7?KPD_zzzu&0lugeRW3QBms9P$DgQxkxxLIIdq(9 z*p!MYurgw}($ICKE|x2QVBY64R*u1n=J>K)4MJ^2#+`jP5I`}Tc$~2B-+kox&;muM zR*w(;+}MAa(v>$K*uVc*8@_<&xCIMSer+$m5!>qkdnpGgZZyv$+zIP^++WZK>IdmXTp?7w2G;>_sE=I41sjE z(FJ_;z;bD``<96wu!$BgDe8e1&6OA^sKj5hiI@%E&HdEHr?&1!WK!#)y``&ZxU140 zQ<=}e8r%Iw!RPKrgxjevTHmQu8L(@kokPYG(pPXzJY)1q?dB%IZS(UnktgdS;H~tj z&lNEmFiMUpnPA7u&57jHyWl|T`gC*mtq+GyC9Y0>fbl#Qms^ud=#~|s?NXN9;)hB^fp>rto6#5h?T(IUJ$TM@94maYTDhg7n5m$(mkv_W{)?6)U0{D zyx-wk5t&}XP-4bcMOx&Sk0(RsS+3R`8B_x8YJdwgE5LWeaKshlmV8%bjoUbS1vRK$ zQX-bgw)P3;+iz74s4XpPEYf^pnZvGSprwU~uqTq-U7d0-9c3-14O8qpJNMT5b?*`i z-^y;tBscE3M?z0>vPq;M?>^ZFPKd7K+FPF@EsOai`}}PHs$SWB|Kdar`@O8;UNX7; zaY|vwwQYNc@-yU1r=&V~2hc@^^M&j71oG>PDNmbKkz*@oh#M|twQbz(AU=^-v(=nY z9~?eJO)d*dST4JY4EWC>-1M$$>X0#Qnz))GWAP8+LgMhy_ zZvkPGF_tyGzqY|so&t+Go17Q}vD~(Hd3Y`O+ng@XMsDR1LNqCnoD0aC`s?%Zrut7A zQZKq6X3wtINN%-ki&F}h#fm*Sa=w=#;KiY5W>RiXeu&E&YuV5P>k{itiuAfwUiy?K}CQh3zFCE;Ej=j4_6-PfC!dkB{#c6eK)pk@Q!nVV+6iB#$kSKuXQi z{~u-j~*C*ZkE6!ZI7Tgxte0&#m_A;l0)(Blw(FDVc7wS8q- zIOWq0DFdrChi%;p=DS4Hq2&v^Ym^|6nOjzrnxbI^Zrlt$L}1VGmX|6ZwOOHLt5r}2285k{(w%ZUirY*`K5i* zrQ4YBQCB?jr|+(iE|^*h}gMrtQs z8#2NpUh;oKxp2O zG~Ytw+npg60Tb<4oV`=N9c0iW!?Ko)GO5GFFVSu~ytQjw!lc~7HeI79CmIV@E&PvW z7}G67$EqJ~QrHZucd$8rM0b9XUTPhZ)TG;3qeYw+cwiuc^VOC2<)!?%`fZw3fxkQp1;+S?eL1BWKq0)K{8yhc-8w|ZSw45xpCMzK)*p~HNQq`yYqzD8ba zNWtgdAdxJjEMNWj_+J*J2R{7rs!7D3Y7Ti`-86}82xt-mCxQTffK9so=)ckO2Hrcz z>$d)fNFD!6oloE^%3XiOpGNJ!m%|3o53kX5*JJVbI(-26)t{RdgWUA;pG4^ZG~0DQ z{*-?HkFsh&hyRME7y$ccr6>UODuLTe=LgvzTrQYy#NjU80w$;>Mzw10O+4- z@c`%-iO~R3{JO*bkp>GOe@)Tl2KdkF9RT##xD_|RzfYe7uz!&s1HM+`)gZr0Qz6--2muckKzF^!}Wvwe_`$XvWEwB-M``M z0`PxB0tMjzf>j8h{}o#ez)xLuA^cCA#Vbvc-~EUm?6%+cz5#T{HM;S7MY=vYAb#ow zjrJQC?(b+%S8e~_&;Rwn|LPw2H}4d%ZIJdm4>+_=)(c2kDE44eb8C80>1HLTQK#iL zvcC9I&t@t`cVQIMzcGj3k%5a%W-BbL&uguY! zjxF-wB}oM)iNQYVEN1?Yxr0ip+rzYE5OWp=++EKnLPKNm*(@uetI}l)+&=8R`cd-* z^IzD>fU~U1BHDF>5XOagHsq$zehTEo}Mcr9zGQHy}d^=NMKbo`~n{k)PxYy>P!;}0Tm9}tiOxu_ZR zLNRd)R2n+EPL+U15#O*8$yvRq5YU*h(cj-44p+iPFsMsve-}CUK}^ZrL+bJoc_jO5 z6b_wKM@w+)OGKSm-;?fY=tg&~213oP`(8^hq;Phgz7*W#^`p2!k~9sGCyy!Ww>P9YYWQ0;Q(}E7KoUMcnlaImgUg|Ui$YvLM95@GBAUpu z8x-VN(^pOM94rO1#KwM>&a|9_qtY7cMq=#SsOIE w5U(ToUd8DA>cE*hgjZJ}E8sWr>Plr821eEC>dF>i)7jBpT_w>t+>3zue}54El>h($ literal 0 HcmV?d00001 From 577c1626a94a711429a0ec51cd84b2856513b17c Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 6 Oct 2023 10:32:11 +0200 Subject: [PATCH 3/7] feat: Python bindings for geantino hypothesis (#2508) --- .../Acts/EventData/ParticleHypothesis.hpp | 6 ++++++ Examples/Python/src/EventData.cpp | 17 ++++++++++++++--- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/Core/include/Acts/EventData/ParticleHypothesis.hpp b/Core/include/Acts/EventData/ParticleHypothesis.hpp index 994f5179e85..3fb33b27664 100644 --- a/Core/include/Acts/EventData/ParticleHypothesis.hpp +++ b/Core/include/Acts/EventData/ParticleHypothesis.hpp @@ -110,6 +110,9 @@ class NonNeutralChargedParticleHypothesis pion().mass(), absQ); } + static NonNeutralChargedParticleHypothesis chargedGeantino() { + return chargedGeantino(Acts::UnitConstants::e); + } static NonNeutralChargedParticleHypothesis chargedGeantino(float absQ) { return NonNeutralChargedParticleHypothesis(PdgParticle::eInvalid, 0, absQ); } @@ -154,6 +157,9 @@ class ParticleHypothesis : public GenericParticleHypothesis { static ParticleHypothesis geantino() { return NeutralParticleHypothesis::geantino(); } + static ParticleHypothesis chargedGeantino() { + return chargedGeantino(Acts::UnitConstants::e); + } static ParticleHypothesis chargedGeantino(float absQ) { return ParticleHypothesis(PdgParticle::eInvalid, 0, absQ); } diff --git a/Examples/Python/src/EventData.cpp b/Examples/Python/src/EventData.cpp index ffbcaba41be..455b6c541cc 100644 --- a/Examples/Python/src/EventData.cpp +++ b/Examples/Python/src/EventData.cpp @@ -36,9 +36,20 @@ void addEventData(Context& ctx) { [](py::object /* self */) { return Acts::ParticleHypothesis::pion(); }) - .def_property_readonly_static("electron", [](py::object /* self */) { - return Acts::ParticleHypothesis::electron(); - }); + .def_property_readonly_static( + "electron", + [](py::object /* self */) { + return Acts::ParticleHypothesis::electron(); + }) + .def_property_readonly_static( + "geantino", + [](py::object /* self */) { + return Acts::ParticleHypothesis::geantino(); + }) + .def_property_readonly_static( + "chargedGeantino", [](py::object /* self */) { + return Acts::ParticleHypothesis::chargedGeantino(); + }); } } // namespace Acts::Python From 9f0b3e4367d2133e4f1cbf35771f18e321d4d7b3 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 6 Oct 2023 14:14:12 +0200 Subject: [PATCH 4/7] feat: Python bindings for particle hypothesis in `PropagationAlgorithm` (#2510) --- Examples/Python/src/Propagation.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Examples/Python/src/Propagation.cpp b/Examples/Python/src/Propagation.cpp index 15f04f74d6e..9876481d553 100644 --- a/Examples/Python/src/Propagation.cpp +++ b/Examples/Python/src/Propagation.cpp @@ -116,9 +116,9 @@ void addPropagation(Context& ctx) { propagatorImpl, randomNumberSvc, mode, sterileLogger, debugOutput, energyLoss, multipleScattering, recordMaterialInteractions, ntests, d0Sigma, z0Sigma, phiSigma, thetaSigma, qpSigma, tSigma, phiRange, - etaRange, ptRange, ptLoopers, maxStepSize, propagationStepCollection, - propagationMaterialCollection, covarianceTransport, covariances, - correlations); + etaRange, ptRange, particleHypothesis, ptLoopers, maxStepSize, + propagationStepCollection, propagationMaterialCollection, + covarianceTransport, covariances, correlations); py::class_>( From e1323f7bfb7bb8d185c4dccbe52efd2123519171 Mon Sep 17 00:00:00 2001 From: felix-russo <72298366+felix-russo@users.noreply.github.com> Date: Fri, 6 Oct 2023 16:17:20 +0200 Subject: [PATCH 5/7] feat: time vertex seeding (#2460) The `AdaptiveGridTrackDensity` is adapted to accommodate time vertex seeding. Unit tests are added. In the unit tests, the analytical maximum and width of a 3D Gaussian along the `z`-axis is used, see [Gaussian_Track_Density.pdf](https://github.com/acts-project/acts/files/12713019/Gaussian_Track_Density.pdf) for a derivation. --- .../AdaptiveGridDensityVertexFinder.hpp | 2 +- .../AdaptiveGridDensityVertexFinder.ipp | 20 +- .../Vertexing/AdaptiveGridTrackDensity.hpp | 204 +++++--- .../Vertexing/AdaptiveGridTrackDensity.ipp | 305 +++++++---- .../AdaptiveGridTrackDensityTests.cpp | 495 ++++++++++++------ 5 files changed, 676 insertions(+), 350 deletions(-) diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp index 4e0085517f1..4896c7c6757 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp @@ -77,7 +77,7 @@ class AdaptiveGridDensityVertexFinder { /// /// Only needed if cacheGridStateForTrackRemoval == true struct State { - // Map from z from the z bin values to the corresponding track density + // Map from the z bin values to the corresponding track density DensityMap mainDensityMap; // Map from input track to corresponding track density map diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp index bc0737d0ad7..f33cf2947a7 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp @@ -60,22 +60,22 @@ auto Acts::AdaptiveGridDensityVertexFinder::find( if (not state.mainDensityMap.empty()) { if (not m_cfg.estimateSeedWidth) { // Get z value of highest density bin - auto maxZres = m_cfg.gridDensity.getMaxZPosition(state.mainDensityMap); + auto maxZTRes = m_cfg.gridDensity.getMaxZTPosition(state.mainDensityMap); - if (!maxZres.ok()) { - return maxZres.error(); + if (!maxZTRes.ok()) { + return maxZTRes.error(); } - z = *maxZres; + z = (*maxZTRes).first; } else { // Get z value of highest density bin and width - auto maxZres = - m_cfg.gridDensity.getMaxZPositionAndWidth(state.mainDensityMap); + auto maxZTResAndWidth = + m_cfg.gridDensity.getMaxZTPositionAndWidth(state.mainDensityMap); - if (!maxZres.ok()) { - return maxZres.error(); + if (!maxZTResAndWidth.ok()) { + return maxZTResAndWidth.error(); } - z = (*maxZres).first; - width = (*maxZres).second; + z = (*maxZTResAndWidth).first.first; + width = (*maxZTResAndWidth).second; } } diff --git a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp index c6befda103f..31c53af0d64 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp @@ -13,41 +13,77 @@ #include +#include + namespace Acts { /// @class AdaptiveGridTrackDensity -/// @brief Implements a 1-dim density grid to be filled with -/// track Gaussian distributions. Each single track is modelled -/// as a 2-dim Gaussian distribution grid in the d0-z0 plane, -/// but only the overlap with the z-axis (i.e. a 1-dim density -/// vector) needs to be calculated. +/// @brief Implements a 1D (no time seeding) / 2D (time seeding) +/// grid that is filled with track densities. +/// Each track is modelled by a 2D / 3D Gaussian distribution in the +/// d0-z0 / d0-z0-t0 plane, which is evaluated at d0=0. Therefore, +/// each track effectively lives in 1D / 2D. /// The position of the highest track density (of either a single /// bin or the sum of a certain region) can be determined. /// Single tracks can be cached and removed from the overall density. -/// Unlike the GaussianGridTrackDensity, the overall density vector -/// grows adaptively with the tracks densities being added to the grid. +/// Unlike in the GaussianGridTrackDensity, the overall density map +/// grows adaptively when tracks densities are added to the grid. /// -/// @tparam trkGridSize The 2-dim grid size of a single track, i.e. -/// a single track is modelled as a (trkGridSize x trkGridSize) grid -/// in the d0-z0 plane. Note: trkGridSize has to be an odd value. -template +/// @tparam spatialTrkGridSize Number of bins per track in z direction +/// @tparam temporalTrkGridSize Number of bins per track in t direction +/// @note In total, a track is represented by a grid of size +/// spatialTrkGridSize * temporalTrkGridSize +template class AdaptiveGridTrackDensity { - // Assert odd trkGridSize - static_assert(trkGridSize % 2); + // Assert odd spatial and temporal track grid size + static_assert(spatialTrkGridSize % 2); + static_assert(temporalTrkGridSize % 2); public: - using DensityMap = std::unordered_map; + // The first (second) integer indicates the bin's z (t) position + using Bin = std::pair; + // Mapping between bins and track densities + using DensityMap = std::unordered_map>; + // Coordinates in the z-t plane; the t value will be set to 0 if time + // vertex seeding is disabled + using ZTPosition = std::pair; + // z-t position of a maximum and its width + using ZTPositionAndWidth = std::pair; /// The configuration struct struct Config { - /// @param binSize_ The binSize in mm - Config(float binSize_ = 0.1) : binSize(binSize_) {} - - // Z size of one single bin in grid - float binSize; // mm + /// @param spatialBinExtent_ The spatial extent of a bin in mm + Config(float spatialBinExtent_) : spatialBinExtent(spatialBinExtent_) { + if constexpr (temporalTrkGridSize > 1) { + throw std::invalid_argument( + "temporalBinExtent must be provided if temporalTrkGridSize > 1 " + "(i.e., if time vertex seeding is enabled)."); + } + } + + /// @param spatialBinExtent_ The spatial extent of a bin in mm + /// @param temporalBinExtent_ The temporal extent of a bin in mm + /// @note The speed of light is set to 1, hence the unit. + Config(float spatialBinExtent_, float temporalBinExtent_) + : spatialBinExtent(spatialBinExtent_), + temporalBinExtent(temporalBinExtent_) { + if constexpr (temporalTrkGridSize == 1) { + throw std::invalid_argument( + "temporalBinExtent must not be provided if temporalTrkGridSize == " + "1 (i.e., if time vertex seeding is disabled)."); + } + } + + // Spatial extent of a bin in d0 and z0 direction, should always be set to a + // positive value + float spatialBinExtent = 0.; // mm + + // Temporal extent of a bin, should be set to 0 if time vertex seeding is + // disabled (i.e., if temporalTrkGridSize = 1) + float temporalBinExtent = 0.; // mm // Do NOT use just the z-bin with the highest - // track density, but instead check the (up to) + // track density, but instead check (up to) // first three density maxima (only those that have // a maximum relative deviation of 'relativeDensityDev' // from the main maximum) and take the z-bin of the @@ -64,100 +100,120 @@ class AdaptiveGridTrackDensity { /// @brief Calculates the bin center from the bin number /// @param bin Bin number + /// @param binExtent Bin extent /// @return Bin center - float getBinCenter(int bin) const; + static float getBinCenter(int bin, float binExtent); - /// @brief Calculates the bin number corresponding to a d or z value - /// @param value d or z value + /// @brief Calculates the bin number corresponding to a d, z, or time value + /// @param value d, z, or time value + /// @param binExtent Bin extent /// @return Bin number - int getBin(float value) const; + static int getBin(float value, float binExtent); /// @brief Finds the maximum density of a DensityMap - /// @param densityMap Map between z bins and corresponding density value - /// @return Iterator of the map entry with the highest density value + /// @param densityMap Map between bins and corresponding density + /// values + /// @return Iterator of the map entry with the highest density DensityMap::const_iterator highestDensityEntry( const DensityMap& densityMap) const; - /// @brief Returns the z position of maximum (surrounding) track density + /// @brief Returns the z and t coordinate of maximum (surrounding) + /// track density + /// @note if time vertex seeding is not enabled, the t coordinate + /// will be set to 0. /// - /// @param densityMap Map from z bins to corresponding track density + /// @param densityMap Map between bins and corresponding density + /// values /// - /// @return The z position of maximum track density - Result getMaxZPosition(DensityMap& densityMap) const; + /// @return The z and t coordinates of maximum track density + Result getMaxZTPosition(DensityMap& densityMap) const; - /// @brief Returns the z position of maximum track density and - /// the estimated width + /// @brief Returns the z-t position of maximum track density + /// and the estimated z-width of the maximum /// - /// @param densityMap Map from z bins to corresponding track density + /// @param densityMap Map between bins and corresponding density + /// values /// - /// @return The z position of maximum track density and width - Result> getMaxZPositionAndWidth( + /// @return The z-t position of the maximum track density and + /// its width + Result getMaxZTPositionAndWidth( DensityMap& densityMap) const; /// @brief Adds a single track to the overall grid density /// - /// @param trk The track to be added. - /// @param mainDensityMap Map from z bins to corresponding track density. + /// @param trk The track to be added + /// @param mainDensityMap Map between bins and corresponding density /// /// @return The density map of the track that was added DensityMap addTrack(const BoundTrackParameters& trk, DensityMap& mainDensityMap) const; - /// @brief Removes a track from the overall grid density + /// @brief Removes a track from the overall grid density. /// - /// @param trackDensityMap Map from z bins to corresponding track density. - /// @note The track density comes from a single track. - /// @param mainDensityMap Map from z bins to corresponding track density. - /// @note The track density comes an arbitrary number of tracks. + /// @param trackDensityMap Map between bins and corresponding density + /// @note The track density comes from a single track + /// @param mainDensityMap Map between bins and corresponding density + /// @note The track density comes from an arbitrary number of tracks void subtractTrack(const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const; private: - /// @brief Function that creates a track density map, i.e., a map of z bins - /// to corresponding density values coming from a single track. - /// - /// @param d0 Transverse impact parameter - /// @param z0 Longitudinal impact parameter - /// @param centralZBin Central z bin of the track (where its density is the highest) - /// @param cov 2x2 impact parameter covariance matrix - DensityMap createTrackGrid(float d0, float z0, int centralZBin, - const Acts::SquareMatrix2& cov) const; - - /// @brief Function that estimates the seed width based on the full width - /// at half maximum (FWHM) of the maximum density peak + /// @brief Function that creates a track density map, i.e., a map from bins + /// to the corresponding density values for a single track. + /// + /// @param impactParams vector containing d0, z0, and t0 of the track + /// @param centralBin Central z and t bin of the track (where its + /// density is the highest) + /// @param cov 3x3 impact parameter covariance matrix + DensityMap createTrackGrid(const Acts::Vector3& impactParams, + const Bin& centralBin, + const Acts::SquareMatrix3& cov) const; + + /// @brief Function that estimates the seed width in z direction based + /// on the full width at half maximum (FWHM) of the maximum density peak /// @note This only works if the maximum is sufficiently isolated since /// overlapping neighboring peaks might lead to an overestimation of the /// seed width. /// - /// @param densityMap Map from z bins to corresponding track density - /// @param maxZ z-position of the maximum density value + /// @param densityMap Map from bins to corresponding track density + /// @param maxZT z-t position of the maximum density value /// /// @return The width Result estimateSeedWidth(const DensityMap& densityMap, - float maxZ) const; + const ZTPosition& maxZT) const; - /// @brief Helper to retrieve values according to a 2-dim normal distribution - /// @note This function is defined in coordinate system centered around d0 and z0 - float normal2D(float d, float z, const SquareMatrix2& cov) const; - - /// @brief Checks the (up to) first three density maxima (only those that have - /// a maximum relative deviation of 'relativeDensityDev' from the main - /// maximum) and take the z-bin of the maximum with the highest surrounding - /// density + /// @brief Helper to retrieve values of an nDim-dimensional normal + /// distribution + /// @note The constant prefactor (2 * pi)^(- nDim / 2) is discarded + /// + /// @param args Coordinates where the Gaussian should be evaluated + /// @note args must be in a coordinate system with origin at the mean + /// values of the Gaussian + /// @param cov Covariance matrix + /// + /// @return Multivariate Gaussian evaluated at args + template + static float multivariateGaussian(const Acts::ActsVector& args, + const Acts::ActsSquareMatrix& cov); + + /// @brief Checks (up to) first three density maxima that have a + /// maximum relative deviation of 'relativeDensityDev' from the + /// global maximum. Returns the bin of the maximum that has the + /// highest surrounding density in z direction. /// - /// @param densityMap Map between z bins and corresponding density value + /// @param densityMap Map between bins and corresponding density values /// - /// @return The z-bin position - int highestDensitySumBin(DensityMap& densityMap) const; + /// @return The bin corresponding to the highest surrounding density + Bin highestDensitySumBin(DensityMap& densityMap) const; - /// @brief Calculates the density sum of a z-bin and its two neighboring bins - /// as needed for 'highestDensitySumBin' + /// @brief Calculates the density sum of a bin and its two neighboring bins + /// in z direction /// - /// @param densityMap Map between z bins and corresponding density value - /// @param zBin The center z-bin whose neighbors we want to sum up + /// @param densityMap Map between bins and corresponding density values + /// @param bin Bin whose neighbors in z we want to sum up /// - /// @return The sum - float getDensitySum(const DensityMap& densityMap, int zBin) const; + /// @return The density sum + float getDensitySum(const DensityMap& densityMap, const Bin& bin) const; Config m_cfg; }; diff --git a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp index 8f527132213..8f106eb5ac5 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp @@ -10,20 +10,25 @@ #include -template -float Acts::AdaptiveGridTrackDensity::getBinCenter(int bin) const { - return bin * m_cfg.binSize; +template +float Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::getBinCenter(int bin, + float binExtent) { + return bin * binExtent; } -template -int Acts::AdaptiveGridTrackDensity::getBin(float value) const { - return static_cast(std::floor(value / m_cfg.binSize - 0.5) + 1); +template +int Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::getBin(float value, + float binExtent) { + return static_cast(std::floor(value / binExtent - 0.5) + 1); } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap::const_iterator -Acts::AdaptiveGridTrackDensity::highestDensityEntry( - const DensityMap& densityMap) const { +template +typename Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::DensityMap::const_iterator +Acts::AdaptiveGridTrackDensity:: + highestDensityEntry(const DensityMap& densityMap) const { auto maxEntry = std::max_element( std::begin(densityMap), std::end(densityMap), [](const auto& densityEntry1, const auto& densityEntry2) { @@ -32,118 +37,173 @@ Acts::AdaptiveGridTrackDensity::highestDensityEntry( return maxEntry; } -template -Acts::Result -Acts::AdaptiveGridTrackDensity::getMaxZPosition( - DensityMap& densityMap) const { +template +Acts::Result::ZTPosition> +Acts::AdaptiveGridTrackDensity:: + getMaxZTPosition(DensityMap& densityMap) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } - int zBin = -1; + Bin bin; if (!m_cfg.useHighestSumZPosition) { - zBin = highestDensityEntry(densityMap)->first; + bin = highestDensityEntry(densityMap)->first; } else { // Get z position with highest density sum // of surrounding bins - zBin = highestDensitySumBin(densityMap); + bin = highestDensitySumBin(densityMap); } // Derive corresponding z value - return getBinCenter(zBin); + float maxZ = getBinCenter(bin.first, m_cfg.spatialBinExtent); + + ZTPosition maxValues = std::make_pair(maxZ, 0.); + + // Get t value of the maximum if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + float maxT = getBinCenter(bin.second, m_cfg.temporalBinExtent); + maxValues.second = maxT; + } + + return maxValues; } -template -Acts::Result> -Acts::AdaptiveGridTrackDensity::getMaxZPositionAndWidth( - DensityMap& densityMap) const { - // Get z maximum value - auto maxZRes = getMaxZPosition(densityMap); - if (not maxZRes.ok()) { - return maxZRes.error(); +template +Acts::Result::ZTPositionAndWidth> +Acts::AdaptiveGridTrackDensity:: + getMaxZTPositionAndWidth(DensityMap& densityMap) const { + // Get z value where the density is the highest + auto maxZTRes = getMaxZTPosition(densityMap); + if (not maxZTRes.ok()) { + return maxZTRes.error(); } - float maxZ = *maxZRes; + ZTPosition maxZT = *maxZTRes; // Get seed width estimate - auto widthRes = estimateSeedWidth(densityMap, maxZ); + auto widthRes = estimateSeedWidth(densityMap, maxZT); if (not widthRes.ok()) { return widthRes.error(); } float width = *widthRes; - std::pair returnPair{maxZ, width}; - return returnPair; + ZTPositionAndWidth maxZTAndWidth{maxZT, width}; + return maxZTAndWidth; } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::addTrack( - const Acts::BoundTrackParameters& trk, DensityMap& mainDensityMap) const { - SquareMatrix2 cov = trk.spatialImpactParameterCovariance().value(); - float d0 = trk.parameters()[0]; - float z0 = trk.parameters()[1]; +template +typename Acts::AdaptiveGridTrackDensity::DensityMap +Acts::AdaptiveGridTrackDensity:: + addTrack(const Acts::BoundTrackParameters& trk, + DensityMap& mainDensityMap) const { + ActsVector<3> impactParams = trk.impactParameters(); + ActsSquareMatrix<3> cov = trk.impactParameterCovariance().value(); // Calculate bin in d direction - int centralDBin = getBin(d0); + int centralDBin = getBin(impactParams(0), m_cfg.spatialBinExtent); // Check if current track affects grid density - if (std::abs(centralDBin) > (trkGridSize - 1) / 2.) { + if (std::abs(centralDBin) > (spatialTrkGridSize - 1) / 2.) { DensityMap emptyTrackDensityMap; return emptyTrackDensityMap; } // Calculate bin in z direction - int centralZBin = getBin(z0); + int centralZBin = getBin(impactParams(1), m_cfg.spatialBinExtent); + + // If we don't do time vertex seeding, the time index is set to 0 + Bin centralBin = std::make_pair(centralZBin, 0.); + + // Calculate bin in t direction if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + int centralTBin = getBin(impactParams(2), m_cfg.temporalBinExtent); + centralBin.second = centralTBin; + } - DensityMap trackDensityMap = createTrackGrid(d0, z0, centralZBin, cov); + DensityMap trackDensityMap = createTrackGrid(impactParams, centralBin, cov); for (const auto& densityEntry : trackDensityMap) { - int zBin = densityEntry.first; + Bin bin = densityEntry.first; float trackDensity = densityEntry.second; // Check if z bin is already part of the main grid - if (mainDensityMap.count(zBin) == 1) { - mainDensityMap.at(zBin) += trackDensity; + if (mainDensityMap.count(bin) == 1) { + mainDensityMap.at(bin) += trackDensity; } else { - mainDensityMap[zBin] = trackDensity; + mainDensityMap[bin] = trackDensity; } } return trackDensityMap; } -template -void Acts::AdaptiveGridTrackDensity::subtractTrack( - const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const { +template +void Acts::AdaptiveGridTrackDensity:: + subtractTrack(const DensityMap& trackDensityMap, + DensityMap& mainDensityMap) const { for (auto it = trackDensityMap.begin(); it != trackDensityMap.end(); it++) { mainDensityMap.at(it->first) -= it->second; } } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::createTrackGrid( - float d0, float z0, int centralZBin, const Acts::SquareMatrix2& cov) const { +template +typename Acts::AdaptiveGridTrackDensity::DensityMap +Acts::AdaptiveGridTrackDensity:: + createTrackGrid(const Acts::Vector3& impactParams, const Bin& centralBin, + const Acts::SquareMatrix3& cov) const { DensityMap trackDensityMap; - int halfTrkGridSize = (trkGridSize - 1) / 2; - int firstZBin = centralZBin - halfTrkGridSize; - // Loop over columns - for (int j = 0; j < trkGridSize; j++) { - int zBin = firstZBin + j; - float z = getBinCenter(zBin); - trackDensityMap[zBin] = normal2D(-d0, z - z0, cov); + int halfSpatialTrkGridSize = (spatialTrkGridSize - 1) / 2; + int firstZBin = centralBin.first - halfSpatialTrkGridSize; + + // If we don't do time vertex seeding, firstTBin will be 0. + int halfTemporalTrkGridSize = (temporalTrkGridSize - 1) / 2; + int firstTBin = centralBin.second - halfTemporalTrkGridSize; + + // Loop over bins + for (int i = 0; i < temporalTrkGridSize; i++) { + int tBin = firstTBin + i; + // If we don't do vertex time seeding, we set the time to 0 since it will be + // discarded in the for loop below anyways + float t = 0; + if constexpr (temporalTrkGridSize > 1) { + t = getBinCenter(tBin, m_cfg.temporalBinExtent); + } + for (int j = 0; j < spatialTrkGridSize; j++) { + int zBin = firstZBin + j; + float z = getBinCenter(zBin, m_cfg.spatialBinExtent); + // Bin coordinates in the d-z-t plane + Acts::Vector3 binCoords(0., z, t); + // Transformation to coordinate system with origin at the track center + binCoords -= impactParams; + Bin bin = std::make_pair(zBin, tBin); + if constexpr (temporalTrkGridSize == 1) { + trackDensityMap[bin] = multivariateGaussian<2>( + binCoords.head<2>(), cov.topLeftCorner<2, 2>()); + } else { + trackDensityMap[bin] = multivariateGaussian<3>(binCoords, cov); + } + } } return trackDensityMap; } -template -Acts::Result -Acts::AdaptiveGridTrackDensity::estimateSeedWidth( - const DensityMap& densityMap, float maxZ) const { +template +Acts::Result Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, + temporalTrkGridSize>::estimateSeedWidth(const DensityMap& densityMap, + const ZTPosition& maxZT) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } - // Get z bin of max density and max density value - int zMaxBin = getBin(maxZ); - const float maxValue = densityMap.at(zMaxBin); + // Get z bin of max density + int zMaxBin = getBin(maxZT.first, m_cfg.spatialBinExtent); + int tMaxBin = 0; + // Fill the time bin with a non-zero value if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + tMaxBin = getBin(maxZT.second, m_cfg.temporalBinExtent); + } + const float maxValue = densityMap.at(std::make_pair(zMaxBin, tMaxBin)); int rhmBin = zMaxBin; float gridValue = maxValue; @@ -152,21 +212,21 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( bool binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(rhmBin + 1) == 0) { + if (densityMap.count(std::make_pair(rhmBin + 1, tMaxBin)) == 0) { binFilled = false; break; } rhmBin += 1; - gridValue = densityMap.at(rhmBin); + gridValue = densityMap.at(std::make_pair(rhmBin, tMaxBin)); } // Use linear approximation to find better z value for FWHM between bins float rightDensity = 0; if (binFilled) { - rightDensity = densityMap.at(rhmBin); + rightDensity = densityMap.at(std::make_pair(rhmBin, tMaxBin)); } - float leftDensity = densityMap.at(rhmBin - 1); - float deltaZ1 = m_cfg.binSize * (maxValue / 2 - leftDensity) / + float leftDensity = densityMap.at(std::make_pair(rhmBin - 1, tMaxBin)); + float deltaZ1 = m_cfg.spatialBinExtent * (maxValue / 2 - leftDensity) / (rightDensity - leftDensity); int lhmBin = zMaxBin; @@ -174,25 +234,25 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(lhmBin - 1) == 0) { + if (densityMap.count(std::make_pair(lhmBin - 1, tMaxBin)) == 0) { binFilled = false; break; } lhmBin -= 1; - gridValue = densityMap.at(lhmBin); + gridValue = densityMap.at(std::make_pair(lhmBin, tMaxBin)); } // Use linear approximation to find better z value for FWHM between bins - rightDensity = densityMap.at(lhmBin + 1); + rightDensity = densityMap.at(std::make_pair(lhmBin + 1, tMaxBin)); if (binFilled) { - leftDensity = densityMap.at(lhmBin); + leftDensity = densityMap.at(std::make_pair(lhmBin, tMaxBin)); } else { leftDensity = 0; } - float deltaZ2 = m_cfg.binSize * (rightDensity - maxValue / 2) / + float deltaZ2 = m_cfg.spatialBinExtent * (rightDensity - maxValue / 2) / (rightDensity - leftDensity); - float fwhm = (rhmBin - lhmBin) * m_cfg.binSize - deltaZ1 - deltaZ2; + float fwhm = (rhmBin - lhmBin) * m_cfg.spatialBinExtent - deltaZ1 - deltaZ2; // FWHM = 2.355 * sigma float width = fwhm / 2.355f; @@ -200,76 +260,89 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( return std::isnormal(width) ? width : 0.0f; } -template -float Acts::AdaptiveGridTrackDensity::normal2D( - float d, float z, const Acts::SquareMatrix2& cov) const { - float det = cov.determinant(); - float coef = 1 / std::sqrt(det); - float expo = - -1 / (2 * det) * - (cov(1, 1) * d * d - (cov(0, 1) + cov(1, 0)) * d * z + cov(0, 0) * z * z); - return coef * safeExp(expo); +template +template +float Acts::AdaptiveGridTrackDensity:: + multivariateGaussian(const Acts::ActsVector& args, + const Acts::ActsSquareMatrix& cov) { + float expo = -0.5 * args.transpose().dot(cov.inverse() * args); + float gaussianDensity = safeExp(expo) / std::sqrt(cov.determinant()); + return gaussianDensity; } -template -int Acts::AdaptiveGridTrackDensity::highestDensitySumBin( - DensityMap& densityMap) const { - // Checks the first (up to) 3 density maxima, if they are close, checks which - // one has the highest surrounding density sum (the two neighboring bins) - +template +typename Acts::AdaptiveGridTrackDensity::Bin +Acts::AdaptiveGridTrackDensity:: + highestDensitySumBin(DensityMap& densityMap) const { // The global maximum auto firstMax = highestDensityEntry(densityMap); - int zBinFirstMax = firstMax->first; + Bin binFirstMax = firstMax->first; float valueFirstMax = firstMax->second; - float firstSum = getDensitySum(densityMap, zBinFirstMax); + float firstSum = getDensitySum(densityMap, binFirstMax); + // Smaller maxima must have a density of at least: + // valueFirstMax - densityDeviation float densityDeviation = valueFirstMax * m_cfg.maxRelativeDensityDev; // Get the second highest maximum - densityMap.at(zBinFirstMax) = 0; + densityMap.at(binFirstMax) = 0; auto secondMax = highestDensityEntry(densityMap); - int zBinSecondMax = secondMax->first; + Bin binSecondMax = secondMax->first; float valueSecondMax = secondMax->second; float secondSum = 0; if (valueFirstMax - valueSecondMax < densityDeviation) { - secondSum = getDensitySum(densityMap, zBinSecondMax); + secondSum = getDensitySum(densityMap, binSecondMax); + } else { + // If the second maximum is not sufficiently large the third maximum won't + // be either + densityMap.at(binFirstMax) = valueFirstMax; + return binFirstMax; } // Get the third highest maximum - densityMap.at(zBinSecondMax) = 0; + densityMap.at(binSecondMax) = 0; auto thirdMax = highestDensityEntry(densityMap); - int zBinThirdMax = thirdMax->first; + Bin binThirdMax = thirdMax->first; float valueThirdMax = thirdMax->second; float thirdSum = 0; if (valueFirstMax - valueThirdMax < densityDeviation) { - thirdSum = getDensitySum(densityMap, zBinThirdMax); + thirdSum = getDensitySum(densityMap, binThirdMax); } // Revert back to original values - densityMap.at(zBinFirstMax) = valueFirstMax; - densityMap.at(zBinSecondMax) = valueSecondMax; + densityMap.at(binFirstMax) = valueFirstMax; + densityMap.at(binSecondMax) = valueSecondMax; - // Return the z-bin position of the highest density sum + // Return the z bin position of the highest density sum if (secondSum > firstSum && secondSum > thirdSum) { - return zBinSecondMax; + return binSecondMax; } if (thirdSum > secondSum && thirdSum > firstSum) { - return zBinThirdMax; + return binThirdMax; } - return zBinFirstMax; + return binFirstMax; } -template -float Acts::AdaptiveGridTrackDensity::getDensitySum( - const DensityMap& densityMap, int zBin) const { - float sum = densityMap.at(zBin); +template +float Acts::AdaptiveGridTrackDensity:: + getDensitySum(const DensityMap& densityMap, const Bin& bin) const { + // Add density from the bin. + float sum = densityMap.at(bin); // Check if neighboring bins are part of the densityMap and add them (if they - // are not part of the map, we assume them to be 0) Note that each key in a - // map is unique; the .count() function therefore returns either 0 or 1 - if (densityMap.count(zBin - 1) == 1) { - sum += densityMap.at(zBin - 1); + // are not part of the map, we assume them to be 0). + // Note that each key in a map is unique; the .count() function therefore + // returns either 0 or 1. + Bin binShifted = bin; + // Add density from the neighboring bin in -z direction. + binShifted.first -= 1; + if (densityMap.count(binShifted) == 1) { + sum += densityMap.at(binShifted); } - if (densityMap.count(zBin + 1) == 1) { - sum += densityMap.at(zBin + 1); + + // Add density from the neighboring bin in +z direction. + binShifted.first += 2; + if (densityMap.count(binShifted) == 1) { + sum += densityMap.at(binShifted); } return sum; } diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp index 2b33c5d789c..87f2e141070 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp @@ -37,21 +37,30 @@ namespace Test { using Covariance = BoundSquareMatrix; +Covariance makeRandomCovariance(int seed = 31415) { + std::srand(seed); + Covariance randMat((Covariance::Random() + 1.5 * Covariance::Identity()) * + 0.05); + + // symmetric covariance matrix + Covariance covMat = 0.5 * (randMat + randMat.transpose()); + + return covMat; +} + BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { using Vector2 = Eigen::Matrix; using Matrix2 = Eigen::Matrix; // Using a large track grid so we can choose a small bin size - const int trkGridSize = 4001; + const int spatialTrkGridSize = 4001; // Arbitrary (but small) bin size - const float binSize = 3.1e-4; + const float binExtent = 3.1e-4; // Arbitrary impact parameters const float d0 = 0.4; const float z0 = -0.2; Vector2 impactParameters{d0, z0}; - Covariance covMat(Covariance::Identity() * 0.05); - covMat(0, 1) = -0.02; - covMat(1, 0) = -0.02; + Covariance covMat = makeRandomCovariance(); Matrix2 subCovMat = covMat.block<2, 2>(0, 0).cast(); BoundVector paramVec; paramVec << d0, z0, 0, 0, 0, 0; @@ -63,11 +72,11 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { BoundTrackParameters params1(perigeeSurface, paramVec, covMat, ParticleHypothesis::pion()); - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Add track auto trackDensityMap = grid.addTrack(params1, mainDensityMap); @@ -85,10 +94,10 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { for (auto const& it : mainDensityMap) { // Extract variables for better readability - int zBin = it.first; + int zBin = it.first.first; float density = it.second; // Argument for 2D gaussian - Vector2 dzVec{0., grid.getBinCenter(zBin)}; + Vector2 dzVec{0., grid.getBinCenter(zBin, binExtent)}; // Compute correct density... float correctDensity = gaussian2D(dzVec, impactParameters, subCovMat); // ... and check if our result is equivalent @@ -106,11 +115,11 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { subCovMat.determinant() / subCovMat(0, 0)); // Estimate maximum z position and seed width - auto res = grid.getMaxZPositionAndWidth(mainDensityMap); + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); BOOST_CHECK(res.ok()); // Extract variables for better readability... - float maxZ = res.value().first; + float maxZ = res.value().first.first; float fwhm = res.value().second * 2.355f; // ... and check if they are correct (note: the optimization is not as exact @@ -120,15 +129,116 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); } -BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { +BOOST_AUTO_TEST_CASE( + compare_to_analytical_solution_for_single_track_with_time) { + // Number of bins in z- and t-direction + const int spatialTrkGridSize = 401; + const int temporalTrkGridSize = 401; + // Bin extents + const float spatialBinExtent = 3.1e-3; + const float temporalBinExtent = 3.1e-3; + // Arbitrary impact parameters + const float d0 = -0.1; + const float z0 = -0.2; + const float t0 = 0.1; + Vector3 impactParameters{d0, z0, t0}; + + // symmetric covariance matrix + Covariance covMat = makeRandomCovariance(); + + BoundVector paramVec; + paramVec << d0, z0, 0, 0, 0, t0; + // Create perigee surface + std::shared_ptr perigeeSurface = + Surface::makeShared(Vector3(0., 0., 0.)); + + BoundTrackParameters params(perigeeSurface, paramVec, covMat, + ParticleHypothesis::pion()); + + ActsSquareMatrix<3> ipCov = params.impactParameterCovariance().value(); + + AdaptiveGridTrackDensity::Config cfg( + spatialBinExtent, temporalBinExtent); + AdaptiveGridTrackDensity grid(cfg); + + // Empty map + AdaptiveGridTrackDensity::DensityMap + mainDensityMap; + + // Add track + auto trackDensityMap = grid.addTrack(params, mainDensityMap); + + float relTol = 1e-5; + float small = 1e-5; + + auto gaussian3D = [&](const Vector3& args, const Vector3& mus, + const SquareMatrix3& sigmas) { + Vector3 diffs = args - mus; + float coef = 1 / std::sqrt(sigmas.determinant()); + float expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); + return coef * std::exp(expo); + }; + + for (auto const& it : mainDensityMap) { + // Extract variables for better readability + float z = grid.getBinCenter(it.first.first, spatialBinExtent); + float t = grid.getBinCenter(it.first.second, temporalBinExtent); + float density = it.second; + // Argument for 3D gaussian + Vector3 dztVec{0., z, t}; + + // Compute correct density... + float correctDensity = gaussian3D(dztVec, impactParameters, ipCov); + + // ... and check if our result is equivalent + CHECK_CLOSE_OR_SMALL(density, correctDensity, relTol, small); + } + + // The analytical calculations of the following can be found here: + // https://github.com/acts-project/acts/pull/2460. + // TODO: upload reference at a better place. + // Analytical maximum of the Gaussian + ActsSquareMatrix<3> ipWeights = ipCov.inverse(); + ActsScalar denom = + ipWeights(1, 1) * ipWeights(2, 2) - ipWeights(1, 2) * ipWeights(1, 2); + + ActsScalar zNom = + ipWeights(0, 1) * ipWeights(2, 2) - ipWeights(0, 2) * ipWeights(1, 2); + ActsScalar correctMaxZ = zNom / denom * d0 + z0; + + ActsScalar tNom = + ipWeights(0, 2) * ipWeights(1, 1) - ipWeights(0, 1) * ipWeights(1, 2); + ActsScalar correctMaxT = tNom / denom * d0 + t0; + + // Analytical FWHM of the Gaussian + ActsScalar correctFWHM = 2. * std::sqrt(2 * std::log(2.) / ipWeights(1, 1)); + + // Estimate maximum z position and seed width + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); + BOOST_CHECK(res.ok()); + + // Extract variables for better readability... + float maxZ = res.value().first.first; + float maxT = res.value().first.second; + float fwhm = res.value().second * 2.355f; + + // ... and check if they are correct (note: the optimization is not as exact + // as the density values). + float relTolOptimization = 1e-1; + CHECK_CLOSE_REL(maxZ, correctMaxZ, relTolOptimization); + CHECK_CLOSE_REL(maxT, correctMaxT, relTolOptimization); + CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); +} + +BOOST_AUTO_TEST_CASE(seed_width_estimation) { // Dummy track grid size (not needed for this unit test) - const int trkGridSize = 1; - float binSize = 2.; - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + const int spatialTrkGridSize = 1; + float binExtent = 2.; + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // z-position of the maximum density float correctMaxZ = -2.; @@ -138,16 +248,16 @@ BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { // of 20 mm. The linear approximation we use during the seed width estimation // should be exact in this case. for (int i = -6; i <= 4; i++) { - mainDensityMap[i] = - 1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i)); + mainDensityMap[std::make_pair(i, 0)] = + 1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i, binExtent)); } // Get maximum z position and corresponding seed width - auto res = grid.getMaxZPositionAndWidth(mainDensityMap); + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); BOOST_CHECK(res.ok()); // Check if we found the correct maximum - float maxZ = res.value().first; + float maxZ = res.value().first.first; BOOST_CHECK_EQUAL(maxZ, correctMaxZ); // Calculate full width at half maximum from seed width and check if it's @@ -156,13 +266,13 @@ BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { BOOST_CHECK_EQUAL(fwhm, 10.); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { - const int trkGridSize = 15; +BOOST_AUTO_TEST_CASE(track_adding) { + const int spatialTrkGridSize = 15; - double binSize = 0.1; // mm + double binExtent = 0.1; // mm - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Create some test tracks in such a way that some tracks // e.g. overlap and that certain tracks need to be inserted @@ -192,7 +302,7 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { ParticleHypothesis::pion()); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Track is too far away from z axis and was not added auto trackDensityMap = grid.addTrack(params0, mainDensityMap); @@ -200,38 +310,49 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); + BOOST_CHECK_EQUAL(mainDensityMap.size(), spatialTrkGridSize); // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params2, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 2 * trkGridSize); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 2 * spatialTrkGridSize); // Track 3 has overlap of 2 bins with track 1 trackDensityMap = grid.addTrack(params3, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * trkGridSize - 2); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); // Add first track again, should *not* introduce new z entries trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * trkGridSize - 2); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_max_z_and_width_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(max_z_t_and_width) { + const int spatialTrkGridSize = 29; + const int temporalTrkGridSize = 29; + + // spatial and temporal bin extent + double binExtent = 0.05; - double binSize = 0.05; // mm + // 1D grid of z values + AdaptiveGridTrackDensity::Config cfg1D(binExtent); + AdaptiveGridTrackDensity grid1D(cfg1D); - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + // 2D grid of z and t values + AdaptiveGridTrackDensity::Config + cfg2D(binExtent, binExtent); + AdaptiveGridTrackDensity grid2D( + cfg2D); // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); float z0Trk1 = 0.25; + float t0Trk1 = 0.05; float z0Trk2 = -10.95; + float t0Trk2 = 0.1; BoundVector paramVec1; - paramVec1 << 0.02, z0Trk1, 0, 0, 0, 0; + paramVec1 << 0.02, z0Trk1, 0, 0, 0, t0Trk1; BoundVector paramVec2; - paramVec2 << 0.01, z0Trk2, 0, 0, 0, 0; + paramVec2 << 0.01, z0Trk2, 0, 0, 0, t0Trk2; // Create perigee surface std::shared_ptr perigeeSurface = @@ -242,37 +363,65 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_max_z_and_width_test) { BoundTrackParameters params2(perigeeSurface, paramVec2, covMat, ParticleHypothesis::pion()); - // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; - - auto trackDensityMap = grid.addTrack(params1, mainDensityMap); - auto res1 = grid.getMaxZPosition(mainDensityMap); - BOOST_CHECK(res1.ok()); - // Maximum should be at z0Trk1 position - BOOST_CHECK_EQUAL(*res1, z0Trk1); - - trackDensityMap = grid.addTrack(params2, mainDensityMap); - auto res2 = grid.getMaxZPosition(mainDensityMap); - BOOST_CHECK(res2.ok()); - // Trk 2 is closer to z-axis and should yield higher density values - // New maximum is therefore at z0Trk2 - BOOST_CHECK_EQUAL(*res2, z0Trk2); - - auto resWidth1 = grid.getMaxZPositionAndWidth(mainDensityMap); - BOOST_CHECK(resWidth1.ok()); - BOOST_CHECK_EQUAL((*resWidth1).first, z0Trk2); - BOOST_CHECK((*resWidth1).second > 0); + // Empty maps + AdaptiveGridTrackDensity::DensityMap mainDensityMap1D; + AdaptiveGridTrackDensity::DensityMap + mainDensityMap2D; + + // Add first track to spatial grid + auto trackDensityMap = grid1D.addTrack(params1, mainDensityMap1D); + auto firstRes = grid1D.getMaxZTPosition(mainDensityMap1D); + BOOST_CHECK(firstRes.ok()); + // Maximum should be at z0Trk1 position ... + BOOST_CHECK_EQUAL((*firstRes).first, z0Trk1); + // ... and the corresponding time should be set to 0 + BOOST_CHECK_EQUAL((*firstRes).second, 0.); + + // Add first track to 2D grid + auto trackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); + auto firstRes2D = grid2D.getMaxZTPosition(mainDensityMap2D); + BOOST_CHECK(firstRes2D.ok()); + // Maximum should be at z0Trk1 position ... + BOOST_CHECK_EQUAL((*firstRes2D).first, z0Trk1); + // ... and the corresponding time should be at t0Trk1 + BOOST_CHECK_EQUAL((*firstRes2D).second, t0Trk1); + + // Add second track to spatial grid + trackDensityMap = grid1D.addTrack(params2, mainDensityMap1D); + // Calculate maximum and the corresponding width + auto secondRes = grid1D.getMaxZTPositionAndWidth(mainDensityMap1D); + BOOST_CHECK(secondRes.ok()); + // Trk 2 is closer to z-axis and should thus yield higher density values. + // Therefore, the new maximum is at z0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes).first.first, z0Trk2); + // ... the corresponding time should be set to 0... + BOOST_CHECK_EQUAL((*secondRes).first.second, 0.); + // ... and it should have a positive width + BOOST_CHECK((*secondRes).second > 0); + + // Add second track to 2D grid + trackDensityMap = grid2D.addTrack(params2, mainDensityMap2D); + // Calculate maximum and the corresponding width + auto secondRes2D = grid2D.getMaxZTPositionAndWidth(mainDensityMap2D); + BOOST_CHECK(secondRes2D.ok()); + // Trk 2 is closer to z-axis and should thus yield higher density values. + // Therefore, the new maximum is at z0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes2D).first.first, z0Trk2); + // ... the corresponding time should be at t0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes2D).first.second, t0Trk2); + // ... and it should have approximately the same width in z direction + CHECK_CLOSE_OR_SMALL((*secondRes2D).second, (*secondRes).second, 1e-5, 1e-5); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_highest_density_sum_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(highest_density_sum) { + const int spatialTrkGridSize = 29; - double binSize = 0.05; // mm + double binExtent = 0.05; // mm - AdaptiveGridTrackDensity::Config cfg(binSize); + AdaptiveGridTrackDensity::Config cfg(binExtent); cfg.useHighestSumZPosition = true; - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity grid(cfg); // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); @@ -294,56 +443,69 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_highest_density_sum_test) { ParticleHypothesis::pion()); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Fill grid with track densities auto trackDensityMap = grid.addTrack(params1, mainDensityMap); - auto res1 = grid.getMaxZPosition(mainDensityMap); + auto res1 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res1.ok()); // Maximum should be at z0Trk1 position - BOOST_CHECK_EQUAL(*res1, z0Trk1); + BOOST_CHECK_EQUAL((*res1).first, z0Trk1); // Add second track trackDensityMap = grid.addTrack(params2, mainDensityMap); - auto res2 = grid.getMaxZPosition(mainDensityMap); + auto res2 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res2.ok()); // Trk 2 is closer to z-axis and should yield higher density values // New maximum is therefore at z0Trk2 - BOOST_CHECK_EQUAL(*res2, z0Trk2); + BOOST_CHECK_EQUAL((*res2).first, z0Trk2); // Add small density values around the maximum of track 1 const float densityToAdd = 0.5; - mainDensityMap.at(4) += densityToAdd; - mainDensityMap.at(6) += densityToAdd; + mainDensityMap.at(std::make_pair(4, 0)) += densityToAdd; + mainDensityMap.at(std::make_pair(6, 0)) += densityToAdd; - auto res3 = grid.getMaxZPosition(mainDensityMap); + auto res3 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res3.ok()); // Trk 2 still has the highest peak density value, however, the small // added densities for track 1 around its maximum should now lead to // a prediction at z0Trk1 again - BOOST_CHECK_EQUAL(*res3, z0Trk1); + BOOST_CHECK_EQUAL((*res3).first, z0Trk1); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_removing_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(track_removing) { + const int spatialTrkGridSize = 29; + const int temporalTrkGridSize = 29; + + int trkGridSize = spatialTrkGridSize * temporalTrkGridSize; - double binSize = 0.05; // mm + // bin extent in z and t direction + double binExtent = 0.05; - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + // 1D grid + AdaptiveGridTrackDensity::Config cfg1D(binExtent); + AdaptiveGridTrackDensity grid1D(cfg1D); + + // 2D grid + AdaptiveGridTrackDensity::Config + cfg2D(binExtent, binExtent); + AdaptiveGridTrackDensity grid2D( + cfg2D); // Create some test tracks - Covariance covMat(Covariance::Identity()); + Covariance covMat = makeRandomCovariance(); // Define z0 values for test tracks float z0Trk1 = -0.45; + float t0Trk1 = -0.15; float z0Trk2 = -0.25; + float t0Trk2 = t0Trk1; BoundVector paramVec0; - paramVec0 << 0.1, z0Trk1, 0, 0, 0, 0; + paramVec0 << 0.1, z0Trk1, 0, 0, 0, t0Trk1; BoundVector paramVec1; - paramVec1 << 0.1, z0Trk2, 0, 0, 0, 0; + paramVec1 << 0.1, z0Trk2, 0, 0, 0, t0Trk2; // Create perigee surface std::shared_ptr perigeeSurface = @@ -354,84 +516,119 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_removing_test) { BoundTrackParameters params1(perigeeSurface, paramVec1, covMat, ParticleHypothesis::pion()); - // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; - - // Add track 0 - auto trackDensityMap0 = grid.addTrack(params0, mainDensityMap); - BOOST_CHECK(not mainDensityMap.empty()); - // Grid size should match trkGridSize - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - - // Calculate total density - float densitySum0 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum0 += it->second; - } + // Empty maps + AdaptiveGridTrackDensity::DensityMap mainDensityMap1D; + AdaptiveGridTrackDensity::DensityMap + mainDensityMap2D; + + // Lambda for calculating total density + auto densitySum = [](const auto& densityMap) { + float sum = 0.; + for (auto it = densityMap.begin(); it != densityMap.end(); it++) { + sum += it->second; + } + return sum; + }; - // Add track 0 again - trackDensityMap0 = grid.addTrack(params0, mainDensityMap); - BOOST_CHECK(not mainDensityMap.empty()); + // Add track 0 to 1D grid + auto firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); + BOOST_CHECK(not mainDensityMap1D.empty()); + // Grid size should match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + float firstDensitySum1D = densitySum(mainDensityMap1D); + + // Add track 0 to 2D grid + auto firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); + BOOST_CHECK(not mainDensityMap2D.empty()); + // Grid size should match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + float firstDensitySum2D = densitySum(mainDensityMap2D); + + // Add track 0 again to 1D grid + firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); + BOOST_CHECK(not mainDensityMap1D.empty()); + // Grid size should still match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + // Calculate new total density ... + float secondDensitySum1D = densitySum(mainDensityMap1D); + // ... and check that it's twice as large as before + BOOST_CHECK(2 * firstDensitySum1D == secondDensitySum1D); + + // Add track 0 again to 2D grid + firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); + BOOST_CHECK(not mainDensityMap2D.empty()); // Grid size should still match trkGridSize - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + // Calculate new total density ... + float secondDensitySum2D = densitySum(mainDensityMap2D); + // ... and check that it's twice as large as before + BOOST_CHECK(2 * firstDensitySum2D == secondDensitySum2D); + + // Remove track 0 from 1D grid + grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); // Calculate new total density - float densitySum1 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum1 += it->second; - } - - BOOST_CHECK(2 * densitySum0 == densitySum1); - - // Remove track 0 - grid.subtractTrack(trackDensityMap0, mainDensityMap); + float thirdDensitySum1D = densitySum(mainDensityMap1D); + // Density should be old one again + BOOST_CHECK(firstDensitySum1D == thirdDensitySum1D); + // Grid size should still match spatialTrkGridSize (removal does not change + // grid size) + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + // Remove track 0 from 2D grid + grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); // Calculate new total density - float densitySum2 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum2 += it->second; - } - + float thirdDensitySum2D = densitySum(mainDensityMap2D); // Density should be old one again - BOOST_CHECK(densitySum0 == densitySum2); - // Grid size should still match trkGridSize (removal does not touch grid size) - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - - // Add track 1, overlapping track 0 - auto trackDensityMap1 = grid.addTrack(params1, mainDensityMap); - - int nNonOverlappingBins = int(std::abs(z0Trk1 - z0Trk2) / binSize + 1); - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize + nNonOverlappingBins); - - float densitySum3 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum3 += it->second; - } - - // Remove second track 0 - grid.subtractTrack(trackDensityMap0, mainDensityMap); - - float densitySum4 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum4 += it->second; - } - + BOOST_CHECK(firstDensitySum2D == thirdDensitySum2D); + // Grid size should still match trkGridSize (removal does not change grid + // size) + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + + // Add track 1 to 1D grid (overlaps with track 0!) + auto secondTrackDensityMap1D = grid1D.addTrack(params1, mainDensityMap1D); + int nNonOverlappingBins1D = int(std::abs(z0Trk1 - z0Trk2) / binExtent + 1); + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), + spatialTrkGridSize + nNonOverlappingBins1D); + float fourthDensitySum1D = densitySum(mainDensityMap1D); + + // Add track 1 to 2D grid (overlaps with track 0!) + auto secondTrackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); + int nNonOverlappingBins2D = nNonOverlappingBins1D * temporalTrkGridSize; + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), + trkGridSize + nNonOverlappingBins2D); + float fourthDensitySum2D = densitySum(mainDensityMap2D); + + // Remove second track 0 from 1D grid + grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); + float fifthDensitySum1D = densitySum(mainDensityMap1D); // Density should match differences of removed tracks - CHECK_CLOSE_ABS(densitySum4, densitySum3 - densitySum0, 1e-5); + CHECK_CLOSE_REL(fifthDensitySum1D, fourthDensitySum1D - firstDensitySum1D, + 1e-5); - // Remove track 1 - grid.subtractTrack(trackDensityMap1, mainDensityMap); + // Remove second track 0 from 2D grid + grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); + float fifthDensitySum2D = densitySum(mainDensityMap2D); + // Density should match differences of removed tracks + CHECK_CLOSE_REL(fifthDensitySum2D, fourthDensitySum2D - firstDensitySum2D, + 1e-5); + // Remove track 1 from 1D grid + grid1D.subtractTrack(secondTrackDensityMap1D, mainDensityMap1D); // Size should not have changed - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize + nNonOverlappingBins); - - float densitySum5 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum5 += it->second; - } - - // Grid is now empty after all tracks were removed - CHECK_CLOSE_ABS(densitySum5, 0., 1e-5); + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), + spatialTrkGridSize + nNonOverlappingBins1D); + float sixthDensitySum1D = densitySum(mainDensityMap1D); + // 1D grid is now empty after all tracks were removed + CHECK_CLOSE_ABS(sixthDensitySum1D, 0., 1e-4); + + // Remove track 1 from 2D grid + grid2D.subtractTrack(secondTrackDensityMap2D, mainDensityMap2D); + // Size should not have changed + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), + trkGridSize + nNonOverlappingBins2D); + float sixthDensitySum2D = densitySum(mainDensityMap2D); + // 2D grid is now empty after all tracks were removed + CHECK_CLOSE_ABS(sixthDensitySum2D, 0., 1e-4); } } // namespace Test From 020a076bebfea8f9c01a4fb0813f5bcdb495590c Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 6 Oct 2023 18:41:46 +0200 Subject: [PATCH 6/7] ci: Drop GitHub Actions LCG builds (#2516) We keep the LCG jobs on the GitLab side --- .github/workflows/builds.yml | 109 ----------------------------------- 1 file changed, 109 deletions(-) diff --git a/.github/workflows/builds.yml b/.github/workflows/builds.yml index 32072613347..b42793a85b8 100644 --- a/.github/workflows/builds.yml +++ b/.github/workflows/builds.yml @@ -22,115 +22,6 @@ env: CCACHE_KEY_SUFFIX: r1 jobs: - lcg: - runs-on: ubuntu-latest - container: ghcr.io/acts-project/${{ matrix.image }}:v41 - strategy: - matrix: - image: - - centos7-lcg100-gcc10 - - centos7-lcg101-gcc11 - - centos8-lcg100-gcc10 - - centos8-lcg101-gcc11 - env: - SETUP: source /opt/lcg_view/setup.sh - INSTALL_DIR: ${{ github.workspace }}/install - ACTS_LOG_FAILURE_THRESHOLD: WARNING - steps: - - uses: actions/checkout@v3 - - - name: Restore ccache - uses: actions/cache/restore@v3 - id: ccache-restore - with: - path: ${{ github.workspace }}/ccache - key: ${{ runner.os }}-ccache-${{ matrix.image }}_${{ env.CCACHE_KEY_SUFFIX }}_${{ github.sha }} - restore-keys: | - ${{ runner.os }}-ccache-${{ matrix.image }}_${{ env.CCACHE_KEY_SUFFIX }}_ - - - name: Configure - # setting CMAKE_CXX_STANDARD=17 is a workaround for a bug in the - # dd4hep CMake configuration that gets triggered on recent CMake - # versions - run: > - ${SETUP} && - ccache -z && - cmake -B build -S . - -GNinja - -DCMAKE_CXX_COMPILER_LAUNCHER=ccache - -DCMAKE_BUILD_TYPE=Release - -DCMAKE_CXX_FLAGS=-Werror - -DCMAKE_CXX_STANDARD=17 - -DCMAKE_INSTALL_PREFIX="${INSTALL_DIR}" - -DACTS_LOG_FAILURE_THRESHOLD=WARNING - -DACTS_BUILD_EXAMPLES_PYTHON_BINDINGS=ON - -DACTS_FORCE_ASSERTIONS=ON - -DACTS_BUILD_EXAMPLES=ON - -DACTS_BUILD_PLUGIN_DD4HEP=OFF - -DACTS_BUILD_PLUGIN_TGEO=ON - -DACTS_BUILD_PLUGIN_IDENTIFICATION=ON - -DACTS_BUILD_PLUGIN_JSON=ON - -DACTS_BUILD_FATRAS=ON - -DACTS_BUILD_PLUGIN_LEGACY=ON - -DACTS_BUILD_PLUGIN_AUTODIFF=OFF - -DACTS_BUILD_BENCHMARKS=ON - -DACTS_BUILD_UNITTESTS=ON - -DACTS_BUILD_INTEGRATIONTESTS=ON - -DACTS_BUILD_EXAMPLES_DD4HEP=OFF - -DACTS_BUILD_PLUGIN_EDM4HEP=OFF - -DACTS_BUILD_EXAMPLES_GEANT4=ON - -DACTS_BUILD_EXAMPLES_HEPMC3=ON - -DACTS_BUILD_EXAMPLES_PYTHIA8=ON - -DACTS_BUILD_FATRAS_GEANT4=ON - -DACTS_BUILD_FATRAS=ON - -DACTS_BUILD_ALIGNMENT=ON - -DACTS_BUILD_ANALYSIS_APPS=ON - -DACTS_BUILD_PLUGIN_ACTSVG=ON - - - name: Build - run: ${SETUP} && cmake --build build - - - name: ccache stats - run: ${SETUP} && ccache -s - - - name: Save ccache - uses: actions/cache/save@v3 - if: always() - with: - path: ${{ github.workspace }}/ccache - key: ${{ steps.ccache-restore.outputs.cache-primary-key }} - - - name: Unit tests - run: ${SETUP} && cmake --build build --target test - - - name: Integration tests - run: ${SETUP} && cmake --build build --target integrationtests - - - name: Install - run: ${SETUP} && cmake --build build --target install - - - uses: actions/upload-artifact@v3 - with: - name: acts-${{ matrix.image }} - path: ${{ env.INSTALL_DIR }} - - - name: Downstream configure - run: > - ${SETUP} && - cmake -B build-downstream -S Tests/DownstreamProject - -GNinja - -DCMAKE_BUILD_TYPE=Release - -DCMAKE_CXX_FLAGS=-Werror - -DCMAKE_CXX_STANDARD=17 - -DCMAKE_PREFIX_PATH="${INSTALL_DIR}" - -DDD4HEP=OFF - - - name: Downstream build - run: ${SETUP} && cmake --build build-downstream - - - name: Downstream run - run: ${SETUP} && ./build-downstream/bin/ShowActsVersion - linux_ubuntu: runs-on: ubuntu-latest container: ghcr.io/acts-project/ubuntu2204:v41 From d8c447c139eab534a30fef6003570c5003134fdf Mon Sep 17 00:00:00 2001 From: felix-russo <72298366+felix-russo@users.noreply.github.com> Date: Sat, 7 Oct 2023 11:54:53 +0200 Subject: [PATCH 7/7] perf: use analytical solution for straight tracks during impact point estimation (#2506) We use the Newton method to find the 3D PCA for helical tracks. Straight tracks are currently approximated by setting the helix radius to a very large value. In this PR, the approximation is replaced by an analytical solution (which exists for straight tracks). This should lead to a better performance and (slightly) more precise results. To check the analytical solution, a unit test is added. A derivation of the analytical solution can be found in the updated reference: [Track_Linearization_and_3D_PCA.pdf](https://github.com/acts-project/acts/files/12794276/Track_Linearization_and_3D_PCA.pdf) --- .../Acts/Vertexing/ImpactPointEstimator.hpp | 8 +- .../Acts/Vertexing/ImpactPointEstimator.ipp | 80 ++++++++++++++----- .../Vertexing/ImpactPointEstimatorTests.cpp | 60 +++++++++++--- 3 files changed, 111 insertions(+), 37 deletions(-) diff --git a/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp b/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp index f12a3054e91..30d385859ad 100644 --- a/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp +++ b/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp @@ -40,7 +40,7 @@ struct ImpactParametersAndSigma { /// /// @brief Estimator for impact point calculations /// A description of the underlying mathematics can be found here: -/// https://github.com/acts-project/acts/pull/2414 +/// https://github.com/acts-project/acts/pull/2506 /// TODO: Upload reference at a better place template > @@ -81,10 +81,6 @@ class ImpactPointEstimator { int maxIterations = 20; /// Desired precision of deltaPhi in Newton method double precision = 1.e-10; - /// Minimum q/p value - double minQoP = 1e-15; - /// Maximum curvature value - double maxRho = 1e+15; }; /// @brief Constructor @@ -147,6 +143,8 @@ class ImpactPointEstimator { /// corresponding 3D PCA. Returns also the momentum direction at the 3D PCA. /// The template parameter nDim determines whether we calculate the 3D /// distance (nDim = 3) or the 4D distance (nDim = 4) to the 3D PCA. + /// @note For straight tracks we use an analytical solution; for helical + /// tracks we use the Newton method. /// /// @tparam nDim Number of dimensions used to compute compatibility /// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we diff --git a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp index 6e023f384d9..a5a698e0295 100644 --- a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp +++ b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp @@ -241,13 +241,65 @@ Acts::ImpactPointEstimator:: // Reference point R Vector3 refPoint = trkParams.referenceSurface().center(gctx); - // Perigee parameters (parameters of 2D PCA) + // Extract charge-related particle parameters + double absoluteCharge = trkParams.particleHypothesis().absoluteCharge(); + double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP]; + + // Z-component of the B field at the reference position. + // Note that we assume a constant B field here! + auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache); + if (!fieldRes.ok()) { + return fieldRes.error(); + } + double bZ = (*fieldRes)[eZ]; + + // The particle moves on a straight trajectory if its charge is 0 or if there + // is no B field. In that case, the 3D PCA can be calculated analytically, see + // Sec 3.2 of the reference. + if (absoluteCharge == 0. or bZ == 0.) { + // Momentum direction (constant for straight tracks) + Vector3 momDirStraightTrack = trkParams.direction(); + + // Current position on the track + Vector3 positionOnTrack = trkParams.position(gctx); + + // Distance between positionOnTrack and the 3D PCA + ActsScalar distanceToPca = + (vtxPos.template head<3>() - positionOnTrack).dot(momDirStraightTrack); + + // 3D PCA + ActsVector pcaStraightTrack; + pcaStraightTrack.template head<3>() = + positionOnTrack + distanceToPca * momDirStraightTrack; + if constexpr (nDim == 4) { + // Track time at positionOnTrack + double timeOnTrack = trkParams.parameters()[BoundIndices::eBoundTime]; + + ActsScalar m0 = trkParams.particleHypothesis().mass(); + ActsScalar p = std::abs(absoluteCharge / qOvP); + + // Speed in units of c + ActsScalar beta = p / std::hypot(p, m0); + + pcaStraightTrack[3] = timeOnTrack + distanceToPca / beta; + } + + // Vector pointing from the vertex position to the 3D PCA + ActsVector deltaRStraightTrack = pcaStraightTrack - vtxPos; + + return std::make_pair(deltaRStraightTrack, momDirStraightTrack); + } + + // Charged particles in a constant B field follow a helical trajectory. In + // that case, we calculate the 3D PCA using the Newton method, see Sec 4.2 in + // the reference. + + // Spatial Perigee parameters (i.e., spatial parameters of 2D PCA) double d0 = trkParams.parameters()[BoundIndices::eBoundLoc0]; double z0 = trkParams.parameters()[BoundIndices::eBoundLoc1]; + // Momentum angles at 2D PCA double phiP = trkParams.parameters()[BoundIndices::eBoundPhi]; double theta = trkParams.parameters()[BoundIndices::eBoundTheta]; - double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP]; - double tP = trkParams.parameters()[BoundIndices::eBoundTime]; // Functions of the polar angle theta for later use double sinTheta = std::sin(theta); double cotTheta = 1. / std::tan(theta); @@ -256,22 +308,8 @@ Acts::ImpactPointEstimator:: // Note that phi corresponds to phiV in the reference. double phi = phiP; - // B-field z-component at the reference position. - // Note that we assume a constant B field here! - auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache); - if (!fieldRes.ok()) { - return fieldRes.error(); - } - double bZ = (*fieldRes)[eZ]; - // Signed radius of the helix on which the particle moves - double rho = 0; - // Curvature is infinite w/o b field - if (bZ == 0. || std::abs(qOvP) < m_cfg.minQoP) { - rho = m_cfg.maxRho; - } else { - rho = sinTheta * (1. / qOvP) / bZ; - } + double rho = sinTheta * (1. / qOvP) / bZ; // Position of the helix center. // We can set the z-position to a convenient value since it is not fixed by @@ -308,9 +346,11 @@ Acts::ImpactPointEstimator:: helixCenter + rho * Vector3(-sinPhi, cosPhi, -cotTheta * phi); if constexpr (nDim == 4) { + // Time at the 2D PCA P + double tP = trkParams.parameters()[BoundIndices::eBoundTime]; + ActsScalar m0 = trkParams.particleHypothesis().mass(); - ActsScalar p = - std::abs(trkParams.particleHypothesis().absoluteCharge() / qOvP); + ActsScalar p = std::abs(absoluteCharge / qOvP); // Speed in units of c ActsScalar beta = p / std::hypot(p, m0); diff --git a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp index d812db2f0f6..ada4d758d16 100644 --- a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp @@ -25,6 +25,7 @@ #include "Acts/MagneticField/NullBField.hpp" #include "Acts/Propagator/EigenStepper.hpp" #include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/StraightLineStepper.hpp" #include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -54,10 +55,13 @@ using namespace Acts::UnitLiterals; using Acts::VectorHelpers::makeVector4; using MagneticField = Acts::ConstantBField; +using StraightPropagator = Acts::Propagator; using Stepper = Acts::EigenStepper<>; using Propagator = Acts::Propagator; using Estimator = Acts::ImpactPointEstimator; +using StraightLineEstimator = + Acts::ImpactPointEstimator; const Acts::GeometryContext geoContext; const Acts::MagneticFieldContext magFieldContext; @@ -95,7 +99,7 @@ Estimator makeEstimator(double bZ) { Estimator::Config cfg(field, std::make_shared( std::move(stepper), detail::VoidNavigator(), - getDefaultLogger("Prop", Logging::Level::VERBOSE))); + getDefaultLogger("Prop", Logging::Level::WARNING))); return Estimator(cfg); } @@ -132,7 +136,7 @@ BOOST_AUTO_TEST_SUITE(VertexingImpactPointEstimator) // Check `calculateDistance`, `estimate3DImpactParameters`, and // `getVertexCompatibility`. -BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0, +BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3D, tracks, d0, l0, t0, phi, theta, p, q) { auto particleHypothesis = ParticleHypothesis::pion(); @@ -193,12 +197,24 @@ BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0, BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, q, vx0, vy0, vz0, vt0) { using Propagator = Acts::Propagator; + using StraightPropagator = Acts::Propagator; + + // Set up quantities for constant B field auto field = std::make_shared(Vector3(0, 0, 2_T)); Stepper stepper(field); auto propagator = std::make_shared(std::move(stepper)); Estimator::Config cfg(field, propagator); Estimator ipEstimator(cfg); - Estimator::State state(magFieldCache()); + Estimator::State ipState(magFieldCache()); + + // Set up quantities for B = 0 + auto zeroField = std::make_shared(Vector3(0, 0, 0)); + StraightLineStepper straightLineStepper; + auto straightLinePropagator = + std::make_shared(straightLineStepper); + StraightLineEstimator::Config zeroFieldCfg(zeroField, straightLinePropagator); + StraightLineEstimator zeroFieldIPEstimator(zeroFieldCfg); + StraightLineEstimator::State zeroFieldIPState(magFieldCache()); // Vertex position and vertex object Vector4 vtxPos(vx0, vy0, vz0, vt0); @@ -240,7 +256,7 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, // Perigee surface at vertex position auto refPerigeeSurface = Surface::makeShared(refPoint); - // Propagate to the 2D PCA of the reference point + // Set up the propagator options (they are the same with and without B field) PropagatorOptions pOptions(geoContext, magFieldContext); auto intersection = refPerigeeSurface ->intersect(geoContext, params.position(geoContext), @@ -248,13 +264,21 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, .closest(); pOptions.direction = Direction::fromScalarZeroAsPositive(intersection.pathLength()); + + // Propagate to the 2D PCA of the reference point in a constant B field auto result = propagator->propagate(params, *refPerigeeSurface, pOptions); BOOST_CHECK(result.ok()); - const auto& refParams = *result->endParameters; + // Propagate to the 2D PCA of the reference point when B = 0 + auto zeroFieldResult = + straightLinePropagator->propagate(params, *refPerigeeSurface, pOptions); + BOOST_CHECK(zeroFieldResult.ok()); + const auto& zeroFieldRefParams = *zeroFieldResult->endParameters; + BOOST_TEST_CONTEXT( - "Check time at 2D PCA (i.e., function getImpactParameters)") { + "Check time at 2D PCA (i.e., function getImpactParameters) for helical " + "tracks") { // Calculate impact parameters auto ipParams = ipEstimator .getImpactParameters(refParams, vtx, geoContext, @@ -270,14 +294,14 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, 1e-3); } - BOOST_TEST_CONTEXT( - "Check time at 3D PCA (i.e., function getDistanceAndMomentum)") { + auto checkGetDistanceAndMomentum = [&vtxPos, &corrMomDir, corrTimeDiff]( + const auto& ipe, const auto& rParams, + auto& state) { // Find 4D distance and momentum of the track at the vertex starting from // the perigee representation at the reference position - auto distAndMom = - ipEstimator - .getDistanceAndMomentum<4>(geoContext, refParams, vtxPos, state) - .value(); + auto distAndMom = ipe.template getDistanceAndMomentum<4>( + geoContext, rParams, vtxPos, state) + .value(); ActsVector<4> distVec = distAndMom.first; Vector3 momDir = distAndMom.second; @@ -293,6 +317,18 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, // Momentum direction should correspond to the momentum direction at the // vertex CHECK_CLOSE_OR_SMALL(momDir, corrMomDir, 1e-5, 1e-4); + }; + + BOOST_TEST_CONTEXT( + "Check time at 3D PCA (i.e., function getDistanceAndMomentum) for " + "straight tracks") { + checkGetDistanceAndMomentum(zeroFieldIPEstimator, zeroFieldRefParams, + zeroFieldIPState); + } + BOOST_TEST_CONTEXT( + "Check time at 3D PCA (i.e., function getDistanceAndMomentum) for " + "helical tracks") { + checkGetDistanceAndMomentum(ipEstimator, refParams, ipState); } }