% Tento soubor nahraďte vlastním souborem s obsahem práce.
%=========================================================================
% Autoři: Michal Bidlo, Bohuslav Křena, Jaroslav Dytrych, Petr Veigend a Adam Herout 2019
\chapter{Úvod}

Teorie velkého třesku, nejpopulárnější kosmologická teorie o vzniku a prvních okamžicích existence našeho vesmíru tvrdí, že vesmír, jak ho známe vznikl mohutnou \uv{explozí} a po miliardy let postupně chladl a rozpínal se do podoby, kterou známe dnes. Přestože současné teorie odhadují čas od velkého třesku na přibližně 13,8 miliardy let, o tom, jak vypadal náš vesmír v prvních okamžicích od svého vzniku víme jen velmi málo a přesný průběh velkého třesku nám není znám. Předpokládáme, že velký třesk probíhal ve fázích, kdy se jeho podoba drasticky lišila od podoby, jakou má dnes. Mezi prvními z těchto fází (řádově mikrosekundy od velkého třesku) byla pravděpodobně fáze, kdy veškerá hmota v mladém vesmíru se vyskytovala ve formě takzvaného \textbf{kvark-gluonového plazmatu}\footnote{z anglického QGP -- Quark-Gluon Plasma}, kdy veškerá hmota byla natolik hustá a horká, že se jednotlivé částice hmoty mohly volně a samostatně pohybovat (jev, který například u kvarků za normálních okolností nelze sledovat), aniž by je jakákoliv síla dokázala spoutat. Narozdíl od jiných forem plazmatu se tedy volně pohybovaly všechny elementární částice, aniž by se sdružovaly do větších celků, jako je tomu například u běžně sledových forem plazmatu, kdy sledujeme převážně volný pohyb elektronů a atomových jader, odtržených od jejich elektronových obalů. Kvark-gluonové plazma, jak již název napovídá, posouvá tento jev o dva kroky dál, kdy pozorujeme nejen rozpad jader na jednotlivé hadrony (např. neutrony, protony nebo mezony), ale i ty se dále rozpadají na jednotlivé kvarky a gluony, které se v QGP volně pohybují. Předpokládáme, že stav hmoty s podobnými vlastnostmi dnes můžeme ve vesmíru sledovat pouze za extrémních podmínek, například v jádrech neutronových hvězd.

Kvantová chromodynamika\footnote{z anglického QCD -- Quantum ChromoDynamics} je disciplína teoretické fyziky, zabývající se silnou interakcí mezi \textbf{barevně nabitými částicemi}, tedy interakcí mezi kvarky, zprostředkovávanou částicí zvanou \textit{gluon}, podobně jako je elektromagnetická interakce zprostředkovávána \textit{fotony}. Barevně nabité částice jsou specifické tím, že je za \uv{normálních} teplot a tlaků nelze pozorovat samostatně: vždy se vyskytují po trojicích (hadrony, například \textit{proton} a \textit{neutron}) nebo po dvojicích (mezony, například \textit{pion}), zřídkakdy také po čtveřicích, pěticích a šesticích (tetrakvarky, pentakvarky, hexakvarky). To je způsobeno jevem zvaným \textit{barevné uvěznění}\footnote{z anglického \textit{colour confinement}}, který mimo jiné říká, že dodáme-li kvarkům dostatek energie na to, aby se dokázali rozdělit, je energeticky výhodnější vytvořit nový pár kvark-antikvark. V případě, že stlačíme hmotu do tak extrémních tlaků, že jednotlivé \textit{nukleony} nemají dost prostoru k tomu, aby mohly samostatně existovat, přestávají být kvarky drženy pohromadě gluony a mohou se volně pohybovat. Takové fázi hmoty říkáme \textit{kvark-gluonové plazma}.

Jeden ze způsobů, jak můžeme laboratorně pozorovat a studovat strukturu hmoty za extrémních podmínek je pomocí řízených kolizí těžkých iontů s vysokou energií, jako jsou například jádra zlata či olova. Taková jádra urychlujeme na rychlosti velmi blízké rychlosti světla a necháváme s obrovskou energií kolidovat s terči detektorů (v případě experimentu CBM vyrobený z tenkého proužku zlata), čímž jádra přivedeme k přechodnému stavu s extrémně vysokou hustotou a teplotou, který fyzikové nazývají \textit{fireball} (ohnivá koule). Aby bylo možné vytvořit vhodné podmínky pro zkoumání hmoty za tak extrémních podmínek, musíme použít urychlovače s velmi vysokou energií paprsku. Experiment CBM, připravovaný experiment při moderním výzkumném středisku FAIR v německém Darmstadtu, bude hrát důležitou roli při zkoumání fázového diagramu kvantové chromodynamiky. Schopnost sledovat takto extrémní stav hmoty je ale vykoupena obrovským množstvím experimentálních dat (cca. 1 TB/s), které je potřeba zpracovávat \textit{on-line}.

Rychlé zpracování velkého objemu hrubých experimentálních dat je netriviální problém, jehož vhodné řešení je stěžejní pro úspěšnost fyzikálního experimentu. Ve svých pracích se tímto problémem zabývali například prof. Dr. Ivan Kisel a Dr. Valentina Akishina, kteří se zabývali vývojem algoritmu rekonstrukce trajektorií částic pro potřeby detektoru CBM a jeho možnostmi akcelerace, v tomto případě pro platformu \textit{Intel Xeon Phi}. Tato práce si klade za cíl prozkoumat alternativní způsoby akcelerace rekonstrukce trajektorií pro experiment CBM, zejména pak za pomocí protokolu MPI pro paraleliyzaci na mnohaprocesorových farmách. Kapitola \ref{CBM} poskytuje stručný úvod do výzkumného střediska FAIR a experimentu CBM, kapitola \ref{track_reconstruction} se poté dopodrobna zabývá jednotlivými algoritmy, použitými pro rekonstrukci trajektorií a aproximaci vlastností jejich částic a jejich zhodnocením proveditelnosti jejich akcelerace na hardwarových platformách. Kapitola \ref{implementation} se zabývá vývojem vhodného simulátoru pohybu částic uvnitř detektoru CBM a jeho použitím v rámci vlastního algoritmu sledování částic. Nakonec Kapitola \ref{reconstruction} popisuje vlastní implementaci rekonstrukčního algoritmu, jeho akceleraci pomocí metod distribuovaného počítání a diskutuje další možné způsoby akcelerace výpočtu.

\chapter{Experiment CBM a výzkumné středisko FAIR}
\label{CBM}

Tato kapitola je věnována stručnému popisu pracoviště FAIR\footnote{FAIR -- Facility for Antiproton and Ion Research (středisko pro výzkum antiprotonů a iontů)} a jeho výpočetní infrastruktury, fyzikálního experimentu CBM\footnote{CBM -- Compressed Baryonic Matter (stlačená baryonová hmota)} a jeho cílů. Kapitola se velmi stručně zabývá detektorem STS\footnote{STS -- Silicon Tracking System (křemíkový sledovací systém)} a supravodivým dipólovým magnetem, ve kterém se STS nachází. Zvláštní pozornost je pak věnována výpočetní farmě při GSI \textbf{Green IT Cube}. 

\section{Pracoviště FAIR}
Facility for Antiproton and Ion Research (FAIR) je mezinárodní výzkumné centrum, v současné době ve výstavbě v německém Darmstadtu při Helmholtzově centru pro výzkum těžkých iontů GSI\footnote{GSI -- Gesellschaft für Schwerionenforschung (sdružení pro výzkum těžkých iontů)}. Projekt je realizován partnery z Finska, Francie, Indie, Německa, Polska, Rumunska, Ruska, Slovinska, Švédska a Velké Británie. Česká republika se stala aspirujícím partnerem v roce 2018. \cite{gsi_about, fair_overview} FAIR se v budoucnosti stane jedním z největších a nejsložitějších urychlovačů na světě. Výsadou urychlovače bude jeho schopnost produkovat svazky částic všech prvků periodické tabulky (nebo jejich iontů), případně pak svazky antiprotonů, urychlené až na 99\% rychlosti světla. \cite{fair_accelerator}

\begin{figure*}[ht]\centering
  \centering
  \includegraphics[width=0.43\linewidth]{obrazky-figures/FAIR_diagram.jpg}\hfill
  \includegraphics[width=0.53\linewidth]{obrazky-figures/FAIR_construction2021.jpg}
  \caption{Rozvržení objektu FAIR (vlevo). Stávající urychlovače GSI UNILAC a SIS18 (vyobrazené modře) budou kromě obsluhy stávajícího komplexu GSI (vyobrazeno šedě) sloužit jako injektor pro dva velké synchrotrony SIS100 a SIS300 (červeně). Částicové svazky budou rozváděny do různých experimentálních stanic nebo k terčovým stanovištím pro experimenty vyžadující svazky antiprotonů nebo radioaktivních izotopů. Převzato z \cite{AUMANN20073}. Fotografie areálu GSI a staveniště FAIR (vpravo), převzatá z \cite{fair_construction_site}, pořízená v dubnu 2021.}
  \label{fig:fair_layout}
\end{figure*}

\subsection*{Stávající infrastruktura GSI}
Stávající zařízení GSI se v budoucnu stanou součástí infrastruktury FAIR a jeho dva urychlovače budou sloužit jako první urychlovací stupeň pro připravované experimenty. Lineární urychlovač UNILAC\footnote{UNILAC -- UNIversal Linear ACcelerator (univerzální lineární urychlovač)}, dlouhý 120 metrů, urychluje částice až na 20\% rychlosti světla. Odsud částice putují k prvním experimentům a prvnímu prstencovému urychlovači SIS18\footnote{SIS -- Schwerionensynchrotron (synchrotron pro těžké ionty)} o obvodu 120 metrů, který částice dále urychluje až na 90\% rychlosti světla a dále vysílá ke stávajícím experimentům GSI a úložným prstencům. \cite{fair_accelerator}

\subsection*{Prstencové urychlovače SIS100 a SIS300}
FAIR bude pracovat s částicovými svazky s vysokými intenzitami iontů ($10^9$ částic za sekundu v případě iontů zlata) nebo pozitronů ($10^{13}$ částic za sekundu). \cite{akishina4D}
V současné době jsou ve výstavbě nové prstencové urychlovače SIS100 a SIS300, uloženy v podzemních tunelech v hloubce až 17 metrů s obvodem 1100 metrů. Částice budou do urychlovačů SIS100 a SIS300 vháněny z urychlovače SIS18 a dále urychleny až na 99\% rychlosti světla. Takto urychlené částice jsou pak buďto použity přímo k experimentům, k přípravě tzv. sekundárních částic, nebo vehnány do tzv. úložných prstenců k pozdějšímu použití. \cite{fair_accelerator} 

\subsection*{Úložné prstence}
K prstencovým urychlovačům je připojen komplexní systém úložných prstenců a dalších experimentálních stanic. Úložné prstence vytvářejí pomocí silného magnetického pole cyklickou trajektorii příchozích svazků, na rozdíl od prstencových urychlovačů ale svazkům již nedodávají žádnou energii, tedy je již dále neurychlují. Tento přístup je výhodný, protože příprava vysokoenergetických svazků je náročná a bez použití úložných prstenců jsou připravené svazky po průletu experimentálním stanovištěm ztraceny. Použitím úložných prstenců zajistíme, že vědecké týmy mohou na stejných částicích pracovat při každém průletu experimentálním stanovištěm a dále manipulovat s vlastnostmi svazku. Navíc použití úložných prstenců de facto dále zvyšuje intenzitu svazku bez nutnosti dalšího použití urychlovače.

Pro skladování a manipulaci se svazky se počítá s několika prstenci, včetně vysokoenergetického prstence HESR\footnote{HESR -- High Energy Storage Ring (vysokoenergetický úložný prstenec)}, určeného k práci s antiprotonovými svazky, sběrného prstence CR\footnote{CR -- Collector Ring (sběrný prstenec)} a experimentálního prstence NESR\footnote{NESR -- New Experimental Storage Ring (nový experimentální úložný prstenec)}. Komplex úložných prstenců umožňuje chlazení svazků a zahrnuje terčová stanoviště pro rozptylové experimenty. \cite{AUMANN20073}

\subsection*{Green IT Cube}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/green_cube_photo.jpg}
    \caption{Fotografie datového centra Green IT Cube, Administrativní budova GSI je viditelná vpravo. Obrázek byl převzat z webových stránek GSI \cite{gsi_green_cube_covid}.}
    \label{fig:green_it_cube}
\end{figure}

V roce 2008 se odhadovalo, že informační a komunikační technologie byly zodpovědné za 2\% světové produkce skleníkových plynů, což je množství srovnatelné s emisemi letadlové dopravy. \cite{green_cube_powerpoint}. Firma Intel v roce 2012 odhadovala rozložení emisí své ICT\footnote{ICT -- Information and Communication Technologies (informační a komunikační technologie)} infrastruktury následovně: \cite{intel_sustainability}

\begin{itemize}
    \item 6\% připadá na uživatelskou výpočetní techniku (kancelářské počítače, monitory, tiskárny apod.),
    \item 24\% připadá na ICT infrastrukturu mimo datová centra,
    \item 70\% připadá na datová centra.
\end{itemize}

Podle výpočtů institutu Borderstep činila již v roce 2008 spotřeba elektrické energie datových center v Německu 10,1 TWh, což odpovídá 1,8\% celkové německé spotřeby toho roku a celkovým nákladům přibližně 1,1 miliardy eur. Přitom výpočetní výkon (a s ním i nároky na výkon chlazení) datových center se neustále zvyšuje: Studie stejného institutu a autora o deset let později uvádí, že v roce 2017 činila spotřeba energie datových center v Německu 13,2 Terawatthodin.

Podle EYP Mission Critical Facilities Inc., New York \cite{green_cube_powerpoint}, typicky ne víc než 50\% energie potřebné k provozu datových center je využito na vlastní napájení IT vybavení. Uveďme si metriku pro určení energetické efektivity Datového centra: \textbf{PUE}\footnote{PUE -- Power Usage Effectiveness (efektivita spotřeby energie)} je metrika, která zkoumá, jak velkou část celkové energetické spotřeby datového centra vyžaduje napájení výpočetní techniky:

\begin{equation}\label{eq:pue}
    \textbf{PUE} = \frac{\text{Celková spotřeba objektu}}{\text{Spotřeba výpočetní techniky}}
\end{equation}

Přičemž PUE typického datového centra se může pohybovat kdekoli mezi $1,6$ a $3,5$, a teoretické $\text{PUE} = 1$ by mělo takové datové centrum, které využije veškerou svoji spotřebu energie pouze k napájení výpočetní techniky.

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{obrazky-figures/green_cube_cooling.png}\\
    \includegraphics[width=0.6\textwidth]{obrazky-figures/green_cube_plan.png}
    \caption{Diagram chladícího systému datového centra Green IT Cube (nahoře), půdorys přízemního patra budovy (dole). Obrázky byly přejaty z \cite{green_cube_powerpoint}.}
    \label{fig:green_cube_plans}
\end{figure}

Součástí pracoviště GSI je nové datové centrum \textbf{Green IT Cube}, které si mimo jiné klade za cíl prozkoumat způsoby, jak snížit zátěž datových center na životní prostředí. Green IT Cube je postaveno tak, aby mohlo být osazeno 768 standardními devatenáctipalcovími rozvaděči v šesti patrech budovy. Zařízení disponuje kapacitou až 15 MW plně redundantního chladícího výkonu s PUE menším než 1,07. Takové efektivity bylo dosaženo pomocí speciálního dvouokruhového vodního chlazení, kdy sekundární okruh uvnitř datového centra chladí jednotlivé rozvaděče (nikoliv jednotlivá zařízení v nich, díky čemuž je možné rozvaděče osadit výrazně širším výběrem výpočetní techniky), výměník tepla v přízemí budovy pak odpadní teplo převádí do primárního okruhu, které je pak částečně využito k vytápění kanceláří a jídelny budovy GSI (viz obrázek \ref{fig:green_cube_plans}). Po uvedení do plného výkonu bude mít pracoviště FAIR k dispozici výpočetní výkon tří set tisíc procesorových jader, až 135 Petabajtů (35 PB pro stávající infrastrukturu GSI a 100 PB v rámci rozšíření pro FAIR) online úložiště na pevných discích a robota s 4280ti sloty pro magnetické pásky pro archivaci až 62 PB nekomprimovaných experimentálních dat. Pro experimenty FAIR se dále počítá s vysokou přenosovou rychlostí přes 1 Terabajt za sekundu pro přenos velkého objemu nezpracovaných experimentálních dat z pracovišť jednotlivých experimentů\cite{gsi_green_cube_overview}.

Při uvedení do provozu v roce 2014 obdrželo centrum Green IT Cube mnohá ocenění, včetně prvního místa v žebříčku superpočítačů Green500, které řadí nejefektivnější superpočítače podle GFLOPS/W.

\section{Experiment CBM}

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{obrazky-figures/cbm_variants.png}
    \caption{Rozvržení experimentu CBM v jeho dvou operačních nastaveních: nastavení pro detekci elektronů (obrázek nahoře) a pro detekci mionů (obrázek dole). Pořadí součástí detektoru v elektronové konfiguraci je následující: Micro Vertex Detector (MVD) a Silicon Tracking System (STS), oba usazené uvnitř silného supravodivého magnetu, Ring Imaging Cherenkov Detector (RICH), Transicion Radiation Detector (TRD, dva za sebou) Time of Flight (ToF) detektor, elektromagnetický kalorimetr (ECAL) a Projectile Spectator Detector (PSD) funguje jako hadronový kalorimetr. V mionové konfiguraci je odstraněn elektromagnetický kalorimetr a RICH detektor je nahrazen systémem mionových komnat (Muon Chambers system -- MUCH) Převzato z \cite{akishina4D}.}
    \label{fig:cbm_layout}
\end{figure}

CBM je jeden z vědeckých pilířů pracoviště FAIR. Jeho cílem bude studovat strukturu a stavovou rovnici baryonové hmoty o hustotách, které bychom mohli najít v jádrech neutronových hvězd. Experiment bude zkoumat srážky svazků protonů a těžkých iontů se statickým terčem při svazkových energiích od 2 do 45 AGeV (GeV na nukleon). \cite{akishina4D} Experiment bude zaměřen mimo jiné na \cite{cbm_physbook}:

\begin{itemize}
  \item{studium stavové rovnice jaderné hmoty o hustotách, které bychom mohli najít v neutronových hvězdách}
  \item{hledání fázové hranice mezi hadronovou fází a kvark-gluonovým plazmatem, případě oblasti fázové koexistence a koncového kritického bodu kvantové chromodynamiky.}
  \item{pátrání po jednoduchých a dvojitých hyperjádrech: těžkých objektech s nenulovým kvantovým číslem podivnosti}    
\end{itemize}

Experiment CBM je velmi rozsáhlý a složitý detektor, stávající se z mnoha podsystémů, my se však v rámci této práce budeme zabývat pouze jeho hlavním detektorem \textbf{STS}.

\subsection*{Supravodivý magnet}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/cbm_magnet.png}
    \caption{Nákres supravodivého magnetu \textbf{(A)} a znázornění intenzity magnetického pole v magnetu \textbf{(B)}. V bodě (0,0) je znázorněn terč detektoru, graf intenzity je vykreslen podél roviny YZ, $x=0$. Převzato z \cite{akishina4D}.}
    \label{fig:cbm_magnet}
\end{figure}

Účel silného supravodivého magnetu je \textit{ohýbat} nabité částice, u kterých pak podle jejich zakřivení dokážeme určit jejich náboj a energii. K tomu, aby bylo možné uvnitř magnetu uložit detektor STS bylo potřeba, aby měl magnet rozměry minimálně $1,3\times1,3m^2$ a generovat dostatečně silné magnetické pole, aby bylo magnetické pole uvnitř prvního detektoru \uv{co nejhomogennější}. Magned funguje na bázi velké Helmholzovy cívky, která opravdu uvnitř sebe sama vytváří magnetické pole, které se za určitých podmínek dá považovat za homogenní.

\subsection*{Silicon Tracking System}

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{obrazky-figures/STS.png}
    \caption{3D pohled na detektor STS a rozvržení jeho stanovišť \textbf{(A)}, diagram křemíkového proužkového detektoru \textbf{(B)} a schéma rozvržení šesté hladiny detektoru STS (světlejší barva proužku značí vytíženější oblasti detektoru) \textbf{(C)}. Převzato z \cite{akishina4D}}
    \label{fig:sts}
\end{figure}

Silicon Tracking System (STS) je primární detektor experimentu CBM a druhý detektor v pořadí po detektoru Minimum Vertex Detector (MVD). Cílem detektoru STS je vhodný odhad vlastností jednotlivých trajektorií částic, kterých může detektorem proletět během jedné \uv{události} při vzorkovací frekvenci 10 MHz až 700. Detektor je uspořádán do osmi, na podpůrných lešení z uhlíkových vláken upevněných \uv{hladin}: stanovišť, kdy první je od terče detektoru vzdáleno 30 centimetrů, mezi každou hladinou detektoru je pak vždy deset centimetrů prostor. Detektorové hladiny jsou uspořádány tak, aby pokrývaly polární úhly od 2,5\textdegree do 25\textdegree. Pro polární úhly menší než 2,5\textdegree prochází detektorem silný \uv{paprsek}\footnote{z angl. \textbf{Jet}} hadronů, u kterého je vzhledem k vysokým energiím a množství částic detekce nepraktická. 

STS má citlivé elektrody na každé straně detektoru, kde každá strana se stará o odečítání jedné ze dvou souřadnic, přitom proužky nesvírají pravý úhel, nýbrž úhel 15\textdegree, čímž se sice zmenší rozlišení detektoru podél jedné osy (y), zvýší se ale rozlišení druhé osy (x). Detektor je takto zařízen z toho důvodu, že nabité částice jsou ovlivněny magnetickým polem silného supravodivého magnetu, ve kterém se detektor nachází, v důsledku čehož mají částice uvnitř magnetu mnohem větší rozptyl ve směru osy $x$, než ve směru osy $y$. Algoritmus, který je v této práci použit pro rekonstrukci trajektorií částic, tento fakt zohledňuje upraveným modelem měření $H$:

\begin{equation}\label{eq:measurement_model}
    H = \begin{bmatrix}
        \cos{\alpha} & 0 & 0 & 0 & 0 \\
        0 & \sin{\alpha} & 0 & 0 & 0
    \end{bmatrix}
\end{equation}

kde $\alpha$ je úhel mezi jednotlivými osami detektoru, v případě STS $\alpha = 15$\textdegree. Model měření Kalmanova Filtru diskutujeme podrobně v kapitole \ref{track_reconstruction}.

\chapter{Teorie rekonstrukce trajektorií částic}\label{track_reconstruction}
Začátek následující podkapitoly byl převzat z disertační práce Dr. Akishiny \cite{akishina4D}. Díky neustálému vývoji v oborech \textbf{urychlovačové fyziky}\footnote{z angl. accelerator physics -- odnož aplikované fyziky, zabývající se návrhem, stavbou a provozem urychlovačů částic} a výpočetní techniky je dnes možné provozovat HEP experimenty při energiích a četnostech kolizí, které dříve nebyly možné. Tento fakt má ale také za důsledek, že zpracování stále se zvětšujících objemů experimentálních dat je čím dál komplikovanější úkol, na který již běžně dostupná výpočetní technika není vhodná.

Výrazně vyšší kolizní energie částic v HEP experimentech způsobují komplikovanější vzorce událostí a výrazně větší množství sledovaných částic v detektoru. Zároveň se zvyšuje také frekvence těchto událostí, což má za následek omezení času, který máme k dispozici na vyhodnocení každé takové události, tedy k rekonstrukci trajektorií jednotlivých částic a srážek mezi nimi. 

Vzorce srážek rekonstruujeme tak, že nejprve zrekonstruujeme trajektorie jednotlivých částic, které při průletu detektorem zanechávají měřitelné stopy na jednotlivých detektorových hladinách, tzv. \textbf{zásahy}. Proces rekonstrukce každé události probíhá v experimentu CBM ve Čtyřech fázích:

\begin{itemize}
    \item hledání trajektorií,
    \item prokládání trajektorií,
    \item identifikace částic podle proložených trajektorií,
    \item hledání částic s krátkou dobou života.
\end{itemize}

Přestože hledání a prokládání trajektorií částic spolu do značné míry souvisí, tradičně na ně pohlížíme jako na dvě odlišné fáze rekonstrukce. Zatímco hledání trajektorií je problém rozpoznávání vzorů, kde jednotlivé zásahy rozdělujeme do skupin, o kterých předpokládáme, že by mohly být vyvolány stejnou částicí. Po tom, co jednotlivé zásahy takto seskupíme, lze pro každou skupinu odhadnout nejpravděpodobnější parametry dráhy částice, a to i za přítomnosti velkého množství šumu. Tomuto postupu říkáme prokládání trajektorií.

\begin{figure*}[ht]\centering \label{fig:track_reconstruction}
  \centering
  \includegraphics[width=0.8\linewidth]{obrazky-figures/Track-reconstruction.png}
  \caption{První dvě fáze rekonstrukce trajektorií. Převzato z \cite{akishina4D}.}
  \label{fair_layout}
\end{figure*}

Proces hledání trajektorií pracuje zpravidla nad nezpracovanými, silně zašuměnými daty, která nemohla být předem nijak předzpracována ani zredukována. Z tohoto důvodu je proces hledání trajektorií považován za nejsložitější a časově nejnáročnější fázi rekonstrukce trajektorií.

Jednotlivé fáze rekonstrukce trajektorií se ale postupně stávají čím dál provázanějšími a čím silnější jsou souvislosti mezi nimi, tím jsou hranice mezi nimi mlhavější. V poslední době je běžné během první rekonstrukční fáze odhadovat jednotlivé parametry dráhy částice v každém kroku, díky čemuž získáváme s každým dalším měřením lepší povědomí o parametrech proložené dráhy, což nám dále umožňuje již v této první fázi vylučovat zásahy, u kterých je nepravděpodobné, že budou patřit k právě rekonstruované dráze. Tento trend je primárně určován existencí rekurzivních algoritmů užívajících metody Kalmanovy Filtrace. Tyto algoritmy jsou hlavním tématem této práce a metodami Bayesovských Filtrů se budeme dopodrobna zabývat v příští podkapitole. 

\section{Rekurzivní Bayesovský odhad}

Rekurzivní Bayesovský odhad, též Bayesovský Filtr, je obecný algoritmus pro výpočet odhadu stavu systému na základě jednotlivých měření systému:

\begin{equation}\label{eq:bayes_believed}
    bel\left(x_t\right) = p\left(x_t | z_{1:t}, u_{1:t}\right)
\end{equation}

kde $x_t$ značí stav systému v čase $t$, $bel(x_t)$ značí předpokládaný stav v čase $t$, $z_t$ značí měření v čase $t$ a $u_t$ značí \textbf{ovládání systému}\footnote{z angl. system control commands} v čase $t$. Stěžejní zákon pro odvození Bayesovského Filtru je \textit{Bayesův teorém}, který zní následovně:

\begin{equation}\label{eq:bayes_theorem}
    p\left(A|B\right) = \frac{p\left(B|A\right)p\left(A\right)}{p\left(B\right)}
\end{equation}

kde $A$ a $B$ jsou náhodné jevy, $p\left(A\right)$ je pravděpodobnost jevu $A$ a $p\left(A|B\right)$ je pravděpodobnost jevu $A$ za předpokladu, že je jev $B$ pravdivý. Pro účely odvození Bayesovských filtrů je užitečnější aplikace Bayesova teorému pro 3 náhodné jevy:

\begin{equation}\label{eq:bayes_theorem_3}
    p\left(A|B,C\right) = \frac{p\left(B|A,C\right)p\left(A|C\right)}{p\left(B|C\right)}
\end{equation}

Když dosadíme definici Bayesovského filtru (rovnice \ref{eq:bayes_believed}) do rovnice výše, dostaneme:

\begin{equation}\label{eq:bayes_2}
    p\left(x_t| z_t, z_{1:t-1}, u_{1:t}\right) = \frac{p\left(z_t|x_t, z_{1:t-1},u_{1:t}\right)p(x_t|z_{1:t-1},u_{1:t})}{p\left(z_t|z_{1:t-1}, u_{1:t}\right)}
\end{equation}

Pro lepší čitelnost zavedeme za jmenovatel rovnice \ref{eq:bayes_2} normalizační konstantu $\eta$:

\begin{equation}\label{eq:bayes_eta}
        \eta = \frac{1}{p\left(z_t|z_{1:t-1}, u_{1:t}\right)}
\end{equation}

\begin{equation}\label{eq:bayes_eta_whole}
        bel(x_t) = \eta \; p\left(z_t|x_t, z_{1:t-1}, u_{1:t}\right) \; p\left(x_t | z_{1:t-1}, u_{1:t}\right)
\end{equation}

Pro odvození rekurzivního vzorce Bayesovského filtru použijeme fakt, že skutečný stav systému $x$ je \textbf{Markovův proces se skrytými stavy}. Na rozdíl od klasických \textit{Markovových řetězců} neznáme vnitřní stav systému, o tomto stavu získáváme informace pouze za pomoci měření $z_k$. Protože je vnitřní stav systému dán Markovovým procesem, víme, že bude systém vykazovat \textit{Markovovu vlastnost}, která říká, že vývoj Markovova procesu je závislý pouze na jeho předchozím stavu. Díky této vlastnosti víme, že pokud známe stav $x_t$, k určení pravděpodobnosti měření $z_t$ nepotřebujeme žádné znalosti, protože pravděpodobnost $p(z_t|x_t)$ nijak neovlivní: 

\begin{equation}\label{eq:markov_property_1}
    p\left(z_t|x_t, z_{1:t-1}, u_{1:t}\right) = p\left(z_t|x_t\right)
\end{equation}

Rovnici \ref{eq:bayes_eta_whole} tedy můžeme zjednodušit následujícím způsobem:

\begin{equation}\label{eq:bayes_markov_1}
    bel\left(x_t\right) = \eta \; p\left(z_t|x_t\right)\; p\left(x_t|z_{1:t-1}, u_{1:t}\right)
\end{equation}

V dalším kroku aplikujeme \textbf{pravidlo úplné pravděpodobnosti}, které říká, že máme-li úplný systém jevů $\left\{B_n : n = 1,2,3,\dots\right\}$, pak pravděpodobnost libovolného jevu $A$ lze určit jako:

\begin{equation}\label{eq:total_probability_discrete}
    p(A|C) = \sum^n_{i=1} p(A|B_i, C) \; p(B_i| C)
\end{equation}

V našem případě vzorec úplné pravděpodobnosti (v jeho spojité podobě) umožní přidat k druhému členu rovnice \ref{eq:bayes_markov_1} další podmínku, kterou využijeme následně při odvození rekurzivního kroku filtru:

\begin{equation}\label{eq:total_probability_continuous}
    p\left(a|c\right) = \int^\infty_{-\infty} p(a|b, c) \; p(b | c) \; db
\end{equation}

Dosadíme-li do původní rovnice, dostaneme:

\begin{equation}\label{eq:bayes_total_probability}
    bel\left(x_t\right) = \eta \; p\left(z_t|x_t\right) \int p(x_t| x_{t-1}, z_{1:t-1}, u_{1:t}) \; p(x_{t-1} | z_{1:t-1}, u_{1:t}) \; dx_{t-1}
\end{equation}

Nyní můžeme zjednodušit člen $p\left(x_t|x_{t-1}, z_{1:t},u_{1:t}\right)$ pomocí Markovovy vlastnosti, díky které víme, že ke stanovení pravděpodobnosti, že se nacházíme ve stavu $x_t$, nám stačí znát předchozí stav systému $x_{t-1}$ a informaci o tom, jak se systém vyvíjí ze stavu $x_{t-1}$ do stavu $x_t$ (tuto znalost reprezentujeme jako $u_t$). Všechny ostatní znalosti jsou pro tento výpočet irelevantní, takže tento člen můžeme pomocí Markovovy vlastnosti přepsat jako:

\begin{equation}\label{eq:markov_2}
    p\left(x_t|x_{t-1}, z_{1:t}, u_{1:t}\right) = p\left(x_t|x_{t-1},u_t\right)
\end{equation}

Výsledná rovnice má pak následující tvar:

\begin{equation}\label{eq:bayes_markov_2}
    bel\left(x_t\right) = \eta \; p\left(z_t|x_t\right) \int p(x_t| x_{t-1}, u_t) \; p(x_{t-1} | z_{1:t-1}, u_{1:t}) \; dx_{t-1}
\end{equation}

V následujícím kroku zkoumáme člen $p(x_{t-1}|z_{1:t-1},u_{1:t})$: zkoumáme pravděpodobnost, že minulý stav systému byl $x_{t-1}$ za předpokladu, že známe veškerá předešlá měření až do času $t-1$, stejně tak víme, jak se systém vyvíjel v každém okamžiku až do času $t-1$. Zároveň je ale přítomna podmínka $u_t$, která říká, že víme, jakým způsobem se systém vyvine ze stavu $x_{-1}$ do stavu $x_t$. Můžeme předpokládat, že znalost, jakým způsobem se bude systém v budoucnosti vyvíjet, nemá na pravděpodobnost, že v čase $t-1$ bude systém ve stavu $x_{t-1}$ žádný vliv, a tuto podmínku ze členu $p(x_{t-1}|z_{1:t-1},u_{1:t})$ vyloučit:

\begin{equation}\label{eq:independence_assumption}
    p\left(x_{t-1}|z_{1:t-1},u_{1:t}\right) = p\left(x_{t-1}|z_{1:t-1},u_{1:t-1}\right)
\end{equation}

Dosazením získáváme tvar rovnice:

\begin{equation}\label{eq_bayes_identity}
    bel\left(x_t\right) = \eta \; p\left(z_t|x_t\right) \int p(x_t| x_{t-1}, u_t) \; p(x_{t-1} | z_{1:t-1}, u_{1:t-1}) \; dx_{t-1}
\end{equation}

Konečně, z původní definice odhadu $bel\left(x_t\right)$ v rovnici \ref{eq:bayes_believed} vyplývá, že člen 

$$p\left(x_{t-1}|z_{1:t-1},u_{1:t-1}\right)$$ 

můžeme přepsat jako:

\begin{equation}\label{eq:believed_recursion}
    p\left(x_{t-1}|z_{1:t-1},u_{1:t-1}\right) = bel\left(x_{t-1}\right)
\end{equation}

Čímž získáme potřebný rekurzivní krok a můžeme dosazením do rovnice \ref{eq:bayes_markov_2} dokončit odvození Bayesovského filtru:

\begin{equation}\label{eq:bayes_believed_final}
    bel\left(x_t\right) = \eta \; p\left(z_t|x_t\right) \int p(x_t| x_{t-1}, u_t) \; bel\left(x_{t-1}\right) \; dx_{t-1}
\end{equation}

Výsledný filtr můžeme dále rozepsat (a typicky rozepisujeme) na následující dva kroky:

\begin{itemize}
    \item predikční krok
    \begin{equation}\label{eq:bayes_prediction}
        \overline{bel}\left(x_t\right) = \int p(x_t | x_{t-1, u_t}) \; bel\left(x_{t-1}\right) \; dx_{t-1}
    \end{equation}
    \item korekční krok
    \begin{equation}\label{eq:bayes_correction}
        bel\left(x_t\right) = \eta \; p\left(z_t | x_t\right) \; \overline{bel}\left(x_t\right)
    \end{equation}
\end{itemize}

kde $\overline{bel}\left(x_t\right)$ je odhad $x_t$ na základě předchozích měření, $bel(x_t)$ je upravený odhad $x_t$ na základě nového měření, $p\left(k_t|u_t, x_{t-1}\right)$ v predikčním kroku nazýváme \textit{pohybový model} a $p\left(z_t|x_t\right)$ v korekčním kroku nazýváme \textit{model měření}.


\subsection*{Implementace Bayesovského Filtru}
Rekurzivní Bayesův odhad je obecný pravděpodobnostní přístup k rekurzivnímu odhadu neznámého rozložení pravděpodobnosti; v závislosti na typu řešeného problému vybíráme jednu z implementací Bayesovského filtru. Existuje více možných implementací Bayesovských filtrů, zejména pak:

\begin{itemize}
    \item \textbf{Kalmanův filtr (KF)}, optimální Bayesovský filtr pro filtrování nad normálními rozděleními a lineárními modely,
    \item \textbf{Rozšířený Kalmanův filtr (EKF)}, rozšíření Kalmanova filtru pro nelineární modely. Rozšířený Kalmanův filtr je v praxi nejrozšířenější metoda pro aproximaci pohybu částic.
    \item \textbf{Částicové filtry}\footnote{z angl. Particle Filters}, jinak také metody SMC\footnote{SMC -- Sequential Monte Carlo (sekvenční Monte Carlo metody)}, pracují nad nelineárními modely s arbitrárními rozloženími pravděpodobností. Díky tomu jsou ze zde zmíněných metod použitelné pro nejširší spektrum problémů, pro naše účely jsou ale z důvodu nízké rychlosti nevhodné.
\end{itemize}

Právě Kalmanovými filtry a EKF se budeme v další podkapitole podrobně zabývat.

\section{Kalmanův Filtr}\label{kalman_filter}

\textbf{Kalmanův Filtr}, známý také jako metoda \textbf{LQE}\footnote{LQE -- Linear Quadratic Estimation (lineárně-kvadratický odhad)}, je implementace Bayesova Filtru, která poskytuje optimální odhad stavu dynamického systému za předpokladu, že:

\begin{itemize}
    \item pohybový model i model měření jsou lineární
    \item všechny náhodné proměnné v systému mají \textit{normální rozdělení}
\end{itemize}

Metoda slouží k získání optimálního odhadu $\boldsymbol{r}$ neznámého stavového vektoru $\boldsymbol{r}^t$ za pomocí \textit{vektorů měření} (jinak také \textit{pozorování}) $\boldsymbol{m}_k, k = 1 \dots n$ vektoru $\boldsymbol{r}^t$. Metoda je inicializována volbou počátečního odhadu $\boldsymbol{r} = \boldsymbol{r}_0$\footnote{U EKF musíme být narozdíl od klasického KF při volbě $\boldsymbol{r}_0$ opatrní, metoda může při špatně zvolené počáteční hodnotě divergovat} a postupně s každým dalším měřením zpřesňuje odhad $\boldsymbol{r}$ stavového vektoru $\boldsymbol{r}^t$. \cite{simd_kalman_filter}.

Hodnota vektoru $\boldsymbol{r}^t$ se může v čase měnit: jeho změnu u lineárních modelů reprezentujeme maticí $F_k$, která obsahuje informaci o tom, jak se v čase mění každá složka stavového vektoru $\boldsymbol{r}^t$. Stav systému se tak mezi jednotlivými měřeními mění, nemusíme ale přesně vědět, jak a/nebo kdy k těmto změnám dochází. Veškeré náhodné jevy mezi měřeními $\boldsymbol{m}_{k-1}$ a $\boldsymbol{m}_k$, které nedokážeme modelovat v rámci propagační matice $F_k$, modelujeme pomocí takzvaného \textit{procesního šumu}, který je reprezentován vektorem $\boldsymbol{\nu}_k$. Výsledná rovnice změny stavového vektoru $\boldsymbol{r}^t$ vypadá následovně:

\begin{equation}\label{eq:kalman_rt_prediction}
    \boldsymbol{r}^t_k = F_k\boldsymbol{r}^t_{k-1}+\boldsymbol{\nu}_k
\end{equation}

kde $F_k$ je lineární propagační operátor a $\boldsymbol{\nu}_k$ je vektor procesního šumu mezi měřeními $\boldsymbol{m}_{k-1}$ a $\boldsymbol{m}_k$. V případě Kalmanova Filtru předpokládáme, že měření $\boldsymbol{m}_k$ je lineárně závislé na stavu systému $\boldsymbol{r}^t_k$, tedy říkáme, že \textbf{model měření} $H_k$ je lineární:

\begin{equation}\label{eq:kalman_measurement_model}
    \boldsymbol{m}_k = H_k\boldsymbol{r}^t_k + \boldsymbol{\eta}_k
\end{equation}

kde $\boldsymbol{\eta}_k$ je vektor chyby měření $\boldsymbol{m}_k$, $H_k$ je \textit{model měření}. Model měření je matice s dimenzí $m \times n$, kde $m$ je rovno velikosti stavového vektoru $\boldsymbol{r}^t$ a $n$ je rovno velikosti vektoru měření $\boldsymbol{m}_k$. Model měření určuje, jaký je vztah mezi jednotlivými měřeními systému a jeho reálným stavem.
Kalmanův Filtr předpokládá, že chyba měření a procesní šum nijak nekorelují a že známe jejich \textbf{kovarianční matice}:

\begin{equation}\label{eq:kalman_covariance}
    \begin{split}
        \left\langle \boldsymbol{\eta}_i\right\rangle = \left\langle\boldsymbol{\nu}_j\right\rangle &= \boldsymbol{0},\\
        \left\langle \boldsymbol{\eta}_k \cdot \boldsymbol{\eta}_k^T \right\rangle &= V_k,\\
        \left\langle \boldsymbol{\nu}_k \cdot \boldsymbol{\nu}_k^T \right\rangle &= Q_k.
    \end{split}
\end{equation}

kde notace $\left\langle \boldsymbol{X} \right\rangle = E\left(\boldsymbol{X}\right)$ značí střední hodnotu náhodné veličiny $\boldsymbol{X}$.
\textbf{Kovarianční matice} (také variačně-kovarianční matice) reálné vícerozměrné náhodné veličiny $X$ o $n$ složkách je čtvercová matice $\left(n \times n\right)$, která obsahuje \textbf{kovarianci} $i$-té a $j$-té složky náhodné veličiny $\boldsymbol{X}$:

$$
    X = 
    \begin{bmatrix}
    \sigma_1^2 & \sigma_1\sigma_2 & \cdots & \sigma_1\sigma_n \\
    \sigma_2\sigma_1 & \sigma_2^2 & \cdots & \sigma_2\sigma_n \\
    \vdots & \vdots & \ddots & \vdots \\
    \sigma_n\sigma_1 & \sigma_n\sigma_2 & \cdots & \sigma_n^2\\
    \end{bmatrix}
$$

kde $\sigma_i\sigma_j = cov(X_i, X_j) = E\left[\left(X_i - E[X_i]\right)\left(X_j - E[X_j]\right)\right]$ je kovariance veličin $X_i$ a $X_j$ a $E[X_i] = \langle X_i \rangle$ je střední hodnota veličiny $X_i$. 

Klasický Kalmanův Filtr běžně členíme do čtyř základních kroků:

\begin{figure}[ht]
   \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/Kalman_filter_diagram.png}
    \caption{Blokové schéma klasického Kalmanova Filtru. Převzato z \cite{simd_kalman_filter}}
    \label{fig:kalman_filter_diagram}
\end{figure}

\begin{enumerate}
    \item \textbf{inicializace}, kdy zvolíme počáteční hodnoty odhadu stavového vektoru $\boldsymbol{r}^+_0$ a kovarianční matice $C^+_0$,
    \item \textbf{predikce}, kdy pomocí propagačního operátoru $F_k$ stanovíme odhad dalšího stavu stavového vektoru a kovarianční matice,
    \item \textbf{zpracování procesního šumu}, kdy upravíme odhad příští kovarianční matice $C^-_k$ pomocí kovarianční matice procesního šumu $Q_k$,
    \item \textbf{filtrace}, kdy pomocí nového měření upravíme odhad stavového vektoru a kovarianční matice tak, aby příští odhady byly přesnější.
\end{enumerate}

\paragraph{Inicializace} Během inicializace Kalmanova Filtru musíme určit první aproximaci stavového vektoru $\boldsymbol{r}^+_0$ a kovarianční matice $C^+_0$: Mějme reálný stav systému, vyjádřený vektorem $\boldsymbol{r}^t$ a jeho odhad $\boldsymbol{r}$. Odhad stavu systému $\boldsymbol{r}$ má vždy konečnou přesnost: zavedeme vektor chyby $\boldsymbol{\xi}$, který vyjadřuje rozdíl mezi reálným stavem systému a jeho odhadem a kovarianční matici $C$, které definujeme následovně:

\begin{equation}\label{eq:error_xi}
    \begin{split}
        \boldsymbol{r} &= \boldsymbol{r}^t + \boldsymbol{\xi} \\
        C &= \left\langle \boldsymbol{\xi} \cdot \boldsymbol{\xi}^T \right\rangle
    \end{split}
\end{equation}

Podle dostupných informací zvolíme počáteční odhad vektoru $\boldsymbol{r}^+_0$, a jeho kovarianční matici potřebujeme zvolit tak, aby reprezentovala fakt, že před prvním měřením nemáme o stavu systému žádné znalosti. Kovarianční matici $C_0$ proto zvolíme následovně:

\begin{equation}\label{eq:covariance_matrix_C_0}
    C^+_0 = \boldsymbol{I} \cdot \text{inf}^2
\end{equation}

kde $\boldsymbol{I}$ je jednotková matice a \textit{inf} značí velké kladné číslo.

\paragraph{Predikce a zpracování procesního šumu} Pokud očekáváme, že se stav systému mezi jednotlivými měřeními změní, musí se odpovídajícím způsobem změnit i jeho odhad a kovarianční matice: Pomocí propagačního operátoru $F_{k-1}$ získáme predikci stavu systému $\boldsymbol{r}^-_k$ a jeho kovarianční matice $C^-_k$:

\begin{equation}\label{eq:prediction_step}
    \begin{split}
        \boldsymbol{r}^-_k &= F_{k-1} \boldsymbol{r}^+_{k-1} \\
        C^-_k &= F_{k-1}C^+_{k-1}F_{k-1}^T + Q_k
    \end{split}
\end{equation}

kde operátor $F_{k-1}$ modeluje vývoj systému mezi stavy $\boldsymbol{r}_{k-1}$ a $\boldsymbol{r}_k$, $F^T_{k-1}$ je transponovaná matice $F_{k-1}$, $Q_k$ je kovarianční matice procesního šumu, kterou metoda Kalmanova Filtru předpokládá, že předem známe. Zatímco operátor $F_{k-1}$ modeluje deterministické změny uvnitř systému, přidáním kovarianční matice procesního šumu $Q_k$ zohledňujeme stochastické procesy, které popsat v pohybovém modelu je buď nemožné, nebo zbytečně složité.

Zatímco $\boldsymbol{r}_k^-$ a $C^-_k$, které vypočítáváme v rámci predikčního kroku na základě minulých měření $\boldsymbol{m}_i$, $i = 1\dots k-1$, jsou predikce hodnot $\boldsymbol{r}_k$ a $C_k$, $\boldsymbol{r}_{k-1}^+$ a $C^+_{k-1}$ jsou již opravené (resp. filtrované) predikce z předchozí iterace Kalmanova Filtru. Vypočítáváme je v rámci filtračního kroku pomocí hodnoty nejnovějšího měření $m_k$ (viz. níže). 

V případě, že jsme dosud neprovedli žádné měření, propagujeme původní odhad stavu systému $\boldsymbol{r}_0^+$ a jeho kovarianční matici $C_0^+$, které jsme určili během inicializačního kroku:

\begin{equation}\label{eq:prediction_propagation}
    \begin{split}
        \boldsymbol{r}^-_1 &= F_0\boldsymbol{r}^+_0\\
        C^-_1 &= F_0C^+_0F_0^T+Q_1
    \end{split}
\end{equation}

\paragraph{Filtrace} V závěrečném kroku Kalmanova filtru upravíme odhad hodnoty $\boldsymbol{r}^-_{k-}$ (a její kovarianční matice $C_{k-1}^-$) pomocí nejnovějšího měření $\boldsymbol{m}_k$, díky čemuž získáme nové, optimální odhady těchto veličin $\boldsymbol{r}^+_k$ a $C^+_k$. Zároveň v tomto kroku vypočítáme $\chi^2$ odchylku odhadu $r^+_k$ od dosud proběhlých měření $\boldsymbol{m}_i$, $i=1\dots k$, díky čemuž dokážeme pro každou iteraci Kalmanova Filtru určit, jak dobře odpovídá odhad měřením, což je pro naše účely velmi užitečné, protože díky $\chi^2_k$ odchylce můžeme rychle vypočítat, které zásahy v detektoru s nejvyšší pravděpodobností patří ke stejné trajektorii.

Pro zjednodušení dalších výpočtů (rovnice \ref{eq:kalman_filtration} a \ref{eq:kalman_chi_square}) zavedeme nejprve pomocný člen -- takzvané \textbf{reziduum} (resp. reziduál, z anglického \textit{residual}, v praxi používáno nejčastěji), korekční člen, který vyjadřuje rozdíl mezi odhadem Kalmanova Filtru $\boldsymbol{r}_k^-$ a vlastním měřením $\boldsymbol{m}_k$.

\begin{equation}\label{eq:residual}
    \boldsymbol{\zeta}_k = \boldsymbol{m}_k - H_k\boldsymbol{r}_k^-
\end{equation}

Abychom mohli aktualizovat odhad stavu $\boldsymbol{r}^+_k$, potřebujeme spočítat takzvanou \textit{matici Kalmanova zesílení}\footnote{z angl. Gain Matrix} $K_k$:

\begin{equation}\label{eq:kalman_gain}
    K_k = C_k^-H_k^T\left(V_k + H_kC_k^-H_k^T\right)^{-1}
\end{equation}

kde $C_k^-$ je v predikčním kroku předpokládaná hodnota kovarianční matice, $V_k$ je kovarianční matice chyby měření a $H_k$ je model měření. Rovnici Kalmanova zesílení můžeme ještě drobně zjednodušit tak, že zavedeme takzvanou \textit{váhovou matici} $W_k$:

\begin{equation}\label{eq:kalman_weight_matrix}
    W_k = \left(V_k + H_kC_k^-H_k^T\right)^{-1}
\end{equation}

Rovnice Kalmanova zrychleni pak bude mít následující tvar:

\begin{equation}\label{eq:kalman_gain_weight}
    K_k = C_k^-H_k^TW_k
\end{equation}

 Jakmile máme Kalmanovo zesílení spočítané, můžeme aktualizovat odhad stavového vektoru $\boldsymbol{r}^+_k$ (a jeho kovarianční matice):

\begin{equation}\label{eq:kalman_filtration}
    \begin{split}
        \boldsymbol{r}^+_k &= \boldsymbol{r}^-_k + K_k \boldsymbol{\zeta}_k\\
        C_k^+ &= \left(\boldsymbol{I} - K_kH_k\right)\cdot C_k^-
    \end{split}
\end{equation}

kde $\boldsymbol{I}$ je jednotková matice a $\boldsymbol{m}_k$ je poslední měření systému. V závěrečné fázi spočítáme novou odchylku $\chi^2_k$ na základě posledního měření $\boldsymbol{m}_k$:

\begin{equation}\label{eq:kalman_chi_square}
    \chi^2_k = \chi^2_{k-1} + \boldsymbol{\zeta}_k^TW_k \boldsymbol{\zeta}_k
\end{equation}


V tuto chvíli získaná hodnota $\boldsymbol{r}_k^+$ s kovarianční maticí $C_k^+$ je optimální odhad současného stavu systému. Tato dvojice je zároveň výstup dokončené iterace a vstup následující iterace Kalmanova Filtru. Hodnota $\chi^2_k$ nám později pomocí celulárního automatu pomůže ke hledání zásahů, které s nejvyšší pravděpodobností patří ke trajektorii stejné částice. 

 

\subsection*{Kalmanův filtr v experimentu CBM}
Následující podsekce byla převzata z disertační práce Dr. Akishiny \cite{akishina4D}. Výsledkem hledání trajektorií (fáze 1: Track finding, viz. kapitola \ref{track_reconstruction}) by měly být sady měření detektoru uskupené do souborů zásahů, které by v ideálním případě byly všechny vytvořené stejnou částicí. Úkolem fáze 2 (Track fitting) rekonstrukce trajektorií je odhadnout parametry trajektorií a jejich chyby, aby bylo možné získat kinematické vlastnosti částic, které nám později umožní rekonstrukci částic s krátkou \textit{střední dobou života} a nakonec i analýzu fyzikálních procesů, které v detektoru proběhly.

Algoritmus ve fázi 2 využívá tzv. \textit{modelu trajektorií}, což je teoretický předpoklad o pohybových rovnicích pro nabité částice v rámci objemu detektoru. Na pohyb částic v detektoru má vliv vícero jevů, které nelze v rámci modelu trajektorií jednoduše zohlednit, zejména se jedná o:

\begin{itemize}
    \item mnohonásobný rozptyl,
    \item ionizace částic a
    \item ztráty energie způsobené radiací.
\end{itemize}
 
 Tyto vlivy zavádějí do rekonstrukce kinematických vlastností částic perturbace a celý proces tak narušují. Použití správné rekonstrukční metody nám umožňuje dopad těchto rušivých jevů na měření minimalizovat. V dnešní době se ve většině moderních HEP experimentů používá pro odhad parametrů trajektorií částic Kalmanův Filtr.
 
 I když se parametry trajektorií dají z jednotlivých měření odvodit pomocí \textit{metody nejmenších čtverců}, v praxi upřednostňujeme Kalmanův filtr z důvodu, že díky své rekurzivní povaze je výpočetně mnohem výhodnější. I když obě metody počítají nad maticemi, zatímco mocnost matice u metody nejmenších čtverců závisí na počtu měření, u Kalmanova filtru na počtu stupňů volnosti systému.

\subsection*{Nedostatky konvenčního Kalmanova Filtru}

Konvenční Kalmanův Filtr je určený pro lineární dynamické systémy s normálními rozděleními pravděpodobnosti, což znamená, že vývoj systému mezi jednotlivými měřeními $\boldsymbol{m}_{k-1}$ a $\boldsymbol{m}_k$ dokážeme popsat pomocí \textbf{lineární} rovnice (viz rovnice \ref{eq:kalman_rt_prediction}). V případě, že systém pracuje s lineárními modely a pouze s normálními rozloženími pravděpodobnosti, Kalmanův Filtr je \textbf{optimální metoda aproximace} stavu takového systému.

Bohužel, v praxi se často setkáváme se situacemi, kdy model závislostí mezi jednotlivými měřeními nelze popsat pomocí lineární funkce, případně stochastické jevy systému nelze modelovat jako normální rozdělení pravděpodobnosti. V případě nelineárních modelů (případ užití, který platí i pro aproximaci trajektorií částic v detektoru CBM) je potřeba nejprve najít takový lineární model, který bude co nejlépe aproximovat chování původního, nelineárního modelu. K nalezení takového modelu nám pomůže speciální případ aproximace pomocí \textbf{Taylorova rozvoje}, kterému říkáme \textbf{linearizace}.

\subsection*{Taylorův rozvoj}
Mějme nepolynomiální funkci $f : \mathbb{R} \to \mathbb{R}$, (například $f(x) = \sin (x)$): Pokud potřebujeme použít takovou funkci a z jakéhokoli důvodu je pro nás nepraktické vypočítávat její výstup přesně (například pokud potřebujeme vypočítat mnoho hodnot pro vzájemně blízká $x$ a pokaždé vypočítávat přesnou hodnotu je příliš výpočetně náročné), můžeme funkci vhodně \textit{aproximovat} a dále používat aproximovanou funkci, která je pro náš případ užití výhodná.

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/taylor_series_graph.png}
    \caption{Funkce $y = \sin (x)$ a její Taylorovy rozvoje prvního, třetího, pátého, sedmého a devátého stupně. Každý vyšší stupeň Taylorova rozvoje aproximuje původní funkci přesněji. Vytvořeno pomocí nástroje Desmos \cite{desmos}.}
    \label{fig:sin_x_taylor_series}
\end{figure}

\textbf{Taylorova řada} je mocninná řada, pomocí které můžeme za určitých předpokladů vyjádřit nepolynomiální funkci v okolí bodu $x_{lin}$ jako polynom. Taylorovu řadu reálné, \textbf{hladké} funkce $f(x)$ zapisujeme následujícím způsobem:

\begin{equation}\label{eq:taylor_series}
    f(x) = f(x_{lin}) + \frac{f'(x_{lin})}{1!}(x - x_{lin}) + \frac{f''(x_{lin})}{2!}(x - x_{lin})^2 + \dots = \sum_{k=0}^{\infty} \frac{f^{(k)}(x_{lin})}{k!}(x - x_{lin})^k
\end{equation}

kde $f^{(n)}$ značí $n$-tou derivaci funkce $f$. Funkce je \textbf{hladká} tehdy, kdy nejenom funkce samotná, ale i všechny její derivace jsou spojité.

Uveďme jako příklad funkci $f(x) = \sin x$, kterou chceme aproximovat Taylorovým polynomem v bodě 0 (v případě Taylorovy řady v bodě 0 mluvíme o tzv. \textit{Maclaurinově řadě}):

\begin{equation}\label{eq:taylor_sin_calculation}
    \begin{split}
        T_1^{f,0} = f(0) + \frac{f'(0)}{1!} (x - 0)^1 = 0 + \cos (0)x &= x \\
        T_2^{f,0} = T_1^{f,0} + \frac{f''(0)}{2!} (x - 0)^2 = x - \frac{\sin (0)}{2} x^2 &= x \\
        T_3^{f,0} = T_2^{f,0} + \frac{f^{(3)}(0)}{3!} (x - 0)^3 = x - \frac{\cos (0)}{6} x^3 &= x - \frac{x^3}{6} \\
        T_4^{f,0} = T_3^{f,0} + \frac{f^{(4)}(0)}{4!} (x - 0)^4 = x - \frac{x^3}{6} + \frac{\sin (0)x^4}{4!} &= x - \frac{x^3}{6}
    \end{split}
\end{equation}

Můžeme si všimnout, že sudé stupně Maclaurinovy řady se u funkce $f(x) = \sin x$ vždy rovnají nule, díky čemuž víme, že nám pro její aproximaci stačí liché stupně rozvoje:

\begin{equation}\label{eq:taylor_sin_simplified}
    \sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \dots = \sum_{n=0}^\infty \left(-1\right)^n \frac{x^{2n+1}}{\left(2n + 1\right)!}
\end{equation}

pro $x \in \left(-\infty, \infty\right)$. Analogicky můžeme zapsat například i Maclaurinův rozvoj funkce $f(x) = \cos x$, kde jsou nenulové naopak sudé členy rozvoje. Jednotlivé polynomy jsou ukázány na obrázku \ref{fig:sin_x_taylor_series}.

\textbf{Linearizace} je speciální případ aproximace Taylorovým rozvojem, kdy aproximujeme pouze prvním stupněm Taylorova polynomu (tj. přímkou):

\begin{equation}\label{eq:linearization}
    L(x)_f = f(x^{lin}) + f'(x^{lin})(x - x^{lin})
\end{equation}

Taková lineární funkce nám poskytuje vyhovující aproximaci pro $x$ blízká $x^{lin}$. Mějme funkci $f(x) = \sqrt{x}$, a hledejme hodnotu pro $\sqrt{4,41}$: Víme, že $\sqrt{4} = 2$, a tak vhodně zvolíme $x^{lin} = 4$:

\begin{equation}\label{eq:linearization_sqrt}
    L(x) = \sqrt{4} + \frac{1}{2\sqrt{4}}(x-4) = 1 + \frac{x}{4}
\end{equation}

Linearizovaná funkce $L(x)$ aproximuje $\sqrt{4,41} \approx 2,1025$, přesné řešení je přitom $\sqrt{4,41} = 2,1$. Linearizace nám tedy poskytuje solidní aproximaci pomocí relativně jednoduchého výpočtu. Pokud ale chceme linearizaci použít pro funkce, které pracují nad vektory, narazíme na problém, jak získat první derivaci matice. K tomuto účelu nám poslouží tzv. \textbf{Jacobiho matice}.

\subsection*{Jacobiho matice}

Mějme funkci $\boldsymbol{f}: \mathbb{R}^n \to \mathbb{R}^m$. Jacobiho matice $J$ je \textbf{matice parciálních derivací} takové funkce ve tvaru $m \times n$:

\begin{equation}\label{eq:jacobian_matrix}
    J = \begin{bmatrix}
        \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ 
        \frac{\partial f_2}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\
        \vdots & \vdots & \ddots & \vdots \\
        \frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n} \\
    \end{bmatrix}
\end{equation}

Jacobiho matice nám pomůže provádět linearizaci nelineárních modelů pro použití v rámci Kalmanova Filtru. Této verzi algoritmu pak říkáme \textbf{Rozšířený Kalmanův Filtr (EKF)}\footnote{EKF -- Extended Kalman Filter (rozšířený Kalmanův Filtr)}.

\subsection*{Rozšířený Kalmanův filtr}

\textbf{Rozšířený Kalmanův Filtr}, jak název napovídá, je rozšíření původní metody Kalmanova Filtru, které pracuje nad nelineárními modely měření a/nebo modely vývoje systému. 

Rozšířený Kalmanův Filtr sice umožňuje pracovat s nelineárními modely, z tohoto faktu ale plyne také jeho zásadní nevýhoda: Kvůli procesu linearizace není EKF optimální algoritmus pro odhad stavu systému, a při nevhodně zvolených hodnotách počátečního odhadu systému nebo linearizačního bodu $\boldsymbol{r}_{k-1}^{lin}$ může metoda divergovat, je tedy velmi důležité tyto hodnoty zvolit vhodně.

V případě rozšířeného Kalmanova Filtru musíme některé rovnice pozměnit tak, aby mohly pracovat nad nelineárními systémy. V případě, že pracujeme s nelineárním modelem vývoje systému $\boldsymbol{f}_k$, potřebujeme rovnici \ref{eq:kalman_rt_prediction} upravit následujícím způsobem:

\begin{equation}\label{eq:ekf_evolution_model}
    \boldsymbol{r}_k^- \equiv \boldsymbol{f}_k(\boldsymbol{r}_{k-1}^+) \approx \boldsymbol{f}_k(\boldsymbol{r}_{k-1}^{lin})+F_k(\boldsymbol{r}_{k-1}^+ - \boldsymbol{r}_{k-1}^{lin})
\end{equation}

kde $\boldsymbol{f_k}$ je \textbf{nelineární} model vývoje systému, $\boldsymbol{r}_{k-1}^{lin}$ je linearizační bod, typicky zvolený podle minulého stavu systému, a $F_k$ je Jacobiho matice nelineárního modelu vývoje $\boldsymbol{f}_k$:

\begin{equation}\label{eq:ekf_evolution_model_jacobian}
    F_{k(ij)} = \left. \frac{\partial \boldsymbol{f}_k(\boldsymbol{r}_{k-1}^+)_{(i)}}{\partial \boldsymbol{r}_{k-1(j)}^+}\right|_{\boldsymbol{r}_{k-1}^+=\,\boldsymbol{r}_{k-1}^{lin}}
\end{equation}

Podobně potřebujeme upravit rovnici vztahu mezi měřením a reálným stavem systému (rovnice \ref{eq:kalman_measurement_model}) v případě nelineárního modelu měření:

\begin{equation}\label{eq:ekf_measurement_model}
    \boldsymbol{m}_k(\boldsymbol{r}_k^t) \equiv \boldsymbol{h}_k(\boldsymbol{r_k^t}) + \boldsymbol{\eta}_k \approx \boldsymbol{h}_k(\boldsymbol{r}_k^{lin}) + H_k(\boldsymbol{r}_k^t - \boldsymbol{r}_k^{lin}) + \boldsymbol{\eta}_k
\end{equation}

Kde $H_{k(ij)}$ je Jacobiho matice modelu měření v následujícím tvaru:

\begin{equation}\label{eq:ekf_measurement_model_jacobian}
    H_{k(ij)} = \left. \frac{\partial\boldsymbol{h}_k(\boldsymbol{r}_k)_{(i)}}{\partial\boldsymbol{r}_{k(j)}}\right|_{\boldsymbol{r}_{k-1}^+=\,\boldsymbol{r}_{k-1}^{lin}}
\end{equation}

V případě nelineárního modelu měření se také změní rovnice pro výpočet rezidua:

\begin{equation}\label{eq:ekf_residual}
    \boldsymbol{\zeta} = \boldsymbol{m}_k - (\boldsymbol{h}_k (\boldsymbol{r}_k^{lin}) + H_k (\boldsymbol{r}_k^- - \boldsymbol{r}_k^{lin}))
\end{equation}

Již jsme zmínili, že rozšířený Kalmanův filtr může mít v některých případech problémy se stabilitou. Z použití rozšířeného Kalmanova filtru vyplývají některé potíže s implementací stabilního a efektivního algoritmu a jeho následnou akcelerací na různých hardwarových platformách. Některými možnými řešeními se budeme dále zabývat v kapitole \ref{implementation}.

\section{Metody rekonstrukce trajektorií}

Časově nejnáročnější část rekonstrukce trajektorií je jejich rekonstrukce, proto je potřeba vhodná volba sledovacího algoritmu. Náročnost této části analýzy událostí v detektoru vychází z faktu, že v této fázi algoritmus pracuje s hrubými, dosud nijak nezpracovanými daty a objem dat, který musí algoritmus zpracovat nelze nijak omezit. Často se tedy stává, že algoritmus musí v této fázi zpracovávat data až o objemu jednotek TB/s. Zároveň je časté, že v závislosti na konkrétním experimentu vyžadujeme, aby byla data zpracována v reálném čase, což úkol vývoje rekonstrukčního programu dále komplikuje a je tedy velmi důležité takový algoritmus vhodně zvolit. Situaci dále komplikuje fakt, že každý detektor je jiný, je tedy v případě každého detektoru potřeba výběr algoritmu přizpůsobit konkrétním vlastnostem detektoru. V následující sekci shrneme několik takových algoritmů, které můžeme pro rekonstrukci trajektorie zvolit:

\begin{itemize}
    \item Houghova transformace,
    \item celulární automaty.
\end{itemize}

Právě metodou rekonstrukce trajektorií pomocí celulárních automatů se pak budeme podrobněji zabývat v kontextu experimentu CBM.

Důležitým faktorem při výběru vhodného algoritmu je otázka, jak moc je možné problém rekonstrukce trajektorií zjednodušit, tedy jaké zjednodušující předpoklady, které můžeme zavést, usnadní (popřípadě umožní zrychlení např. paralelizací algoritmů tam, kde to u obecného problému bylo nemožné nebo nepraktické) řešení daného problému, aniž by (příliš) snížil přesnost výpočtu. Pro nás je pravděpodobně nejdůležitější vědět, jestli pro rekonstrukci trajektorií daná metoda pracuje se všemi vstupními daty, nebo jestli dokáže pracovat i pouze nad nějakou jejich podmnožinou. Pokud metoda vyžaduje přístup ke všem vstupním datům, říkáme, že mluvíme o \textbf{globální} metodě, v opačném případě mluvíme o \textbf{lokální} metodě. Pro naše účely je výhodné pracovat s metodami lokálními, protože se přirozeně nabízí možnost metodu do jisté míry paralelizovat. Níže se budeme zabývat Houghovou transformací jako zástupcem globálních metod rekonstrukce částic a metodou celulárních automatů jako zástupcem lokálních metod. Metodou celulárních automatů se budeme také dále zabývat v kontextu rekonstrukce trajektorií.

\subsection*{Houghova transformace}

Houghova transformace je technika extrakce prvků, původně navržena Paulem V. C. Houghem za účelem strojové analýzy fotografií bublinkových komor \cite{Hough1959MachineAO}. Ve své generalizované formě \cite{duda_generalized_hough} je hojně používaná v analýze obrazu, počítačovém vidění nebo zpracování obrazu. Cílem metody je najít nedokonale vyobrazené objekty (zpravidla objekty stejného charakteru, například přímky nebo kružnice) pomocí takzvaného \textbf{hlasování}.

Mějme nejjednodušší verzi Houghovy transformace, a to takovou, která detekuje v obraze přímky: Přímku definujeme jednoduše pomocí rovnice $$y = mx + b$$ v důsledku čehož můžeme přímku zapsat jako uspořádanou dvojici parametrů $(m, b)$, což je informace, která jednoznačně definuje danou přímku jako bod $(m,b)$ v takzvaném \textbf{parametrovém prostoru}\footnote{z angl. parameter space}. 

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.30\textwidth]{obrazky-figures/hough_image_space_points.png}\hfill
    \includegraphics[width=0.30\textwidth]{obrazky-figures/hough_parameter_space.png}\hfill
    \includegraphics[width=0.30\textwidth]{obrazky-figures/hough_image_space.png}
    \caption{Příklad Houghovy transformace. Body v \textbf{reálném prostoru} (obrázek vlevo) jsou převedeny do \textbf{parametrového prostoru} jako přímky. Pokud se všechny body nachází na jedné přímce v reálném prostoru, budou mít jim odpovídající přímky v parametrovém prostoru právě jeden průsečík (obrázek uprostřed). Průsečík v parametrovém prostoru obsahuje parametry přímky, na které se všechny tyto body nachází (obrázek vpravo). Obrázky převzaty z \cite{hough_watershed_transforms}.}
    \label{fig:hough_transform_line_example}
\end{figure}

Uveďme jako příklad hledání přímek pomocí Houghovy transformace následující situaci, znázorněnou na obrázku \ref{fig:hough_transform_line_example}: mějme přímku v reálném prostoru, definovanou funkcí ${y = 0.5x + 1.5}$, o které máme informaci v podobě čtyř bodů: $A_1 = (1,2)$, $B_1 = (3,3)$, $C_1 = (5,4)$, $D_1 = (7,5)$. Tak, jako můžeme zapsat jakoukoli přímku jako bod v parametrovém prostoru, můžeme obdobným způsobem pracovat s jednotlivými body. Každý bod z reálného prostoru se promítne do parametrového prostoru jako přímka, která obsahuje body s parametry všech přímek, které v reálném prostoru procházejí daným bodem. Všechny čtyři body tak vyneseme jako přímky v parametrovém prostoru a pozorujeme jejich podobu: Z povahy Houghovy transformace vyplývá, že všechny přímky v parametrovém prostoru, které mají v parametrovém prostoru společný průsečík, leží v reálném prostoru na stejné přímce. Průsečík těchto přímek v parametrovém prostoru nám zároveň sděluje, jaké jsou parametry takové přímky, v případě příkladu na obrázku \ref{fig:hough_transform_line_example} bodu ${(0.5, 1.5)}$ v parametrovém prostoru odpovídá přímka ${y = 0.5x + 1.5}$, která prochází všemi body $A_1$ až $D_1$.

Tento přístup má ale nevýhodu v případě, že chceme takto zapsat svislou přímku, což by mělo za následek parametr $m$ jdoucí k nekonečnu. Když se pokusíme převést do parametrového prostoru body, které jsou přímo nad sebou (tj. mají stejnou $x$ souřadnici, viz obrázek \ref{fig:hough_transform_vertical_line_problem}), zjistíme, že jsou všechny odpovídající přímky rovnoběžné a přestože lze všemi body vést jednu přímku, v parametrovém prostoru žádný průsečík nenajdeme.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.45\textwidth]{obrazky-figures/hough_vertical.png}\hfill
    \includegraphics[width=0.45\textwidth]{obrazky-figures/hough_vertical_parameter.png}
    \caption{Příklad nedostatku Houghovy transformace. Obrázek vpravo zobrazuje body na jedné svislé přímce v reálném prostoru, obrázek vpravo zobrazuje tyto body v parametrovém prostoru. Přestože body náleží jedné přímce, v parametrovém prostoru nemají jednotlivé přímky průsečík. Vytvořeno pomocí nástroje Desmos \cite{desmos}.}
    \label{fig:hough_transform_vertical_line_problem}
\end{figure}

Tento nedostatek Houghovy transformace můžeme vyřešit tak, že změníme způsob, jakým přímku popisujeme: Jedním z řešení je použít pro popis přímky takzvanou \textbf{Hesseho normální formu}: Prvním parametrem $\rho$ udáváme vzdálenost přímky od počátku souřadného systému (jinými slovy délku úsečky, která je ohraničena počátkem souřadného systému a bodem na popisované přímce a která svírá s popisovanou přímkou pravý úhel), druhým parametrem $\theta$ pak popisujeme úhel mezi osou x a úsečkou, spojující počátek souřadného systému s nejbližším bodem popisované přímky \cite{duda_generalized_hough}:

\begin{equation}\label{eq:hough_transform_hesse_normal_form}
    \rho = x \cos \theta + y \sin \theta
\end{equation}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/hough_transform.png}
    \caption{Vizualizace parametrů Houghovy transformace v Hesseho normální formě. Délka úsečky $\rho$ značí nejkratší vzdálenost přímky od počátku souřadnicového systému, $\theta$ je úhel mezi úsečkou $\rho$ a osou x (viz rovnice \ref{eq:hough_transform_hesse_normal_form}). Převzato z \cite{duda_generalized_hough}.}
    \label{fig:hesse_hough_transform}
\end{figure}

Uveďme opět příklad se svislou přímkou z obrázku \ref{fig:hough_transform_vertical_line_problem}: Vzhledem k tomu, že používáme jinou transformaci, body se již nebudou do parametrového prostoru promítat jako přímky, nýbrž jako sinusovky. V bodě, kde se všechny sinusovky protnou, můžeme analogicky jako u předchozích příkladů vyčíst parametry přímky, na které se body nachází (viz. obrázek \ref{fig:hough_transform_vertical_line_solution}).

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/hough_vertical_normal_parameter.png}
    \caption{Parametrový prostor Houghovy transformace bodů z obrázku \ref{fig:hough_transform_vertical_line_problem}. Vytvořeno pomocí nástroje Desmos \cite{desmos}.}
    \label{fig:hough_transform_vertical_line_solution}
\end{figure}

Všimněme si u transformace v normální formě, že se průsečíky jednotlivých sinusovek periodicky opakují: Tento jev pozorujeme z toho důvodu, že přímka se pomocí Hesseho normální formy dá vyjádřit vícero způsoby, což vyplývá z faktu, že úsečka $\rho$ může s osou $x$ svírat úhly ${\theta + 2k\pi}$ radiánů (první a třetí průsečík v obrázku \ref{fig:hough_transform_vertical_line_solution}), popřípadě ${\theta + (2k+1)\pi}$ radiánů (druhý průsečík, v takovém případě nám ale délka úsečky $\rho$ vyjde \uv{záporná}).

Generalizovaná forma Houghovy transformace, jak píše Duda a kolektiv \cite{duda_generalized_hough}, může krom přímek v obraze hledat jakýkoli geometrický útvar, který lze popsat dvěma parametry: pro nás může být relevantní například hledání kružnic: mějme například k dispozici čtyři body, které náleží jedné kružnici: ${A_4 = (1,2)}, {B_4 = (2,1)}, {C_4 = (2,3)}, {D_4 = (3,2)}$. Rovnice pro popis kružnice má následující tvar:

\begin{equation}\label{eq:circle}
    r^2 = (x-h)^2 + (y-k)^2
\end{equation}

kde $r$ značí poloměr kružnice a bod $(h,k)$ značí její střed. V tomto případě přichází drobná komplikace, protože kružnici definujeme pomocí tří parametrů. Tento problém obvykle vyřešíme tak, že jeden z hledaných parametrů známe. Pro náš příklad řekněme, že víme, že kružnice má poloměr $r=1$: parametrový prostor tak zůstává dvojrozměrný a můžeme tak promítnout jednotlivé body do parametrového prostoru jako kružnice. V bodě, kde se všechny kružnice v parametrovém prostoru protnou, najdeme parametry pro kružnici v reálném prostoru (konkrétně pak polohu středu kružnice).

\begin{figure}
    \centering
    \includegraphics[width=0.3\textwidth]{obrazky-figures/hough_transform_circle.png}\hfill
    \includegraphics[width=0.3\textwidth]{obrazky-figures/hough_transform_circle_parameters.png}\hfill
    \includegraphics[width=0.3\textwidth]{obrazky-figures/hough_transform_circle_solved.png}\hfill
    \caption{Verze generalizované Houghovy transformace, hledající kružnice. Body, náležící jedné kružnici (obrázek vlevo) promítneme do parametrového prostoru (obrázek uprostřed). Jejich průsečík udává střed hledané kružnice (vpravo). Převzato z \cite{hough_watershed_transforms}.}
    \label{fig:hough_transform_circles}
\end{figure}

Houghovu transformaci můžeme využít pro rekonstrukci trajektorií následujícím způsobem: V parametrovém prostoru vytvoříme histogram, kde zaznamenáváme, kolik obrazů bodů v parametrovém prostoru prochází každým bodem. Poté postupujeme globálně po celém prostoru a hledáme body s největším počtem procházejících přímek. Z takových bodů postupně vyčítáme parametry pravděpodobných trajektorií a body v reálném prostoru, které se do daného bodu v parametrovém prostoru promítly, vyjmeme a proces opakujeme. Tímto způsobem zároveň potlačujeme šum, protože šumem vzniklé \uv{falešné} zásahy zpravidla nebudou tvořit v parametrovém prostoru průsečíky. 

Tato metoda má pro naše účely jednu zásadní nevýhodu: Houghova transformace pracuje nad celým stavovým prostorem, což v případě sledování částic u CBM znamená značnou paměťovou náročnost a výrazné ztížení paralelizace a akcelerace na hardwarových platformách. Z tohoto důvodu je použití Houghovy transformace sice možné pro účely akcelerace výpočtu ale nepraktické.

\subsection*{Celulární automat}

\textbf{Celulární automat} (CA\footnote{CA -- Cellular Automaton (Buněčný automat)}) je souhrnné označení pro typicky diskrétní výpočetní model, studovaný v rámci \textbf{teorie automatů}. Jedná se o překvapivě multifunkční jednoduchý nástroj, hojně používaný pro simulace v oblastech dopravy, epidemiologie, fyziky a mnoho dalších.

Celulární automat můžeme zpravidla definovat pomocí těchto čtyř kritérií:

\begin{itemize}
    \item \textbf{pole buněk}, které určuje rozměr celulárního automatu. Pole buněk může mít libovolný počet rozměrů, obvykle jeden nebo dva. Dále rozlišujeme, jestli je pole nekonečné, nebo ohraničené (v případě, že je ohraničené, rozlišujeme, jak jsou definovány stavy hranic pole). V neposlední řadě určujeme tvar buněk v poli (například čtvercové, šestiúhelníkové, krychlové u 3D atd.), což dále ovlivňuje tvar okolí každé buňky,
    \item \textbf{Množina stavů buňky (S)}, automat může počítat s libovolným \textbf{konečným} počtem možných stavů, kterých mohou jednotlivé buňky nabývat (v nejjednodušším případě například dvojice stavů 0,1),
    \item \textbf{okolí buněk} určuje, které buňky (a jejich stavy) bude každá buňka brát v potaz, když bude rozhodovat o svém příštím stavu. Rozlišujeme například Mooreovo okolí (čtvercové okolí 8 buněk, majících s buňkou společnou stranu nebo vrchol), rozšířené Mooreovo okolí (stejné, jako Mooreovo, rozšířené o Mooreovo okolí všech osmi sousedů, dohromady tedy 24 buněk v okolí), nebo například Von-Neumannovo okolí (pouze čtyři buňky, které mají společnou stranu s buňkou patří do jejího okolí),
    \item \textbf{přechodová funkce}: pravidla, která určují, jak na základě stavu a okolí buňky vypočítat příští stav buňky. Všechny změny stavu probíhají zároveň na základě posledních stavů pole buněk. 
\end{itemize}

Nejznámější příklad celulárního automatu je hra \uv{Život} Johna Conwaye. Následující odstavec čerpá ze článku Martina Gardnera a skript Jarkka Kari \cite{game_of_life, cellular_automata}: Hra \uv{Život} je celulární automat s nekonečným 2D polem buněk o dvou možných stavech (0 -- \uv{mrtvá} a 1 -- \uv{živá}) a jednoduchým Mooreovým okolím. Pravidla pro hru \uv{Život} byla Conwayem navržena tak, aby splňovala následující 3 požadavky:

\begin{enumerate}
    \item Neměla by existovat žádná počáteční konfigurace, pro kterou by existoval jednoduchý důkaz, že populace (počet \uv{živých} buněk, nesoucích hodnotu 1) může neomezeně růst.
    \item Měly by existovat počáteční konfigurace, které \textit{zřejmě} způsobují neomezený růst populace.
    \item Měly by existovat počáteční konfigurace, které rostou a vyvíjejí se po delší dobu, než se doberou k jednomu ze tří možných konců: úplně \uv{vymřou} (tj. v určitý čas nezbude žádná \uv{živá} buňka, ať už kvůli příliš vysoké, nebo příliš nízké hustoty \uv{zalidnění}), ustálí se ve stabilním rozložení, které se dále již nezmění, nebo se ustálí v oscilující konfiguraci, během které se s určitou periodou donekonečna opakuje konečné množství konfigurací.
\end{enumerate}

Tyto požadavky mají zařídit, že se samotný celulární automat bude chovat \textbf{nepředvídatelně}. Conway tedy navrhl pravidla vývoje populace, která Martin Gardner popsal jako \uv{rozkošně jednoduchá} \cite{game_of_life}.

Conwayova hra \uv{Život} je příklad tzv. \textbf{totalistického celulárního automatu}, což znamená, že nezáleží na konkrétní poloze \uv{živých} a \uv{mrtvých} buněk v rámci okolí buňky, což snižuje počet pravidel ze $2^9 = 512$ (počet možných konfigurací Mooreova okolí krát dva, protože pravidla jsou různá pro \uv{živou} a \uv{mrtvou} buňku) na $2 \times 9 = 18$ (výběr pravidla je závislý pouze na stavu buňky a počtu \uv{živých} buněk).

Pravidla celulárního automatu tak můžeme definovat následovně:

\begin{enumerate}
    \item Přežití: Pokud je buňka živá a má právě dva nebo tři živé sousedy, přežívá do další evoluce
    \item Smrt: Každá živá buňka, která má čtyři a více živé sousedy, zemře kvůli přelidnění. Taktéž každá živá buňka, která má pouze jednoho nebo žádného živého souseda, zahyne samotou.
    \item Zrození: Každá \uv{neživá} buňka, která má \textbf{právě} tři živé sousedy, se \uv{narodí} a v příští evoluci algoritmu bude živá.
\end{enumerate}

V průběhu let hra Život přilákala velké množství nadšenců, kteří postupně shromáždili rozměrnou knihovnu vzorů s různými chováními, které lze kategorizovat pomocí mezi nadšenci uznávané terminologie \cite{cellular_automata}:

\begin{itemize}
    \item \textit{still life} (česky \textit{stabilní život}): stabilní vzor, ve kterém se nerodí žádné nové buňky a žádné naopak neumírají, vzor tedy v této konfiguraci, pokud není zvenčí narušen zůstane donekonečna. Nejjednodušší představitel tohoto vzoru je \textbf{blok}, čtyři živé buňky, poskládané do $2 \times 2$ čtverce. Další příklad je ukázán na obrázku \ref{fig:game_of_life}(a).
    \item \textit{oscilátor}: Periodická konfigurace. Konfigurace se může mezi jednotlivými generacemi měnit, avšak po $n$ generacích se vrátí do původní konfigurace, říkáme tedy, že vzor má periodu $n$. Můžeme také tvrdit, že stabilní život je oscilátor s periodou 1. Nejmenší příklad oscilátoru je takzvaná \textit{blikačka} (anglicky \textit{blinker}), která se skládá ze tří sousedících buněk v jedné přímce. S každou další evolucí se změní ze svislé na vodorovnou a naopak. Další příklad oscilátoru (takzvaná \textit{žába}) je vyobrazen na obrázku \ref{fig:game_of_life}(b).
    \item \textit{vesmírné lodě}: Konfigurace, která se s určitou periodou objeví znovu, oproti oscilátorům se liší v tom, že se znovu objeví na jiném místě v poli, působí tak, že \uv{cestuje}. Na obrázku \ref{fig:game_of_life}(c) vidíme nejjednodušší vesmírnou loď, takzvaný \textit{kluzák}.
    \item \textit{dělo}: oscilující konfigurace, která se stejně jako oscilátor s určitou periodou vrátí do původní podoby, na rozdíl od oscilátoru ale během této periody vygeneruje vesmírnou loď. Na obrázku \ref{fig:game_of_life}(d) můžeme vidět takzvané \textit{kluzákové dělo}.
\end{itemize}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/game_of_life_patterns.png}
    \caption{Vybrané objekty Conwayovy hry Život: stabilní život (a), oscilátor (b), kluzák (c), kluzákové dělo (d). Převzato z \cite{cellular_automata}.}
    \label{fig:game_of_life}
\end{figure}

Navzdory své jednoduchosti vykazuje hra Život natolik sofistikovanou vlastnost jako \textit{sebeorganizace}: I když celulární automat inicializujeme v chaotické počáteční konfiguraci, s postupem času začne metoda vykazovat organizovanou strukturu v podobě (zpravidla) oscilátorů a stabilních vzorů. Tuto vlastnost můžeme využít pro rekonstrukci částic, zejména pro potlačení šumu. Můžeme se na hru Život (a celulární automaty obecně) dívat jako na lokální verzi \textbf{rekurentních neuronových sítí (RNN)}\footnote{RNN -- Recurrent Neural Network (rekurentní neuronová síť)}. Navíc, díky jejich lokální povaze, jsou celulární automaty vhodné pro paralelní přístup, kdy v extrémním případě můžeme na každou buňku nahlížet jako na samostatné vlákno programu \cite{akishina4D}. 

Celulární automaty jsou v dnešní době hojně využívány pro rekonstrukci trajektorií v mnoha HEP experimentech, mezi nimi ALICE\footnote{ALICE -- A Large Ion Collider Experiment}, HERA-B\footnote{HERA -- Hadron-Elektron-RingAnlage (Hadron-Elektronový prstencový urychlovač)}, K2K\footnote{K2K -- KEK to Kamioka}, LHCb\footnote{LHCb -- Large Hadron Collider beauty experiment}, NEMO\footnote{NEMO -- Neutrino Ettore Majorana Observatory}, STAR\footnote{STAR -- Solenoidal Tracker} nebo ARES\footnote{ARES -- Automotive Research Experiment Station}, který byl prvním experimentem, který používá celulární automaty pro rekonstrukci trajektorií a konfigurace automatu je velmi podobná Conwayově původní hře Život. Každý z experimentů má jiné požadavky na algoritmus rekonstrukce trajektorií, tak i jednotlivé celulární automaty jsou rozdílné. My se budeme dále zabývat celulárním automatem použitým v experimentu CBM.

\subsection*{Celulární automat pro experiment CBM}
Následující podsekce čerpá z práce Dr. Akishiny \cite{akishina4D}. Celulární automat na podobné bázi jako Conwayova hra Život sice funguje dobře pro relativně nízké množství trajektorií, jako je experiment ARES, pro detektory s výrazně vyšším počtem trajektorií jako je CBM budeme ale nuceni zavést velké změny do tohoto algoritmu. V prvé řadě je potřeba výrazně změnit pohled na pojem buňky automatu: Buňka by měla být schopna vyjadřovat co největší množství parametrů možné trajektorie. Tam, kde detektor experimentu ARES poskytoval nejen informace o pozici, ale i směru zásahu, bude strategie rekonstrukce výrazně odlišná od experimentů, používajících pixelové nebo proužkové detektory (jako v případě CBM), kde máme k dispozici pouze informace o polohách zásahů. Protože neobsahuje informace o směru trajektorie, nenese jeden zásah v hladině detektoru dostatečně velké množství informací, aby mohl popsat trajektorii částice, která ho způsobila. Na druhou stranu, povědomí o parametrech trajektorie můžeme zlepšit, když v rámci jedné buňky automatu zkombinujeme data z vícero po sobě jdoucích zásahů.

\begin{figure}
    \centering
    \includegraphics[height=0.7\textheight]{obrazky-figures/cbm_cellular_automaton.png}
    \caption{Zjednodušená ilustrace celulárního automatu, používaného například v experimentu CBM pro rekonstrukci trajektorií částic. Svislé přerušované čáry představují jednotlivé detektorové hladiny, zelené a modré značky značí zásahy způsobené dvěma různými částicemi, bílá značka je falešný zásah způsobený šumem. Segmenty trajektorií jsou zaznamenány jako černé, plné úsečky. Přejato z \cite{akishina4D}.}
    \label{fig:cbm_cellular_automaton}
\end{figure}

Jak je ukázáno na obrázku \ref{fig:cbm_cellular_automaton}, v případě experimentů s vysokou hustotou trajektorií definujeme buňku automatu odlišně od intuitivního přístupu: Zde je každá buňka potenciální segment trajektorie, skládající se ze dvou sousedních zásahů v detektoru. Vstupní informace algoritmu jsou vyčtené zásahy v jednotlivých detektorových hladinách: Zeleně vybarvené značky reprezentují zásahy způsobené jednou částicí, modré značky způsobila částice jiná. Bílá značka značí falešný zásah způsobený šumem. V první fázi algoritmus spojí všechny zásahy na sousedních detektorových hladinách a takto vytvořené segmenty projde a zahodí takové kombinace, které by nikdy nemohly patřit reálné částici. Od této chvíle už algoritmus nikdy nebude pracovat s jednotlivými zásahy, ale vždy s buňkami v podobě dvojic zásahů na sousedních detektorových hladinách.

Jakmile jsou \uv{buňky} vytvořeny, může začít evoluce celulárního automatu: Metoda hledá pro každou buňku živé sousedy (sousedy v tomto případě definujeme tak, že mají s buňkou společný zásah), kteří mají přibližně stejný směr, přičemž bere v potaz fyzikální model částic (zde nám pomůže výše zmíněný Kalmanův Filtr), protože takoví sousedé mají vyšší pravděpodobnost, že patří ke stejné trajektorii. Zároveň je každé buňce postupně přiřazeno číslo, které reprezentuje počet sousedů vlevo (čímž ohodnocujeme délku tohoto kandidáta na trajektorii, čím delší, tím lepší). Tímto způsobem získáváme stromovou strukturu kandidátů na trajektorie. V tomto okamžiku dochází k soutěži jednotlivých kandidátů: Přežít mohou pouze ti nejdelší kandidáti, u kterých Kalmanův Filtr spočítal nejlepší hodnotu $\chi^2$ a kteří nesdílí žádné společné zásahy s lepšími kandidáty \cite{akishina4D}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kapitola: Návrh a implementace simulátorů částic v detektoru CBM %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Návrh a implementace simulátorů částic v detektoru CBM}\label{implementation}

Tato kapitola se zabývá podrobným rozborem vlastního algoritmu generování trajektorií částic. Podrobně se tato kapitola bude zabývat návrhem, implementací, verifikací a optimalizací dvou simulátorů drah částic za účelem získání testovacích dat pro účely návrhu a verifikace výsledného algoritmu. Zvláštní pozornost je pak věnována demonstraci vlastností jednotlivých simulačních modelů a porovnání jejich praktikality pro generování částic tak, aby vygenerované částice mohly být použity pro verifikaci a měření výkonu sledovacího algoritmu (a jeho součástí). Závěrem se tato kapitola také zabývá možnými problémy jednotlivých modelů, které generování vhodných částic komplikují.

%Na základě získaných znalostí byl vytvořen simulační model, jehož účel bylo generovat data co možná nejpodobnější těm, jaká by generoval funkční detektor. Takový model využijeme k naprogramování simulátoru, který na základě zákonů pro pohyb nabitých částic v homogenním elektromagnetickém poli a zadané konfigurace detektoru vypočítá měření simulovaného detektoru. Při konfiguraci detektoru budeme uvažovat proměnné jako frekvenci srážkových událostí, rychlost a typ částic v urychlovači a v neposlední řadě množství, rozložení, velikost a rozlišení drátových komor detektoru.


Algoritmus sledování částic, diskutovaný v této práci, využívá do značné míry technologie Kalmanova filtru (respektive jeho rozšířené varianty, \textbf{EKF}). Implementace vhodného Kalmanova filtru je komplikovaný proces s mnoha faktory, které mohou drasticky ovlivnit kvalitu jeho výsledku, je proto velmi důležité správně porozumět jeho fungování. Z tohoto důvodu tato práce uvádí za účelem lepšího porozumění algoritmu příklad, který by měl napomoci snazšímu pochopení funkce, kterou Kalmanův filtr vykonává při rekonstrukci trajektorií částic v detektoru.

\section{Simulátor částic v detektoru}\label{simulator}

Simulační model, který byl použit pro generování trajektorií částic v detektoru, je založen na Lorentzových zákonech pohybu částic v elektromagnetickém poli. Problém generování trajektorií částic v detektoru je složitější, než v případě jiných simulačních modelů z toho důvodu, že (zvláště v případě nehomogenního magnetického a elektrického pole) takový model je nepraktické řešit analyticky. Z tohoto důvodu je výhodnější generovat trajektorie částic pomocí spojité (respektive kombinované) simulace.

Pro účely generování testovacích dat se lze dopustit některých zobecnění, která výrazně usnadní jak vlastní generování částic, tak vlastní práci s daty. Konkrétně zjednodušíme fyzikální model následovně:

\begin{itemize}
    \item \textbf{homogenní magnetické pole}: Simulátor si klade za cíl modelovat pouze trajektorii v rámci prvního stanoviště detektoru (STS), který se celý nachází uprostřed velkého supravodivého magnetu. Pro experiment CBM se v současné době používá 3D mapa nehomogenního magnetického pole o velikosti téměř 90MB. Pro naše účely je ale taková implementace zbytečně složitá a vzhledem k tomu, že jsou všechny stanice detektoru STS, který využíváme pro rekonstrukci trajektorií, plně uvnitř tohoto supravodivého magnetu, můžeme (zvláště u prvních detektorových hladin, kde se přesnost magnetického pole projeví nejvíce a které jsou nejblíže středu magnetu, kde je magnetické pole \uv{nejhomogennější}) minimálně uvnitř magnetu pro naše účely předpokládat, že magnetické pole je homogenní. V případě rekonstrukce trajektorií v praxi ale mapa nehomogenního magnetického pole pomůže k dalšímu zpřesnění modelu pro Kalmanův Filtr, díky čemuž je snažší eliminovat šum a celkově zvýšit efektivitu filtrace.
    
    \item \textbf{nedochází k interakcím mezi částicemi}: Ve chvíli, kdy částice vstoupí do detektoru, předpokládáme, že interagují pouze s detektorem, nikoli mezi sebou. Díky tomuto předpokladu si můžeme dovolit simulovat každou částici zvlášť, a v případě potřeby můžeme generování částic paralelizovat, aniž by bylo nutné jednotlivé procesy simulátoru synchronizovat.
    
    \item \textbf{všechny částice jsou primární s vysokou energií}: U reálného detektoru předpokládáme výskyt čtyř typů trajektorií, které vyplývají z vlastností a konfigurace experimentu: Částice se liší v energii, což ovlivňuje jejich zakřivení v detektoru a zároveň míru mnohonásobného rozpadu při průchodu detektorem, což výrazně ovlivňuje nejistotu systému. Dále se částice mohou lišit časem a místem vzniku a poločasem rozpadu, kdy částice mohou vzniknout v místě kolize s terčem (pak říkáme, že je částice \textit{primární}), nebo vzniknout v důsledku rozpadu jiné částice uvnitř detektoru (takové částici říkáme sekundární). Nakonec z nedokonalosti detektoru lze předpokládat, že některé průlety částice detektorovou hladinou nejsou zaznamenány. Pro naše účely budeme předpokládat, že všechny částice v detektoru vznikly v místě terče a že se pohybují s vysokou energií. Nepředpokládáme tedy, že některé částice vznikly během průletu detektorem, protože vytvořit takový model by vyžadovalo výrazně hlubší porozumění problematice kvantové chromodynamiky a kvantové fyziky obecně, nehledě na fakt, že se jedná o \textit{cutting edge} odvětví, které je momentálně intenzivně zkoumáno a výzkum takovýchto modelů je jedna z motivací experimentu CBM.
    
    \item \textbf{všechny částice mají dlouhou střední délku života}: Předpokládáme, že všechny částice, které v rámci detektoru simulujeme, mají dostatečně dlouhou střední délku života, aby stihly ze sledovaného prostoru detektoru vyletět. Toto zjednodušení zavádíme kvůli tomu, aby bylo chování částic konzistentní s ostatními zjednodušeními modelu: Pokud bychom připustili, že některé částice mají kratší délku života, než je potřebné na opuštění prostoru stanoviště STS, a museli bychom v tom případě řešit problém, kdy uvnitř detektoru dochází k sekundárním rozpadům částic, které generují částice nové, sekundární. Tohle chování by bylo ale v rozporu s předchozím předpokladem, proto ho pro naše účely zanedbáváme.
    
    \item \textbf{neuvažujeme mnohonásobný rozpad}: Pro naše účely nebudeme brát v úvahu mnohonásobný rozpad částic vzhledem k tomu, že bychom museli na místě každé detektorové hladiny přepočítávat nejen trajektorii částice, ale i její energii, což je netriviální úkon, který sice zpřesní simulaci částic, toto zpřesnění ale pro nás během implementace sledovacího algoritmu nemá velký přínos. Zanedbání mnohonásobného rozpadu má pro účely simulace a rekonstrukce trajektorií dva výrazné dopady: Zaprvé zjednoduší simulaci průchodu detektorem, kdy odpadá nutnost výpočtu mnohonásobného rozpadu v každé detektorové hladině (při průchodu každou hladinou detektoru STS dochází v závislosti na úhlu dopadu kolem 24000 srážek, které ovlivňují výslednou trajektorii částice a její energii), zároveň zjednoduší proces Kalmanovy Filtrace tím, že můžeme fyzikální model definovat přesně a procesní šum modelu je díky tomu zanedbatelný. Zde je třeba zmínit, že i když sestrojení matice procesního šumu není příliš komplikovaný proces, samotný výpočet procesního šumu nemá principiální vliv na rekonstrukční algoritmus ani jeho složitost.
    
    \item \textbf{nulové elektrické pole}: Protože se částice pohybují relativistickými rychlostmi, je výpočet jejího chování výrazně složitější. Přítomnost nenulového elektrického pole by způsobilo neúměrné ztížení práce se simulátorem, přičemž jeho zohlednění přinese pouze minimální zpřesnění modelu (magnetické pole je v tomto případě výrazně dominantní silou). Díky tomuto předpokladu získáme velkou výhodu z toho důvodu, že na rozdíl od magnetického pole elektrické pole ovlivňuje množství energie částice, což v důsledku zákonů speciální relativity způsobuje, že se mj. hmotnost částice stává proměnnou. V opačném případě je možné vlivy speciální relativity vypočítat již během inicializace simulace a během simulace \textit{de facto} pracujeme s klasickou mechanikou.    
\end{itemize}
 
Všechna tato zobecnění nám umožní vytvořit vhodný simulátor, který bude generovat částice vhodné pro vývoj a validaci sledovacího algoritmu. Takový přístup má nevýhodu v tom, že takto vytvořená testovací data nejsou vhodná pro hodnocení výkonu výsledného algoritmu, vývoj takového simulátoru by ale byl nad rámec této práce, zjednodušený simulátor je proto pro tuto práci vyhovující.

\section{Metody řešení diferenciálních rovnic}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/particle_spiral_in_detector.png}
    \caption{Trajektorie pomalé, kladně nabité teoretické částice v detektoru s homogenním magnetickým polem ve směru osy $Y$. Červeně jsou vyznačeny průniky částice s hladinami detektoru, modře pak vlastní trajektorie částice. Trajektorie byla vytvořena pomocí simulačního modelu a výsledný graf byl vytvořen v programu MATLAB.}
    \label{fig:particle_detector_spiral}
\end{figure}

Abychom mohli úspěšně simulovat průlet částic detektorem, je třeba správně implementovat zvolený fyzikální model, definovaný soustavou diferenciálních rovnic vyšších řádů. Proto je nutné správně zvolit metodu řešení diferenciálních rovnic, na které lze potom celý simulační model postavit. Pro účely této práce bylo zvoleno numerické řešení, za zmínku ale stojí i možnost \textbf{analytických řešení} trajektorie, přestože pro naše účely nejsou zcela vhodné (simulační model je použitý i v rámci vlastního sledovacího algoritmu, kde již s homogenním magnetickým polem nelze počítat, numerické řešení je zde vhodnější): V případě, že se nabitá částice pohybuje ve směru kolmém na směr magnetického pole, její trajektorie pak nabývá tvaru \textbf{kružnice}. V případě, že se směr částice kolmý není (tj. její má nenulovou složku $\dot{y}$, popřípadě $\theta_y$), trajektorie nabývá tvaru \textbf{šroubovice}, kterou můžeme vyjádřit jako \textbf{komplexní exponenciálu} s osou ve směru magnetického pole (v našem případě ve směru osy y):

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/complex_exponential.png}
    \caption{Komplexní exponenciála a její projekce do reálné roviny (kosinusovka), imaginární roviny (sinusovka) a komplexní roviny (kružnice). Převzato z \cite{complex_exponential}.}
    \label{fig:complex_exponential}
\end{figure}

\begin{equation}\label{eq:complex_exponential}
    f(y) = R \times e^{i2\pi k y} + w    
\end{equation}

kde $R$ je poloměr otáčení šroubovice, $e$ je Eulerovo číslo, $i$ je imaginární jednotka, $k$ je perioda exponenciály a $w \in \mathbb{I}$ je posunutí osy otáčení exponenciály v komplexní rovině, která je kolmá na směr magnetického pole (v našem případě rovina $xz$). $R$ můžeme také definovat jako poloměr oblouku trajektorie, jeho znaménko pak určuje směr zakřivení, který je určen nábojem částice. Vzhledem k tomu, že simulační model předpokládá, že všechny částice mají počátek v terči, posunutí $w$ a poloměr $R$ se stávají lineárně závislými, a tak můžeme nabitou částici v homogenním magnetickém poli vyjádřit pomocí jednoho reálného parametru $k$ a jednoho komplexního parametru $w$, kde parametr $k$ udává úhel $\theta_y$ a parametr $w$ osu otáčení komplexní exponenciály. Toto řešení má ale dvě okrajové situace, kvůli kterým je pro účely tohoto simulátoru nevhodné:

\begin{itemize}
    \item \textbf{částice s nulovým nábojem}: v případě, že takto modelovaná částice má nulový náboj, vzniká problém s parametrem $w$, který v tomto případě nemá řešení (u přímé trajektorie nelze určit poloměr, proto nelze určit ani bod na komplexní rovině, kterým by procházela osa komplexní exponenciály).
    \item \textbf{částice s nulovým úhlem $\theta_y$}: částice, jejíž trajektorie je rovnoběžná s rovinou $xz$ taktéž nemá v tomto vyjádření řešení, protože taková komplexní exponenciála by měla teoreticky nulovou periodu, narážíme tedy na analogický problém, jako při Houghově transformaci svislé přímky (viz obrázek \ref{fig:hough_transform_vertical_line_problem}).
\end{itemize}

Z tohoto důvodu se tato práce nebude možnostmi analytického řešení trajektorií v detektoru dále zabývat. 

\subsection*{Numerické metody}\label{numeric methods}
Po vyloučení použití analytických metod je třeba přistoupit k použití numerických metod. V úvahu jsou brány zejména dvě jednokrokové metody: \textbf{Eulerova metoda} a \textbf{metoda Runge-Kutta} čtvrtého řádu (dále RK4):

\textbf{Eulerova metoda} je nejjednodušší metodou numerického řešení obyčejných diferenciálních rovnic s danými počátečními podmínkami. Mějme obyčejnou diferenciální rovnici prvního řádu:

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/euler_method.png}
    \caption{Graf fungování jednoduché Eulerovy metody. Zeleně je vyznačena reálná funkce, modře jednotlivé kroky Eulerovy metody. Obrázek byl převzat ze skript pro předmět Numerická matematika a pravděpodobnost na Fakultě Informačních Technologií VUT \cite{INM}.}
    \label{fig:rk4}
\end{figure}

\begin{equation}\label{eq:numeric_method_benchmark}
    y'(t) = f\left(t, y(t)\right), y(t_0) = y_0 
\end{equation}

jejíž hodnotu chceme numericky spočítat. Zvolíme délku kroku numerické metody $h$ (pro každé $t_n$ platí $t_n = t_0 + nh$) a pomocí rekurentního vzorce Eulerovy metody začneme s numerickou integrací:

\begin{equation}\label{eq:euler_recurrent_equation}
    \begin{split}
        y_{n+1} &= y_n + hf\left(t_n, y_n\right) \\
        t_{n+1} &= t_n + h
    \end{split}
\end{equation}

přičemž platí, že počáteční hodnota $y_0$ v čase $t_0$ je zadána. Myšlenka Eulerovy metody je, že i když neznáme podobu křivky, kterou zkoumáme, známe počáteční bod, ze kterého vycházíme a z diferenciální rovnice dokážeme vypočítat sklon funkce v tomto počátečním bodě, čímž získáme tečnu, po které se posuneme o krátký časový úsek (krok) $h$.

Alternativou Eulerovy metody jsou Runge-Kuttovy metody. Samotná Eulerova metoda odpovídá Runge-Kuttově metodě prvního řádu, pro použití v simulátoru částic byla ale použita nejznámější z Runge-Kuttových metod, konkrétně metoda Runge-Kutty čtvrtého řádu:

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{obrazky-figures/RK4.png}
    \caption{Graf fungování Runge-Kuttovy metody čtvrtého řádu. Modře je vyobrazená reálná funkce, zeleně je vyobrazena výsledná hodnota $y_{n+1}$, červeně jsou zobrazeny jednotlivé koeficienty. Obrázek byl převzat z dokumentace metody RK4 na stránkách \textit{Read the Docs} \cite{rk4_picture}.}
    \label{fig:rk4}
\end{figure}

\begin{equation}\label{eq:rk4}
    \begin{split}
        t_{n+1} &= t_n + h,\\
        y_{n+1} &= y_n + \frac{1}{6} \left( k_1 + 2k_2 + 2_3 + k_4 \right) h, \\
        k_1 &= f\left(t_n, y_n\right), \\
        k_2 &= f\left(t_n + \frac{h}{2}, y_n + h\frac{k_1}{2}\right), \\
        k_3 &= f\left(t_n + \frac{h}{2}, y_n + h \frac{k_2}{2}\right), \\
        k_4 &= f\left( t_n + h, y_n + hk_3\right),
    \end{split}
\end{equation}

Kde $k_1,\dots,k_4$ jsou koeficienty, definující:

\begin{itemize}
    \item $k_1$ je sklon funkce v bodě $y_n$ (Eulerova metoda),
    \item $k_2$ je sklon funkce ve středu intervalu, v bodě získaným pomocí koeficientu $k_1$ (metoda \textbf{Midpoint}),
    \item $k_3$ je opět sklon funkce ve středu intervalu, tentokrát bod získává na základě koeficientu $k_2$,
    \item $k_4$ je sklon na konci intervalu, k jeho výpočtu je potřeba koeficient $k_3$
\end{itemize}

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{obrazky-figures/Runge-kutta_vs_euler.png}
    \caption{Srovnání propagace chyby numerických metod Runge-Kutta 4. řádu (červeně) a Eulerovy metody (modře). Data byla vygenerována simulátorem částic a zpracována do grafu programem MATLAB.}
    \label{fig:rk4_vs_euler}
\end{figure}

\subsection*{Stabilita numerických metod}
Vzhledem k povaze simulačního modelu, kdy se nabitá částice ve směru kolmém na magnetické pole (za předpokladu, že se nemění kinetická energie částice působením elektrického pole) pohybuje po kružnici (respektive po šroubovici za předpokladu, že složka hybnosti částice ve směru magnetického pole je nenulová), je třeba se pozastavit nad problematikou stability vybraných numerických metod. Bohužel v případě nabitých částic bude jak Eulerova metoda, tak metoda RK4 vždy nestabilní (chování Eulerovy metody a RK4 viz obrázek \ref{fig:rk4_vs_euler}). Metoda $n$-tého řádu je schopna přesně určit řešení diferenciální rovnice za předpokladu, že přesné řešení diferenciální rovnice je polynom řádu maximálně $n$. Eulerova metoda tedy přesně určí řešení, pokud je přesné řešení přímka, u RK4 polynom maximálně čtvrtého řádu atd. Vzhledem k tomu, že pohyb po kružnici (respektive po šroubovici) se jako polynom zapsat nedá, lze přesně spočítat pouze tu složku pohybu, která je rovnoběžná se směrem magnetického pole.

V důsledku tohoto faktu je důležité, aby byla vhodně zvolena délka kroku metody, která sice nezaručí stabilitu, ale poskytne vhodný kompromis mezi rychlostí výpočtu a přesností metody. V tomto ohledu je metoda Runge-Kutta výrazně spolehlivější než Eulerova metoda, přitom poskytuje dobrý kompromis mezi přesností a výpočetní náročností: RK4 je nejvyšší řád Runge-Kuttových metod, kde řád metody odpovídá počtu jejích kroků (RK1/Eulerova metoda: 1 krok, \dots RK4: 4 kroky, ale již RK5: 6 kroků), metody vyšších řádů pak rychle nabývají na výpočetní náročnosti a pro nás jsou tak méně vhodné.

\subsection*{Použití numerických metod pro rovnice vyšších řádů}

Přestože jsou Eulerova metoda i metoda Runge-Kutta určeny pro řešení diferenciálních rovnic prvního řádu, můžeme ji použít i pro obyčejné diferenciální rovnice vyšších řádů tak, že převedeme diferenciální rovnici vyššího řádu na soustavu rovnic prvního řádu. Mějme diferenciální rovnici $N$-tého řádu:

\begin{equation}\label{eq:higher_order_equations}
    y^{(N)}(t) = f\left( t, y(t), y'(t), \dots, y^{(N-1)}(t) \right)
\end{equation}

zavedeme pomocné proměnné $z$ následujícím způsobem:

\begin{equation}\label{eq:higher_order_equations_substitutions}
    \begin{split}
        z_1(t) &= y(t) \\
        z_2(t) &= y'(t) \\
        &\enspace \vdots \\
        z_N(t) &= y^{(N-1)}(t)
    \end{split}
\end{equation}

Obyčejnou diferenciální rovnici vyššího řádu pak můžeme přepsat jako následující soustavu diferenciálních rovnic prvního řádu:

\begin{equation}\label{eq:higher_order_ode_transformation}
    \boldsymbol{z}'(t) = \begin{bmatrix}
        z_1'(t) \\
        z_2'(t) \\
        \vdots \\
        z_{N-1}'(t) \\
        z_N'(t)
    \end{bmatrix} = \begin{bmatrix}
        z_2(t) \\
        z_3(t) \\
        \vdots \\
        z_N(t) \\
        f\left( t, z_1(t), \dots , z_N(t) \right)
    \end{bmatrix} = \begin{bmatrix}
        y'(t) \\
        y''(t) \\
        \vdots \\
        y^{(N-1)}(t) \\
        y^{(N)}(t) 
    \end{bmatrix}
\end{equation}

\section{Pohybové rovnice simulačního modelu, první simulační model}

První verze simulátoru dráhy částic byl navržen tak, aby co nejpřesněji modeloval dráhu teoretických, \uv{pomalých} částic. Pro popis dráhy pomalé částice v elektromagnetickém poli slouží jev zvaný \textbf{Lorentzova síla}, která popisuje sílu působící na nabitou částici v elektromagnetickém poli.

\begin{equation}\label{eq:lorentz_force_classical_vector}
    \boldsymbol{F} = q\left(\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} \right)
\end{equation}

kde $\boldsymbol{E}$ je vektor elektrického pole, $\boldsymbol{v}$ je vektor rychlosti částice v elektromagnetickém poli a $\boldsymbol{B}$ je vektor magnetického pole. Jak bylo zmíněno výše, v rámci simulátoru je elektrické pole $\boldsymbol{E}$ zanedbáno, lze tedy rovnici přepsat pro účely simulátoru následujícím způsobem:

\begin{equation}\label{eq:lorentz_force_vector_no_electrical_field}
    \boldsymbol{F} = q\boldsymbol{v} \times \boldsymbol{B}
\end{equation}

Všimněme si, že v tuto chvíli v rovnici nefiguruje hmotnost částice, protože velikost Lorentzovy síly je na ní nezávislá. Hmotnost částice ale potřebujeme znát pro výpočet zrychlení částice, aplikujeme tedy druhý Newtonův zákon $F = ma$, díky čemuž dostáváme rovnici:

\begin{equation}\label{eq:lorentz_force_newton2}
    \boldsymbol{\dot{v}} = \frac{q}{m} \boldsymbol{v} \times \boldsymbol{B}
\end{equation}

Kde $\boldsymbol{\dot{v}} = \boldsymbol{a}$ je vektor zrychlení částice a $m$ je její hmotnost. V této rovnici lze pro potřeby simulátoru rozepsat vektory následujícím způsobem:

\begin{equation}\label{eq:lorentz_force_classical}
    \begin{split}
        \ddot{x} &= \frac{q}{m} \left(\dot{y} B_z - \dot{z} B_y \right) \\
        \ddot{y} &= \frac{q}{m} \left(\dot{z} B_x - \dot{x} B_z\right) \\    
        \ddot{z} &= \frac{q}{m} \left(\dot{x} B_y - \dot{y} B_x \right)
    \end{split}
\end{equation}

 \begin{figure}
     \centering
     \includegraphics[width=0.45\textwidth]{obrazky-figures/nonrel_sim_demo.png}
     \caption{Demonstrace simulace pohybu různých, nerelativistických částic. Červeně je vyznačena dráha těžké / slabě nabité částice s kladným nábojem, modře je vyznačena lehká / silně nabitá částice se záporným nábojem a zeleně je vyznačená dráha částice s nulovým nábojem. Směr magnetického pole je rovnoběžný s osou $y$, všechny částice mají shodnou počáteční rychlost $\boldsymbol{v} = (0,0,1)\,ms^{-1}$. Graf byl vytvořen za pomoci nástroje Desmos \cite{desmos}.}
     \label{fig:simulator_demo}
 \end{figure}
 
kde $\dot{x}$ je rychlost částice ve směru $x$, $\ddot{x} = \dot{v}_x = a_x$ značí zrychlení ve směru dané osy, $B_x$ je velikost intenzita magnetického pole ve směru $x$ (analogicky pro směry $y$ a $z$). Tímto způsobem jsme získali soustavu tří obyčejných diferenciálních rovnic, které nám poslouží jako první simulační model. Za pozornost stojí ještě fakt, že v případě homogenního magnetického pole, které má směr rovnoběžný s některou z os je možné model dále zjednodušit: Uvažujme případ detektoru, ve kterém má magnetické pole směr rovnoběžný s osou $y$, model lze pak zjednodušit následujícím způsobem:

\begin{equation}\label{eq:lorentz_force_y_axis}
    \begin{split}
        \ddot{x} &= -\frac{q}{m} \dot{z} B \\    
        \ddot{z} &= \frac{q}{m} \dot{x} B
    \end{split}
\end{equation}
 
Kde $B = B_y$ je magnetické pole. Všimněme si, že rovnici pro výpočet $\ddot{y}$ můžeme zcela zanedbat, zrychlení v tomto směru se vždy vyhodnotí jako nulové, v tomto směru se tedy částice budou pohybovat rovnoměrným přímočarým pohybem. Zjednodušení simulátoru tímto způsobem je sice možné, pro účely simulátoru je ale výhodné mít možnost velikost a orientaci magnetického pole nastavovat, bylo tedy učiněno rozhodnutí implementovat obecnější variantu z rovnice \ref{eq:lorentz_force_classical}.
 
\section{Speciální relativita v pohybových rovnicích, druhý simulační model}

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{obrazky-figures/Lorentz_factor.png}
    \caption{Graf funkce rychlosti a Lorentzova faktoru. Při rychlostech blížících se rychlosti světla Lorentzův faktor stoupá k nekonečnu. Graf byl vytvořen pomocí nástroje Desmos \cite{desmos}.}
    \label{fig:lorentz_factor}
\end{figure}

Simulační model popsaný výše má zásadní nedostatek, který je patrný ve chvíli, kdy pomocí tohoto simulátoru vygenerujeme první částice s úkolem \uv{modelovat} reálné částice (například protony či tauony): Protože používáme pro výpočet zrychlení klasické zákony mechaniky, nastává problém ve chvíli, kdy simulátor začne generovat částice při relativistických rychlostech. Protože model nebere zákony speciální relativity v potaz, simulátor může generovat částice, které v praxi neexistují a vzhledem k tomu, že u částic v detektoru CBM očekáváme rychlosti přes $0.9c$, spoléhat se na klasickou mechaniku je nevhodné. Další z nedostatků je, že simulátoru teoreticky nic nebrání generovat částice s rychlostmi přesahujícími rychlost světla.

Abychom mohli simulační model upravit tak, aby zohledňoval relativistické chování částic, musíme nejprve zavést takzvaný \textbf{Lorentzův faktor}:

\begin{equation}\label{eq:lorentz_factor}
    \gamma \equiv \frac{\text{d}t}{\text{d}\tau} = \frac{c}{\sqrt{c^2 - v^2}} = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}
\end{equation}

Kde $c$ značí rychlost světla ve vakuu, $\tau$ je takzvaný \textbf{vlastní čas}. Pohybové rovnice předchozího modelu (rovnice \ref{eq:lorentz_force_classical}), po korekci pro speciální relativitu mají pak následující tvar:

\begin{equation}\label{eq:simulator_relativity}
    \begin{split}
        \ddot{x} &= \frac{q}{\gamma m} \left(\dot{y} B_z - \dot{z} B_y \right) \\
        \ddot{y} &= \frac{q}{\gamma m} \left(\dot{z} B_x - \dot{x} B_z\right) \\    
        \ddot{z} &= \frac{q}{\gamma m} \left(\dot{x} B_y - \dot{y} B_x \right)
    \end{split}
\end{equation}

Zavedení Lorentzova faktoru do pohybových rovnic znamená, že dosud konstantní hmotnost částice $m$ se stává proměnnou, závislou na rychlosti částice. Rozhodnutí předpokládat v simulátoru homogenní magnetické pole a nulovou intenzitu elektrického pole se zde projeví jako výhodné z důvodu, že bez vlivu elektrického pole velikost rychlosti částice $v$ zůstává konstantní, proto zůstává konstantní i její celková hmotnost, kterou pouze při inicializaci simulace upravíme na základě klidové hmotnosti a Lorentzova faktoru. V případě druhého, pro speciální relativitu upraveného simulátoru již dokážeme efektivně modelovat částice s relativistickou rychlostí.

\section{Kombinovaný simulátor částic s logováním zásahů v detektoru}

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{obrazky-figures/25_particles_corrected.png}
    \caption{ Graf 25 náhodně vygenerovaných částic, vlastnostmi odpovídajících částicím $p^+, \Lambda_c^+, \Sigma^-$ a $\Sigma_b^-$. Modře jsou naznačeny vygenerované trajektorie, červeně jsou vyznačeny zásahy částic v hladinách detektoru. Data byla vygenerována pomocí simulátoru částic a graf byl vytvořen pomocí programu MATLAB.}
    \label{fig:25_generated_particles}
\end{figure}

Generování cvičných částic bylo implementováno pro demonstrační účely ve formátu Jupyter notebook, původně v jazyce Python 3.9.4 (za účelem rychlého prototypování). V průběhu práce byl simulátor přepsán v jazyce Julia 1.6.1. díky podobnostem s jazykem Python a výrazně vyšší rychlostí zpracování \cite{julia}. Simulátor byl také implementován i v jazyce C pro porovnání rychlostí jednotlivých implementací.

Původně byl simulátor navržen pro spojitou simulaci, která probíhá ve vlastním čase každé částice od času $t_0 = 0s$ do času $t_{max} = \frac{1.25}{v}$ sekund, kde $1.25$ je délka části dráhy částice, která je ovlivněna supravodivým magnetem CBM. Mimo oblast magnetu není nutné trajektorii částice simulovat, protože oproti oblasti uvnitř magnetu, kde lze magnetické pole \textit{téměř} považovat za homogenní, je výhodnější předpokládat, že mimo magnet je intenzita magnetického pole nulová než se pokoušet magnetické pole jinak aproximovat. Navíc algoritmus pro sledování trajektorií v detektoru používá ve své první fázi pouze data z detektoru STS, který je celý uvnitř magnetu, tvar trajektorie ve vzdálenosti od terče větší než 1 metr pro nás tedy již není zajímavá.

Spojitý simulátor používá Runge-Kuttovu metodu pro simulaci trajektorie, jednotlivé kroky pak zapisuje do výstupního \textit{.csv} souboru. Pro tento účel je sice současný návrh simulátoru postačující, takový simulátor ale není schopen přesně určit, kde částice prošla hladinou detektoru. Z tohoto důvodu byl simulátor rozšířen na \textbf{kombinovaný simulátor}, díky čemuž jsme schopni obsluhovat \textit{stavové události} (události, které nelze naplánovat, v případě simulátoru částic okamžiky, kdy částice proletí hladinou detektoru). 

V okamžiku, kdy se změní stavová podmínka (v případě simulátoru částic taková změna souřadnice $z$, která ukazuje na průlet hladinou detektoru mezi jednotlivými měřeními), je potřeba najít čas, kdy ke změně stavové podmínky (a tedy ke stavové události) došlo. Toho lze docílit například pomocí metody \textbf{půlení intervalu}. Mějme numerickou metodu s krokem $h$, minimálním krokem $h_{min}$ a stavovou podmínkou: krok metody opakujeme tak dlouho, než stavová podmínka změní hodnotu (tj. dojde ke stavové události), v takový okamžik krok zmenšíme na poloviční a obnovíme stav systému z okamžiku před provedením kroku, který stavovou událost vyvolal. Krok metody a postupné půlení délky kroku provádíme tak dlouho, dokud je délka kroku než $h_{min}$, v tu chvíli víme, že při vyvolání stavové události byl nalezen okamžik jejího vyvolání s přesností $\pm h_{min}$. Existuje varianta tohoto algoritmu, která neobnovuje původní stav systému, ale provádí zmenšený krok \uv{pozpátku}. Tato verze metody sice konverguje rychleji, z programátorského hlediska bylo ale snazší implementovat původní metodu. Tímto způsobem simulátor nachází okamžiky, ve kterých ke stavovým událostem došlo a umožňuje její další zpracování. Stavová podmínka je v tu chvíli potvrzena, vyhodnotí se stavová událost (například zapsání hodnot stavu systému nebo vyhodnocení procesu mnohonásobného rozptylu, který sice tento simulátor zanedbává, při použití ve vlastním sledovacím algoritmu je ale jeho vyhodnocení na místě).

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{obrazky-figures/bisection.jpg}
    \caption{Metoda půlení intervalu: Metoda provede krok a vyhodnotí, jestli se \textit{změnila hodnota stavové podmínky} a v závislosti na hodnotě stavové podmínky buďto provede další krok, nebo podělí délku kroku dvěma a provede krok znovu. Obrázek byl převzat z webu GeeksforGeeks \cite{gfg_bisection}}
    \label{fig:bisection_method}
\end{figure}

Stavové události jsou vnitřně zpracovávány pomocí jedné funkce, která při každém kroku porovná polohu částice před a po provedení kroku numerické metody. Funkce má uložené $z$ souřadnice všech hladin detektoru, a pokud je některá z těchto hladin \textbf{mezi} dvěma hodnotami $z$ souřadnice před a po provedení kroku, vyvolá stavovou událost. Stavové události jsou, podobně jako data o trajektorii částice, ukládány do samostatného \textit{.csv} souboru, který bude po anonymizaci (program do \textit{.csv} souboru zapisuje jednoznačné \textit{id} částice, která zásah způsobila) a zašumění dat možné použít jako testovací data pro vývoj sledovacího algoritmu.

\subsection*{Struktura částice v detektoru}

Každá částice je v simulátoru definována jako struktura, nesoucí data o všech vlastnostech a stupních volnosti, potřebných k její jednoznačné definici:

\begin{itemize}
    \item jednoznačný identifikátor částice, sloužící ke kontrole správnosti rekonstruovaných částic (lze tímto způsobem kontrolovat, jestli zásahy detektoru, které algoritmus spojil jako trajektorii jedné částice, skutečně byly vyvolány stejnou částicí),
    \item počáteční čas simulace částice $t$, který,
    \item koncový čas simulace $t_{max}$, který říká simulátoru, kdy má ukončit práci s danou částicí (ať už díky tomu, že částice opustila detektor, nebo došlo k sekundárnímu rozpadu),
    \item klidovou hmotnost částice $m_0$ na základě které dopočítáváme Lorentzův faktor částice pro úpravu vlastností při relativistickém pohybu,
    \item náboj částice $q$,
    \item stavový vektor částice s informacemi o poloze a rychlosti částice,
    \item Lorentzův faktor $\gamma$,
    \item Konstantu $C_{\gamma qm} = \frac{q}{\gamma m}$, kterou použijeme při výpočtu zrychlení.
\end{itemize}

\subsection*{simulace průletu částice detektorem}

Průlet částice detektorem je simulován tak, že nejprve inicializujeme stavový vektor částice, počáteční a koncový čas simulace, její náboj a klidovou hmotnost. Pro zjednodušení generování nových částic jsou v rámci souboru \texttt{physics.h} definovány základní potřebné konstanty jako $\pi$, rychlost světla $c$ nebo elementární náboj $e$. Dále je v souboru definován katalog klidových hmotností a kvantových nábojů (náboj částice $q$ získáme tak, že její kvantové číslo náboje vynásobíme konstantou elementárního náboje) známých, pozorovatelných částic na základě kterých lze následně dopočítat Lorentzův faktor. Po zapsání hodnot částice je struktura částice předána funkci \texttt{calculate\_gqm()}, která vypočítá a uloží hodnotu Lorentzova faktoru. V tuto chvíli je částice připravena k simulaci. 

Zmiňme ještě funkci \texttt{simulation\_loop()}, která zprostředkovává vlastní proces numerické simulace. Funkce má na starosti logování dat simulace (za tím účelem jako parametry přijímá dva souborové deskriptory, do kterých data zapisuje: první soubor je určen logování trajektorie, druhý soubor ukládá zásahy v hladinách detektoru), kontrolu (a obsluhu) stavových podmínek a vlastní krokování vybrané numerické metody. V zájmu obecnosti implementace má funkce parametr \texttt{num\_method}, ukazatel na funkci, pomocí kterého můžeme při volání funkce \texttt{simulation\_loop()} určit, jakou metodu si přejeme pro simulaci použít (implementovány jsou Eulerova metoda a metoda RK4).

Implementace simulátoru v Julia je velmi podobná implementaci v C, s tím rozdílem, že Julia poskytuje za účelem fyzikálních simulací modul \texttt{PhysicalConstants.jl}, který obsahuje všechny informace o důležitých fyzikálních konstantách, zároveň definuje jednotky těchto konstant a umožňuje při programování v jazyce Julia pracovat s fyzikálními jednotkami, čímž usnadňuje vývoj a verifikaci fyzikálních výpočtů. Pro práci s logovanými daty trajektorií je pak použit modul \texttt{DataFrames.jl}, který umožňuje rychlou a výpočetně nenáročnou práci s těmito daty (DataFrames.jl implementuje sadu funkcí nad datovými rámci podobným způsobem jako v jazyce Python knihovna \texttt{Pandas} s tím, že všechny \textit{element-wise} operace jsou implicitně vektorizovány pro rychlejší zpracování).

\subsection*{Odvození fyzikálního modelu pro Kalmanův Filtr}

Protože Kalmanův filtr, který je nezbytný k fungování sledovacího algoritmu, potřebuje vhodný fyzikální model k tomu, aby správně odhadoval příští stav systému, bylo by výhodné použít stejný model, který již používáme pro generování testovacích dat. Fyzikální model odvodíme tak, jak bychom model odvozovali (a dále budeme odvozovat, v případě dalších modelů) v případě použití rozšířeného Kalmanova Filtru (EKF):

Mějme soustavu devíti rovnic, vycházejících ze soustavy \ref{eq:simulator_relativity}:

\begin{equation}\label{eq:simulator_relativity_full}
    \begin{split}
        x_k &= x_{k-1} + \Delta t \, \dot{x}_{k-1} + \frac{\Delta t^2}{2}\ddot{x}_{k-1}\\ 
        \dot{x}_k &= \dot{x}_{k-1} + \Delta t \, \ddot{x}_{k-1} \\
        \ddot{x}_k &= \frac{q}{\gamma m} \left(\dot{y}_{k-1} B_z - \dot{z}_{k-1} B_y \right) \\
        y_k &= y_{k-1} + \Delta t \, \dot{y}_{k-1} + \frac{\Delta t^2}{2}\ddot{y}_{k-1}\\
        \dot{y}_k &= \dot{y}_{k-1} + \Delta t \, \ddot{y}_{k-1} \\
        \ddot{y}_k &= \frac{q}{\gamma m} \left(\dot{z}_{k-1} B_x - \dot{x}_{k-1} B_z\right) \\    
        z_k &= z_{k-1} + \Delta t \, \dot{z}_{k-1} + \frac{\Delta t^2}{2}\ddot{z}_{k-1}\\ 
        \dot{z}_k &= \dot{z}_{k-1} + \Delta t \, \ddot{z}_{k-1}\\
        \ddot{z}_k &= \frac{q}{\gamma m} \left(\dot{x}_{k-1} B_y - \dot{y}_{k-1} B_x \right)
    \end{split}
\end{equation}

Nyní můžeme definovat vektorově ohodnocenou funkci $\textbf{f}_k(\boldsymbol{r}_k^+): \mathbb{R}^9 \to \mathbb{R}^9$, která obsahuje pravé strany rovnice \ref{eq:simulator_relativity_full}:

\begin{equation}\label{eq:kalman_model_function_f}
    \boldsymbol{r}_k^{-} = \textbf{f}_k\left(\boldsymbol{r}_{k-1}^+\right) = \begin{bmatrix}
        x_{k-1} + \Delta t \, \dot{x}_{k-1} + 0.5\Delta t^2 \ddot{x}_{k-1} \\ 
        \dot{x}_{k-1} + \Delta t \, \ddot{x}_{k-1} \\
        q\gamma^{-1} m^{-1} \left(\dot{y}_{k-1} B_z - \dot{z}_{k-1} B_y \right) \\
        y_{k-1} + \Delta t \, \dot{y}_{k-1} + 0.5 \Delta t^2 \ddot{y}_{k-1} \\
        \dot{y}_{k-1} + \Delta t \, \ddot{y}_{k-1} \\
        q\gamma^{-1} m^{-1} \left(\dot{z}_{k-1} B_x - \dot{x}_{k-1} B_z\right) \\    
        z_{k-1} + \Delta t \, \dot{z}_{k-1} + 0.5\Delta t^2\ddot{z}_{k-1} \\ 
        \dot{z}_{k-1} + \Delta t \, \ddot{z}_{k-1} \\
        q\gamma^{-1} m^{-1} \left(\dot{x}_{k-1} B_y - \dot{y}_{k-1} B_x \right)
    \end{bmatrix}
\end{equation}

Protože v tomto modelu předpokládáme, že je magnetické pole homogenní a elektrické pole nulové, fyzikální model je lineární (za povšimnutí stojí i fakt, že v důsledku nulového elektrického pole se nemění celková energie částice, tedy ani její hmotnost, můžeme tedy proměnnou hmotnosti $m$ nahradit konstantou klidové hmotnosti částice $m_0$). Díky tomuto faktu není nutné systém linearizovat a můžeme rovnou sestrojit propagační matici $F_k$:

\begin{equation}\label{eq:model_propagation_matrix_kalman}
    F_k = \begin{bmatrix}
        1 & \Delta t & 0.5 \Delta t^2 & 0 & 0 & 0 & 0 & 0 & 0 \\ % x
        0 & 1 & \Delta t & 0 & 0 & 0 & 0 & 0 & 0 \\              % \dot{x}
        0 & 0 & 0 & 0 & \frac{qB_z}{\gamma m_0} & 0 & 0 & -\frac{qB_y}{\gamma m_0} & 0 \\ % \ddot{x}
        0 & 0 & 0 & 1 & \Delta t & 0.5 \Delta t^2 & 0 & 0 & 0 \\ % y
        0 & 0 & 0 & 0 & 1 & \Delta t & 0 & 0 & 0 \\              % \dot{y}
        0 & -\frac{qB_z}{\gamma m_0} & 0 & 0 & 0 & 0 & 0 & \frac{qB_x}{\gamma m_0} & 0 \\ % \ddot{y}
        0 & 0 & 0 & 0 & 0 & 0 & 1 & \Delta t & 0.5 \Delta t^2 \\ % z
        0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & \Delta t \\              % \dot{z}
        0 & \frac{qB_y}{\gamma m_0} & 0 & 0 & -\frac{qB_x}{\gamma m_0} & 0 & 0 & 0 & 0 \\ % \ddot{z}
    \end{bmatrix}
\end{equation}

Zde je na místě poznamenat, že v případě modelu, který by ve své propagační matici používal stavové proměnné, bylo by již nutné pro použití Kalmanovým Filtrem Matici linearizovat podle vhodně zvoleného $\boldsymbol{r}_k^{lin}$.

Jak bylo zmíněno výše, pro homogenní magnetické pole, které je rovnoběžné s některou z os, je možné model ještě mírně zjednodušit. Mějme homogenní magnetické pole, působící ve směru osy $y$. Pohybová složka ve směru $y$ se tím značně zjednoduší, protože se vždy bude jednat jen o rovnoměrný přímočarý pohyb, což mimo jiné znamená, že můžeme zcela vynechat proměnnou $\ddot{y}_k$ a zmenšit tak velikost stavového vektoru na 8. Výsledná predikční matice pak vypadá následovně:

\begin{equation}\label{eq:model_propagation_matrix_kalman_corrected}
    F_k = \begin{bmatrix}
        1 & \Delta t & 0.5 \Delta t^2 & 0 & 0 & 0 & 0 & 0 \\ % x
        0 & 1 & \Delta t & 0 & 0 & 0 & 0 & 0 \\              % \dot{x}
        0 & 0 & 0 & 0 & 0 & 0 & -\frac{qB_y}{\gamma m_0} & 0 \\ % \ddot{x}
        0 & 0 & 0 & 1 & \Delta t & 0 & 0 & 0 \\ % y
        0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\              % \dot{y}
        0 & 0 & 0 & 0 & 0 & 1 & \Delta t & 0.5 \Delta t^2 \\ % z
        0 & 0 & 0 & 0 & 0 & 0 & 1 & \Delta t \\              % \dot{z}
        0 & \frac{qB_y}{\gamma m_0} & 0 & 0 & 0 & 0 & 0 & 0 \\ % \ddot{z}
    \end{bmatrix}
\end{equation}

\subsection*{Nevýhody simulačních modelů s časovým rozměrem}

Simulátor částic je sice funkční a spolehlivě modeluje i speciální teorii relativity, při pokusu o vytvoření Kalmanova filtru pro použití v rekonstrukčním algoritmu ale vyvstává další problém: Ačkoliv pro naše účely je v pořádku předpokládat, že magnetické pole je konstantní a dokonale rovnoběžné s osou $y$, v reálném detektoru tento předpoklad zavést nelze. Situaci dále komplikuje fakt, že reálné magnetické pole je natolik komplexní, že se nevyplatí ho popisovat analyticky. To je v současné době řešeno tak, že rekonstrukční program drží v paměti velmi detailní 3D mapu směru a intenzity magnetického pole, se kterým potom pracuje. Z tohoto důvodu v reálném rekonstrukčním programu není pro predikci příštího stavu použit prostý výpočet $\boldsymbol{r}_k^-=\textbf{f}_k\left(\boldsymbol{r}_{k-1}^+\right)$, ale následující odhad je získán pomocí numerické simulace od detektorové hladiny $a$ k detektorové hladině $b$. Simulátor částic tedy využijeme i v rámci vlastního rekonstrukčního algoritmu. Simulátor ale vzhledem k tomu, že pracuje v čase, nedokáže předem určit, kdy simulaci zastavit a pro získání příští predikce musí čekat na stavovou událost. Kontrola stavových podmínek je výpočetně zbytečně náročná (a v případě devítisložkového stavového vektoru její náročnost dále roste), proto pokud to není nutné, je žádoucí se v rámci vlastního rekonstrukčního programu kombinovaným simulátorům vyhnout. 

Z tohoto důvodu je třeba zavést třetí simulační model, tentokrát postavený tak, aby nesimuloval pohyb částice v závislosti na \textbf{čase}, ale na \textbf{souřadnici z}.

\section{Pohyb částice v závislosti na vzdálenosti od terče detektoru, třetí simulační model}

V reakci na problémy s implementací simulátoru zmíněné výše byl na základě práce Dr. Akishiny \cite{akishina4D} zaveden poslední fyzikální model a implementován simulátor, tentokrát zásadně odlišný od předchozích verzí. Simulátor vychází z faktu, že částici v magnetickém poli lze jednoznačně popsat pomocí pěti proměnných (říkáme, že stavový prostor má pět stupňů volnosti):

\begin{equation}\label{eq:z_state_vector}
    \boldsymbol{r} \left\{ x, y, t_x, t_y, q/p \right\}
\end{equation}

Kde $x$ a $y$ jsou souřadnice částice pro danou souřadnici $z$, $t_x \equiv \tan \theta_x$ a $t_y \equiv \tan \theta_y$ označují sklon směru částice vůči rovinám ZX ($t_y$) a ZY ($t_x$) a $q/p$ označuje poměr mezi nábojem a hybností částice. Rovnice pro výpočet Lorentzovy síly přepíšeme následujícím způsobem: \cite{akishina4D}

\begin{equation}\label{eq:z_model}
    \begin{split}
        \frac{\text{d}t_x}{\text{d}z} &= c (q/p) \, t_r(t_y(B_z + t_xB_x) - B_y (1 + t_x^2))\\
        \frac{\text{d}t_x}{\text{d}z} &= c (q/p) \, t_r( B_x (1 + t_y^2)- t_x(B_z + t_yB_y))\\
        t_r(z) &= \sqrt{t_x^2+t_y^2 + 1}
    \end{split}
\end{equation}

Kde $c$ je rychlost světla a $B$ je intenzita magnetického pole. Všimněme si, že vzhledem k tomu, že tento model nepoužívá derivace v čase, neznáme rychlost částice, v důsledku čehož je obtížné zjišťovat vliv relativity na dráhu konkrétní částice. Výhodou tohoto přístupu ale je, že i když nevíme, jakou rychlostí částice detektorem proletěla, z podstaty modelu ani není třeba rychlost částice znát, protože částice bude mít stejnou trajektorii nezávisle na tom, jestli je částice těžká a pomalá, nebo lehká a rychlá, protože se bude lišit pouze čas, za který částice detektorem proletí, a ten vzhledem k povaze a vzorkovací frekvenci detektoru nejsme schopni změřit (z pohledu detektoru proletí každá rozpoznatelná částice všemi hladinami detektoru STS ve stejný čas).

Výhoda takového modelu je v tom, že simulátor založený na takovém modelu můžeme spustit v libovolném místě v detektoru a v libovolném místě detektoru můžeme simulaci zase zastavit. Tento fakt nám mimo jiné umožní provádět simulace z bodu $a$ do bodu $b$, což je výhodné pro potřeby rekonstrukčního algoritmu, který nyní může provést pomocí Runge-Kuttovy metody odhad dráhy částice od jedné detektorové hladiny k druhé, ale i pro obecný Kalmanův filtr, jehož predikční matice nyní není závislá na časovém intervalu $\Delta t$, který se mezi hladinami detektoru může lišit, ale na místním intervalu
$\Delta z$, který je v experimentu CBM konstantní (hladiny detektoru jsou rozestavěné po intervalech deseti centimetrů).

Tento přístup má svá vlastní omezení, konkrétně zmíníme dvě:

\begin{itemize}
    \item vzhledem k tomu, že používáme souřadnici $z$ jako časovou složku, některé možné dráhy částic nelze tímto způsobem simulovat. Konkrétně se jedná o částice, které z detektoru vyletí pod pravým úhlem k ose z, a o částice, které mají dostatečně nízkou energii na to, aby se v detektoru \uv{točily} (viz částice na obrázku \ref{fig:particle_detector_spiral}). Intuitivně existují tato omezení z toho důvodu, že v rámci simulace se nemůžeme \uv{vracet v čase}, z matematického hlediska je to pak z toho důvodu, že hodnoty stavových proměnných $t_x$ a $t_y$ mají v případě pravého úhlu k ose $z$ funkční hodnotu limitně v $\pm \, \infty$. Pravděpodobnost výskytu natolik pomalé částice je každopádně velmi nízká. Navíc, v důsledku mnohonásobného rozptylu uvnitř detektorových hladin, jehož vliv je silnější pro částice s nižší energií, je prakticky nemožné, aby byl detektor schopen takovou částici zrekonstruovat.
    
    \item protože neznáme rychlost částice, vliv speciální relativity na částici je \uv{skryt} v proměnné $p/q$. To sice pro potřeby filtrace není problém, v případě, že chceme simulovat dráhu reálné částice, výpočet vlivu relativity na ni je netriviální. Nicméně díky faktu, že model nebere v úvahu čas, není toto omezení nijak závažné.
\end{itemize}

Simulátor byl implementován na velmi podobné bázi, jako první verze simulátoru pro nerelativistické částice, a díky jeho povaze bylo možné upustit od kombinované simulace a simulovat místo toho každý segment trajektorie zvlášť (vzhledem k tomu, že derivujeme podle osy $z$, můžeme částici \uv{odstartovat} v terči a jako $z_{max}$ nastavit první hladinu detektoru ve 30cm, vyčíst hodnoty jeho stavového vektoru v $z_{max} = 30$ a inicializovat simulaci znovu, s hodnotami shodnými, jako na konci simulace, počátečním $z$ v první hladině detektoru a $z_{max}$ nastavíme na druhou hladinu detektoru ve 40 cm atd.), což je nejen výhodné pro vlastní rekonstrukční algoritmus, ale i výpočetně efektivní pro vlastní simulaci částic.

\subsection*{Problém nalezení počátečního stavového vektoru třetího simulátoru}

Z použití třetího simulátoru vyplývá další problém, který značně komplikuje práci s ním: V případě, že nenalezneme vhodné hodnoty pro počáteční stavový vektor částice v simulátoru $\boldsymbol{r}_0$, může se stát, že částice během průletu detektorem nabere sklon devadesáti stupňů k ose z, tedy ke směru paprsku. V takovém případě takto sestrojený simulátor selže, protože dále není schopen dráhu částice simulovat. Z experimentování s takovým simulátorem dále vyplývá, že nalézt takový stavový vektor (s vyjímkou triviálních případů, jako například částice s nulovým nábojem nebo s nulovým počátečním úhlem $t_x$) není jednoduché a pro naše potřeby generování velkého množství testovacích částic je nepraktický. Z tohoto důvodu bylo nakonec ze snahy o sjednocení modelů pro Kalmanův Filtr a generování částic upuštěno a zatímco Kalmanův Filtr používá model závislosti na souřadnici z (třetí model), trajektorie částic jsou generovány kombinovaným simulátorem s časovým modelem a úpravou Lorentzových rovnic pro relativistické rychlosti v homogenním elektromagnetickém poli.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Kapitola: Návrh a implementace rekonstrukčního algoritmu %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Návrh a implementace rekonstrukčního algoritmu}\label{reconstruction}
Tato kapitola se podrobně zabývá návrhem a implementací rekonstrukčního algoritmu pro sledování částic v experimentu CBM. Stručně jsou zde shrnuty použité nástroje a programovací jazyky, současné navržené řešení pro platformu \textit{Intel Xeon Phi} profesora Kisela a doktorky Akishiny, velký důraz je kladen na vlastní rekonstrukční algoritmus a jeho rozdíly oproti původnímu návrhu. Dále se v kapitole budeme zabývat možnostmi akcelerace tohoto algoritmu na hardwarových platformách, zejména pak paralelizací výpočtů pomocí protokolu MPI. V závěru kapitoly diskutujeme dosažené zrychlení v důsledku použité metody.

\section{Sekvenční verze rekonstrukčního algoritmu}\label{sequential_algorithm}
V této sekci se budeme zabývat plně sekvenční, nijak neparalelizovanou verzí rekonstrukčního algoritmu. V tomto algoritmu také identifikujeme výpočetně nejnáročnější součásti, na které se následně při paralelizaci v dalších sekcích zaměříme. 

Rekonstrukční algoritmus, jak již bylo zmíněno, sdružuje dohromady dvě první části rekonstrukce částic v detektorech: \textbf{hledání stop} (track finding) a \textbf{prokládání trajektorií křivkou} (track fitting). Tyto dvě části sdružujeme kvůli podstatě Kalmanova Filtru, který k hledání stop používáme: v rámci evoluce Kalmanova Filtru počítáme hodnotu $\chi^2$ (viz rovnice \ref{eq:kalman_chi_square}), která značí, jak moc odpovídá poslední naměřený stav systému poslednímu odhadu Kalmanova Filtru. Pomocí této hodnoty pak měříme, jak moc je pravděpodobné, že série vybraných zásahů patří ke stejné trajektorii. Zároveň však protože hlavní výstup Kalmanova filtru je co nejlepší odhad stavového vektoru, ve chvíli, kdy tímto způsobem nalezneme sadu zásahů, u kterých je nejpravděpodobnější, že byly vytvořeny stejnou trajektorií, jsme zároveň provedli druhou část rekonstrukce, tedy fitování trajektorie. Protože ale fyzikální model pohybu částice v detektoru není lineární, není možné použít klasický Kalmanův Filtr, který je optimální a stabilní, ale jeho rozšířenou podobu (EKF, viz \ref{kalman_filter}), která již optimální ani nutně stabilní není. Navíc, z důvodu šetření paměti a zvýšení rychlosti výpočtů, bylo rozhodnuto v původní implementaci algoritmu pro CBM používat při reprezentaci reálných čísel bitovou délku pouze 32 bitů, což dále potenciálně zvyšuje riziko divergence rozšířeného Kalmanova Filtru.

\subsection*{Inicializace algoritmu}
Vstupem rekonstrukčního algoritmu je neseřazené pole všech zásahů. Každý zásah je struktura dvou čísel s plovoucí desetinnou čárkou a jednoho celého čísla, které určuje hladinu detektoru:

\begin{lstlisting}
struct Hit
    x::Float32
    y::Float32
    z::UInt8
end
\end{lstlisting}

protože algoritmus používá pouze 8 detektorových hladin detektoru STS, které jsou rozestaveny v pravidelných intervalech, nabývá hodnota souřadnice $z$ ve strukturách zásahů hodnoty 3-10, díky čemuž nám stačí ji reprezentovat pouze jako celé číslo s nejmenší bitovou šířkou, jakou nám program dovolí použít, tedy 8 (jeden byte). V další části algoritmu budeme postupně všechny zásahy sdružovat do trojic a vyhodnocovat jejich $\chi^2$ hodnoty, podle kterých budeme dále vyhodnocovat nejpravděpodobnější vlastnosti jednotlivých trajektorií.

Nejdůležitější část této části algoritmu je určení první aproximace stavového vektoru $\boldsymbol{r}^+_0$ a kovarianční matice $C^+_0$. Určení kovarianční matice je triviální (viz rovnice \ref{eq:covariance_matrix_C_0}), vezmeme jednotkovou matici o velikosti 5 a vynásobíme arbitrálním velkým číslem $\mathbf{Inf}^2$. V případě $\boldsymbol{r}^+_0$ již ale postup triviální není. Můžeme například:

\begin{itemize}
    \item inicializovat první hodnotu $\boldsymbol{r}^+_0$ jako přímočarou trajektorii neutrálně nabité částice. V takovém případě inicializujeme $\boldsymbol{r}^+_0 = \{0,0,t_{x}, t_{y}, 0\}$ tak, že $t_{x} = x_0/z$ a $t_y = y_0/z$, kde $x_0, y_0$ jsou souřadnice $x$ a $y$ prvního měření, tedy prvního zásahu potenciální trojice. 

    \item inicializovat počáteční hodnoty $x$ a $y$  na 0 a $t_y$ inicializovat jako v prvním příkladu, $t_x$ pak získat metodou nejmenších čtverců tak, že pohyb trajektorie ve směru x proložíme parabolou. Hodnota $q/p$ je záporná v případě, že je první hodnota $x$ také záporná, jinak je kladná (neuvažujeme případ, kdy první hodnota $x$ je 0, protože středem detektoru prochází paprsek a detektorová hladina je v tomto intervalu prázdná).
\end{itemize}

V první fázi tedy algoritmus pro každý zásah od první po šestou hladinu detektoru (tedy pro hodnoty $z \in \langle3,8\rangle, z \in \mathbb{N}$, další hodnoty $z$ neuvažujeme, protože z nich nelze vytvořit trojice) generujeme první odhad stavového vektoru $\boldsymbol{r}^+_0$. V této práci byl nakonec použit první zmíněný postup, s jednou změnou: Experimentálně bylo zjištěno, že v případě předpokládaného nulového náboje se Kalmanův filtr brzy pokusí aproximovat pouze přímočarými trajektoriemi (tedy ustálit hodnotu $q/p$ jako 0). Tomuto fenoménu zabráníme tak, že určíme počáteční hodnotu $q/p$ jako $\pm 10^{-19}$ v závislosti na tom, jestli první naměřený bod má hodnotu $x$ kladnou, nebo zápodnou.

\subsection*{Hledání trojic zásahů}
Nejdůležitější (a časově nejnáročnější) součást algoritmu je hledání trojic zásahů, které nejpravděpodobněji patří ke stejné částici. Dr. Akishina ve své práci zmiňuje, že v závislosti na počtu částic, které v rámci daného snímku proletělí detektorem, zabírá fáze hledání trojic až $90.4\%$ celkového výpočetního času \cite{akishina4D}, a tedy vykazuje největší potenciál pro akceleraci.

Uveďme pro demonstraci této části algoritmu hypotetickou situaci, kdy detektorem prolétly tři částice, které zanechaly perfektní stopy (žádná částice neproletěla detektorem nezaznamenána) a nezanechaly šum. Pro každý bod, pro který hledáme nejvhodnější trojice, použijeme bod $[0,0]$ jako výchozí a počáteční odhad $\boldsymbol{r}^+_0$, který jsme získali v rámci předchozí fáze programu.

V tento okamžik provede algoritmus první Kalmanovskou Filtraci nad počátečním odhadem $\boldsymbol{r}^+_0$ a s prvním \uv{měřením} $\boldsymbol{m}_1$. Protože Filtr není lineární, je třeba nejprve získat správnou hodnotu propagační matice $F_k$. protože máme sestrojený simulační model, díky kterému dokážeme odsimulovat propagaci částice. Experimentálně bylo zjištěno, že většina částic, které detektorem proletí, je jen velmi málo ovlivněno magnetickým polem a jejich trajektorie je tak víceméně přímá. Z tohoto důvodu příští hodnotu propagační matice $F_k$ získáme tak, že za pomoci bodu $\boldsymbol{m}_1$ a odhadu $\boldsymbol{r}^+_0$ inicializujeme Runge-Kuttovu metodu 4. řádu s krokem 10cm, díky čemuž získáme stav příštího odhadu $\boldsymbol{r}^{lin}_1$ a dopočítáme podle nich $F_0$ podle následující rovnice:

\begin{equation}\label{eq:kalman_1_r_minus}
    \boldsymbol{r}^{lin}_1 = \boldsymbol{r}^+_0 F_0
\end{equation}

Zde je třeba poznamenat, že přestože k získání nového odhadu $\boldsymbol{r}_1^-$ nepotřebujeme linearizovanou propagační matici $F_1$, protože jsme ji získali pomocí Runge-Kuttovy metody zmíněné výše, propagační matici potřebujeme také proto, abychom mohli správně upravit kovarianční matici $C^-_1$ podle rovnice:

\begin{equation}\label{eq:kalman_1_C_minus}
    C^-_1 = F_0 C^+_0 F^T_0 + Q_0
\end{equation}

kde $F_0$ je lineární propagační matice, $C^+_0$ je původní kovarianční matice a $Q_0$ je matice procesního šumu, která je pro naše účely nulová, protože simulátor ani rekonstrukční algoritmus se šumem, vyplývajícím z mnohonásobného rozpadu při průchodu detektorovou hladinou nepočítají. V případě, že bychom chtěli pracovat s procesním šumem, který v praxi vzniká vlivem mnohonásobného rozpadu, ke kterému dochází při průchodu částice jednotlivými proužky detektoru, modelovali bychom ho následujícím způsobem: \cite{akishina4D}

\begin{equation}\label{eq:kalman_process_noise}
    Q_k = \sigma^2_k\left(\theta\right) \begin{bmatrix}
        0 & 0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 & 0 \\
        0 & 0 & (t_x^2 + 1)t_r^2 & t_xt_yt_r^2 & 0 \\
        0 & 0 & t_xt_yt_r^2 & (t_y^2 + 1)t_r^2 & 0 \\
        0 & 0 & 0 & 0 & 0
    \end{bmatrix}
\end{equation}
kde $\theta$ je úhel dopadu částice na hladinu detektoru. $\sigma^2_k$ značí směrodatnou odchylku normálního rozložení, určujícího míru změny směru trajektorie v důsledku mnohonásobného rozpadu. Takovou směrodatnou odchylku určíme následujícím způsobem: \cite{akishina4D}

\begin{equation}\label{eq:kalman_multiple_scattering}
    \sigma\left(\theta\right) = \frac{13.6 \text{MeV}}{\beta cp}q\sqrt{L/X_0}\left[1 + 0.038 \ln{\left(L/X_0\right)}\right]
\end{equation}

Přestože výpočet směrodatné odchylky pro procesní šum není nemožný, narazíme na hned několik problémů, jejichž řešení na první pohled není zřejmé, a pro naše účely zbytečně komplikované:

\begin{itemize}
    \item neznáme hodnotu rychlosti částice $\beta$, protože námi používaný fyzikální model s rychlostními vektory schválně nepočítá, abychom se vyhnuli problémům s časovou proměnnou v Kalmanově Filtraci.
    \item nemáme k dispozici materiálové vlastnosti jednotlivých proužků detektoru, což výpočet mnohonásobného rozpadu dále komplikuje.
\end{itemize}

Z tohoto důvodu bylo od zahrnutí mnohonásobného rozpadu v detektoru jak ze simulačního, tak rekonstrukčního modelu vypuštěno.

Zbytek Kalmanovy filtrace je totožný s postupem zmíněným výše v práci (viz \ref{kalman_filter}). Na konci této filtrace má algoritmus k dispozici \uv{dvojice} bodů, kde první bod je střed terče a druhý bod je výchozí bod pro stavbu trojice zásahů. Tento postup opakujeme tak, že pro každý výchozí bod vyzkoušíme dále Kalmanovu filtraci s každým bodem z následující detektorové hladiny jako další měření. Celý tento proces pak zopakujeme ještě jednou tak, že na konci této fáze programu vznikne $a*b*c$ možných trojic zásahů, kde $a$ je počet zásahů na detektorové hladině $i$, $b$ je celkový počet zásahů na detektorové hladině $i+1$ a $c$ je celkový počet zásahů na detektorové hladině $i+2$. Každá z těchto trojic má v jazyce Julia následující strukturu:

\begin{lstlisting}
    struct Hit
        x::Float32
        y::Float32
    end
    
    struct Triplet
        point1::Hit
        point2::Hit
        point3::Hit
        chi_sq::Float32
    end
\end{lstlisting}

Zde je dobré poznamenat, že ne všechny kombinace zásahů jsou dále používány. Díky hodnotě $\chi^2$ je na první pohled zřejmé, že některé trojice (i dvojice) zásahů není možné, aby tvořily dráhu jedné trajektorie. Z tohoto důvodu lze experimentálně nastavit práh, který nesmí hodnota $\chi^2$ překročit, nebo bude daná dvojice či trojice zahozena a dále s ní není pracováno. Tímto způsobem se vyhneme nutnosti zpracovávat již tak velké množství dat vyloženě kombinatoricky a přispějeme tak k výraznému zrychlení.

\subsection*{Vytvoření a evoluce celulárního automatu}
V okamžiku, kdy program má k dispozici všechna data o možných trojicích zásahů, je možné zkonstruovat celulární automat, díky kterému budeme moci spojit jednotlivé trojice zásahů do kompletních kandidátů na trajektorie.

Narozdíl ale od intuitivního pojetí celulárního automatu, kdy každý pixel reprezentuje jedna buňka automatu a okolí buňky je obdobně složeno z okolních pixelů detektoru, modelujeme celulární automat následujícím způsobem:

\begin{itemize}
	\item každá buňka automatu reprezentuje možnou trojici zásahů, vygenerovanou předchozím krokem programu. Jakmile se program dostane do této fáze, již jednotlivé zásahy nezpracovává, a pouze skládá nejnadějnější trojice zásahů.

	\item Okolí jednotlivých buňek určujeme následovně: pokud dvě trojice mají společné dva zásahy (a neliší se prostředním z nich, potom jsou trojice konkurenti) a podobné hodnoty $q/p$, jsou označeny jako sousední. 
\end{itemize}

Jakmile vytvoříme síť sousedů, vzniká nám jakási stromová struktura, se kterou pracujeme dále: každé buňce (tj. trojici) přiřadíme \uv{počítadlo}, které počítá nejdelší řetězec sousedů směrem k terči (tj. určujeme nejdelší kandidáty trajektorií). Ve chvíli, kdy máme inicializované \uv{buňky} se správně přiřazenými počítadly, můžeme začít vybírat nejlepší kandidáty na částice. Výběr je nyní již jednoduchý: vzhledem k tomu, že jako sousedy trojice uvažujeme pouze tehdy, kdy jsou jejich stavové vektory a hodnoty $\chi^2$ dostatečně podobné, díky čemuž dokážeme vyloučit většinu náhodných kombinací trojic, v tento okamžik vybíráme nejlepší kandidáty podle dvou kritérií:

\begin{itemize}
	\item od nejdelší trajektorie po nejkratší: Pravděpodobnost, že náhodná množina zásahů \textbf{ne}reprezentuje reálný průlet částice exponenciálně klesá s každým dalším zásahem, který by mohl trajektorii odpovídat. Z tohoto důvodu algoritmus upřednostňuje ďelší čtrajektorie před kratšími. 

	\item od nejnižší hodnoty $\chi^2$ po nejvyšší: Čím vyšší hodnota $\chi^2$, tím méně sada vybraných zásahů odpovídá chování reálné potenciální částice. Z tohoto důvodu upřednostňuje program kandidáty s nižšími hodnotami $\chi^2$ před vyššími.

	\item pokud trojice zásahů již byla v trajektorii dříve použita, podruhé již použita být nesmí. Z tohoto důvodu, než program definitivně reportuje nalezenou trajektorii, zkontroluje, jestli některá z použitých trojic nebyla již použita. Pokud ano, tento kandidát není použit a je zahozen.
\end{itemize}

\subsection*{Výběr nejlepších kandidátů a příští evoluce celulárního automatu}
Výstupem předchozí fáze programu je seznam kandidátů na \uv{nejnadějnější} trajektorie, seřazený sestupně podle délky trajektorie a hodnoty $\chi^2$, sestupně v tomto pořadí. V této fázi program projde sekvenčně všechny kandidáty a pokud jsou všechny body, které kandidát obsahuje, dosud nepoužity, zapíše je na výstup, v opačném případě je zahodí. Jakmile program takto projde celým seznamem, tato fáze programu je ukončena, použité zásahy jsou ze seznamu odstraněny a celulární automat je inicializován znovu pro další evoluční fázi, tentokrát s jinými parametry. Dále algoritmus probíhá stejným způsobem, tentokrát pro trajektorie s nízkou energií, poté pro sekundární trajektorie a nakonec pro trajektorie, u kterých kvůli nedokonalosti detektoru STS došlo k \uv{vypadnutí} jednoho zásahu (pravděpodobnost, že v jedné trajektorii vypadne více než jeden zásah je menší než $0.05\%$, což je vzhledem k vlastnostem detektoru zanedbatelné). Vzhledem ale k tomu, že se v rámci této práce zabýváme pouze primárními trajektoriemi s vysokou energií (a pro tyto účely vytvořený simulátor vytváří pouze takové trajektorie) a další iterace celulárního automatu jsou k té první analogické, dále se jimi nebudeme zabývat. 

\section{Rozdíly v implementaci oproti původnímu algoritmu}\label{akishina_vs_roth}
Jak již bylo zmíněno výše, pro tuto práci vyvinutý algoritmus je oproti původnímu algoritmu z práce Dr. Akishiny \cite{akishina4D} zjednodušený tak, aby zachoval podstatu programu spolu s jeho nejnáročnějšími součástmi, byl ale zjednodušen tak, aby se nemusel zabývat pokročilými fyzikálními fenomény, které se jsou z programátorského hlediska sice de facto analogické, ale z fyzikálního hlediska výrazně nad rámec této práce. Jak jsme již zmínili výše, simulátor i rekonstrukční algoritmus vyvinutý v rámci této práce se dopouští zjednodušení zejména v oblasti procesního šumu, který ani simulátor, ani rekonstrukční algoritmus neuvažují, tento fakt ale z důvodu, že je zanedbán již ve fázi generování trajektorií nemá na funkci algoritmu vliv (můžeme říci, že model použitý pro rekonstrukční algoritmus je pouze přesnější). 

Další rozdíl je oproti původní implementaci je v omezení implementace pouze na primární trajektorie s vysokou energií, zčásti proto, že jejich rekonstrukce je de facto analogická, ale také proto, že zejména v případě sekundárních trajektorií nastává velký problém při simulaci takových rozpadů, protože takový simulační model by obnášel výrazně složitější model, který vyžaduje hluboké znalosti z oblasti částicové fyziky, což opět přesahuje rámec této práce.

Jako poslední významný rozdíl oproti původní implementaci zmiňme nástroje a technologie při implementaci použité, a technologický kontext práce. V rámci této práce vznikl rekonstrukční program, který jako hlavní způsob akcelerace využívá zejména distribuované počítání, zejména pomocí protokolu MPI\footnote{MPI -- Message Passing Interface (rozhraní pro předávání zpráv)} a modulu Distributed.jl jazyka Julia, který oproti MPI\footnote{Zde rozumíme standard MPI-1, který neposkytuje mechanismy pro jednostrannou komunikaci RMA} programátorsky zjednodušuje \uv{jednostranné} posílání zpráv pro facilitaci distribuovaných výpočtů. Oproti tomu původní algoritmus pro rekonstrukci částic v experimentu CBM byl vytvořen tak, aby co nejlépe využíval platformy \textit{Intel Xeon Phi}. Silnou stránkou této platformy bylo velké množství CPU jader s instrukční sadou rodiny \textbf{x86} (Intel Xeon Phi 7295, původně vydaný v roce 2017, disponoval 72 jádry o čtyřech vláknech na jádro, tedy 288 vláken \cite{cutress_2017}) a podpora 512-bitových SIMD\footnote{SIMD -- Single Instruction, Multiple Data (jedna instrukce, více dat)} instrukcí AVX-512. Z tohoto důvodu se původní algoritmus zabýval paralelizací na úrovni vláken a vektorizací výpočtů pomocí SIMD instrukcí, které byly silnou stránkou této platformy. Vývoj této platformy byl ale v roce 2019 (částečně kvůli problematické litografické 10nm technologii společnosti Intel a velkou konkurencí ze strany výrobců GPGPU) ukončen a platforma označena jako \textit{End of Life}. \cite{cutress_shilov_2019}

\section{Paralelní rekonstrukční algoritmus}\label{parallel_algorithm}
V této podkapitole se budeme podrobně zabývat paralelizací algoritmu, jehož části byly v této práci již podrobně rozebrány. V této práci se zaměřujeme primárně na paralelizaci náročných částí programu pomocí distribuovaných výpočtů z toho důvodu, že distribuovaný přístup umožňuje nejen paralelizovat výpočty na jednotlivých jádrech procesoru, ale zejména \textit{distribuovat} výpočty mezi různými počítači na síti. Tento přístup je dále podpořen faktem, že datová centra mají zpravidla velmi dobře optimalizovanou infrastrukturu pro distribuované výpočty (dnešní počítače lze například pomocí technologie InfiniBand propojit za přenosových rychlostí až 1200Gb/s s latencí nižší než 600 nanosekund \cite{mellanox}), v důsledku čehož je možné využít velkého množství výpočetního výkonu pro danou aplikaci, což by v extrémních případech, jako je on-line rekonstrukce trajektorií částic v HEP experimentech, často nepraktické kvůli velké režii odesílání dat po síti pomocí protokolů jako MPI.

\subsection*{Analýza paralelizovatelných součástí programu}\label{speedup_potential}
Abychom byli schopni vysoce náročný program efektivně akcelerovat, je nejprve potřeba provést rozbor jeho nejnáročnějších částí. Programovací jazyk Julia pro usnadnění profilování zdrojového kódu interaktivní \textit{plamenové grafy}\footnote{z angl. flame graphs}, graf vykreslující stav zásobníku volání v závislosti na čase. Modul \texttt{FlameGraphs.jl} navíc pomáhá nalézt kritické části programu, které mohou způsobovat problémy například z důvodu neudání explicitních typů proměnných, což může do značné míry omezit schopnost kompilátoru optimalizovat kód před spuštěním.

\begin{figure}
    \centering
    \includegraphics[width=0.8\textwidth]{obrazky-figures/flame_graph.png}
    \caption{Plamenový graf funkce \texttt{run\_reconstruction()} nad datasetem \texttt{cbm\_10\_particles}, vygenerovaným simulátorem s relativistickým časovým modelem. Za povšimnutí stojí červeňe vybarvená oblast z volání funkce nad datovými rámci, kde jsme se nevyhli problémům s dynamickým typováním, časová ztráta způsobená touto chybou je ale vedle ostatních součástí programu zanedbatelná.} 
    \label{fig:flame_graph}
\end{figure}

Z grafu \ref{fig:flame_graph} můžeme vyčíst, že zdaleka nejnáročnější součást programu je vlastní smyčka Kalmanova Filtru. Výrazně tak převáží i časovou náročnost evoluce celulárního automatu i ostatní části programu, jako například operace úprav datasetu pro potřeby programu. Zároveň stojí za povšimnutí, že nejnáročnější část funkce není ani samotné maticové násobení, protože (podobně jako v případě vektorizovaných operací nad datasety) samotné maticové výpočty Julia při dodržení základních pravidel správného návrhu velmi dobře optimalizuje a výrazně tak snižuje atraktivitu manuální vektorizace kódu na tomto místě. Protože je ale funkce navržená tak, aby mohla být spouštěna mnohokrát s na sobě vzájemně nezávislými argumenty, je zřejmé, že paralelizovat bude výhodnější na úrovni o jedna vyšší, tedy ve funkci \texttt{combination\_search()}.

Zároveň je třeba poznamenat, že časová náročnost s počtem zásahů v detektoru stoupá u části vypočítávající odhady Kalmanova Filtru výrazně rychleji, než u ostatních částí programu, jak vyplývá z kombinatorické podstaty algoritmu konstrukce trojic. Dr. Akishina ve své práci zmiňuje, že konstrukce trojic zásahů zabírá přes 90\% výpočetního času programu \cite{akishina4D}, a i když v našem případě tento údaj není tak markantní z důvodů omezení, které blíže popisujeme v podkapitole \ref{simulator} (zejména pak nehomogenní magnetické pole, se kterým současné řešení počítá, výrazně zvyšuje dobu, kterou stráví ve smyčce Kalmanova Filtru, konkrétně pak u propagace $\boldsymbol{r}^t$ Runge-Kuttovou metodou), stále se jedná o nejnáročnější součást výpočtu a tedy i o součást, které věnujeme při paralelizaci největší pozornost.

\begin{figure}
    \centering
    \begin{verbatim}
        BenchmarkTools.Trial: 1880 samples with 10 evaluations.
        Range (min … max):  1.816 ms … 16.123 ms
        Time  (median):     2.085 ms
        Time  (mean ± sig):   2.655 ms ±  1.880 ms
        GC (min … max): 0.00% … 79.34%
        GC (median):    0.00%
        GC (mean ± sig):  9.32% ± 11.36%
        Memory estimate: 2.64 MiB, allocs estimate: 29229.
    \end{verbatim}
    \caption{Výstup benchmarku nástroje \texttt{BenchmarkTools} jazyka Julia pro sekvenční implementaci rekonstrukčního algoritmu. Jako testovací data byl použit dataset \texttt{cbm\_3\_particles}, vygenerovaný simulátorem s relativistickým časovým modelem. Modul BenchmarkTools v rámci makra \texttt{@benchmark} vyhodnotí rychlost provádění měřené funkce a vygeneruje statistiky a histogram jednotlivých iterací.}
    \label{fig:seq_3_particles}
\end{figure}

\subsection*{Paralelizace rekonstrukčního programu}

V předchozí kapitole jsme zjistili, že největší potenciál pro zrychlení má u rekonstrukčního programu smyčka Kalmanova Filtru. Tento fakt nám radí, že pro vhodnou paralelizaci potřebujeme najít způsob, jakým rozdělit \uv{stavový prostor} problému tak, abychom minimalizovali množství dat, se kterými jednotlivé uzly bezpodmínečně potřebují pracovat. V případě paralelizace na úrovní vláken nebo logických procesorů by nám pomohla možnost pracovat se sdílenou pamětí, ať už v podobě vyrovnávací paměti \textit{cache} nebo RAM. V případě distribuovaného počítání je ale nutné se takovému paradigmatu pokud možno vyhnout: jednotlivé uzly (většinou, existují výjimky, viz \cite{mellanox}) nesdílí paměť mezi sebou, datová komunikace probíhá pomocí zpráv, které jsou oproti komunikaci mezi např. jednotlivými jádry procesoru velmi pomalé, a tedy drahé. Z tohoto důvodu je žádoucí, aby každý uzel v rámci výpočetního \textit{clusteru} fungoval pokud možno co nejvíce samostatně a teprve na konci své práce odevzdal vedoucímu uzlu výsledek.

V rámci bližšího zkoumání implementovaného sekvenčního algoritmu byly identifikovány čtyři různé body, potenciálně vhodné k paralelizaci: 

\begin{itemize}
    \item Dva v rámci funkce \texttt{combination\_search()} (dvě náročné, vnořené \texttt{for} smyčky) a
    \item Dvě o stupeň v zásobníku volání níž ve funkci \texttt{reconstruct\_tracks}vnořené \texttt{for} smyčky.
\end{itemize} 

V rámci hledání optimálního řešení byly tedy implementovány všechny čtyři potenciální paralelizace. Pro paralelizaci na bázi distribuovaného počítání bylo třeba provést následující kroky:

\begin{enumerate}
    \item \textbf{Načtení modulu \texttt{Distributed}:} Než nám bude umožněno jakkoli manipulovat s dalšími procesy, musíme nejdřív v mateřském procesu načíst tento modul. Tento modul nám poskytne dvě extrémně pro nás důlkežitá makra: \texttt{@everywhere}, které akci, kterou chceme provést, vynutí na všech uzlech clusteru, a \texttt{@distributed}, makro, které nám umožní bez velkých zásahů do kódu náročné výpočty distribuovat.

    \item \textbf{Vytvoření nových procesů:} Pokud jsme nespustili prostředí Julia pomocí přepínače \texttt{julia -p N}, kde \texttt{N} značí počet procesů, které se k tomu hlavnímu vytvoří \textbf{navíc}, je třeba tak učinit nyní. Učiníme tak prostě pomocí funkce \texttt{addprocs(N)}, kde \texttt{N} opět značí počet dalších procesů, které je třeba vytvořit (zde sice hovoříme o procesech, protože je instanciujeme na stejném počítači, implementace \texttt{Distributed} ale počítá i s komunikací jednotlivých uzlů po síti.

    \item \textbf{Nastavení prostředí nových procesů}: Julia při vytváření nového procesu nekopíruje proměnné prostředí, proces je tak inicializován \textit{de facto} prázdný. Proto je potřeba dceřinným procesům předat všechny důležité informace předtím, než program spustíme a začneme měřit jeho výkon. Zde nám pomůže makro \texttt{@everywhere}, díky kterému standardní operace pro načtení knihovny \texttt{inlude("CBMReconstructor.jl")} nemusí být nijak speciálně upravena ani nemusíme nijak měnit strukturu původního programu, makro problém nových uzlů vyřeší za nás. Kromě knihovny \texttt{CBMReconstructor} je daké potřeba informovat uzly, aby načetly modul \texttt{DistributedArrays}. Tento modul řeší problém, kdy jednotlivé procesy potřebují pracovat nad stejnou proměnnou (například ukládat struktury do stejného pole).

    \item \textbf{Paralelizace náročné smyčky pomocí \texttt{@distributed}}: Po tom, co jsme úspěšně inicializovali cluster a identifikovali nejperspektivnější části programu pro paralelizaci, vlastní rozdělení práce mezi uzly je opět triviální. Pomocí makra \texttt{@distributed} program rozdělí všem uzlům práci a implicitně pokračuje ve vykonávání programu (pokud nepoužijeme explicitně makro vynucující synchronizaci, což pro naše potřeby není nutné), dokud nepotřebuje pracovat s daty, o kterých ví, že zpracovává jiný proces a dosud ji nevrátil (v takovém případě uzel čeká).

    \item \textbf{Vrácení výsledků paralelních bloků a změření rychlosti programu}: V okamžiku, kdy ukončíme provádění paralelní části programu, provedeme zbytek programu stejně, jako v případě sekvenční verze programu. V zájmu dašího zrychlení projektu byl navržen algoritmus pro paralelní/distribuované řazení vybenerovaných kandidátních trajektorií. Navržena byla implementace distribuovaného \textit{merge sortu}, po profilování programu a vygenerování plamenového grafu bylo ale od pokusu implementovat distribuovaný merge sort upuštěno ze dvou důvodů: Zaprvé Julia implementuje velmi dobře optimalizovanou \textit{in-place} funkci, zadruhé u funkce, která není prováděna dostatečně dlouho na to, aby se paralelizace vyplatila, nevhodné vynucení distribuovaného počítání mohou rychlost zpracování naopak dramaticky zpomalit.
\end{enumerate}

Po tom, co jsme úspěšně vytvořili snadno distribuovatelný program, zbývá jen změřit a vyhodnotit výsledky a zjistit tak, která z možností paralelizace je nejvýhodnější.

\subsection*{Porovnání výkonu paralelních algoritmů}
\begin{figure}
    \centering
    \includegraphics[width=0.8\textwidth]{obrazky-figures/speedup table.png}
    \caption{Tabulka měření jednotlivých verzí distribuovaného rekonstrukčního programu. Z dat můžeme vyvést závěr, že distribuce vnitřní smyčky funkce pro rekonstrukci trajektorie je pro zrychlení (a šetření systémovou pamětí). Zatímco některé z implementačních verzí se dle očekávání zrychlovaly, jiné naopak výrazně na výkonu ztrácely i vůči naivní, sekvenční verzi programu.}
    \label{fig:parallel_table}
\end{figure}

Z tabulky \ref{fig:parallel_table} můžeme vyčíst několik zajímavých faktů. Zaprvé, jak vyplývá z plamenového grafu, jednotlivé cykly, které jsme se snažili paralelizovat, jsou vnořené jedna do druhé (i když je tento fakt zakryt faktem, že se objevují v různých funkcích. Zatímco sekvenční algoritmus je neefektivní a pro praktické využití nepoužitelný kvůli své vysoké náročnosti na paměť a malé rychlosti, pokus paralelizovat \uv{nejvnitřnější} smyčku z těchto tří způsobí, že nedostatek výpočetně náročných instrukcí, které by použití distribuovaného počítání ospravedlnilo a relativně vysoká režie posílání zpráv jiným uzlům způsobil, že program se s vyšším počtem jader stává výrazně pomalejším a \uv{hladovějším} po paměti. Není třeba zdůrazňovat, že tento postup není vhodný.

Vnější cyklus v téže funkci naproti tomu už se chová daleko víc, jak bychom si představovali, s výrazně nižšími nároky na paměť a až čtyřnásobné zrychlení oproti sekvenčnímu algoritmu. Nicméně stojí za zmínku, že ve chvíli, kdy uzel přestane být plně utilizovaný, rychle přestává být užitečný. V případě vnější smyčky funkce \texttt{combination\_search()} takový jev sledujeme u měření s dvaceti čtyřmi logickými procesory: řádově vyšší spotřeba paměti a nižší rychlost, než kdybychom použili pouze šest jader. Znovu se tedy prokazuje, že distribuované počítání je výhodné pouze tehdy, máme-li dobrý důvod ho použít.

Vnější cyklus vnější funkce (tedy \uv{nejvnějšejší} cyklus) je další dobrý příklad této teorie: Protože je v této smyčce přesně dané, že bude iterovat pouze v intervalu $[3:8]$. Z toho vyplývá, proč se zvyšujícím se množstvím dostupných jader se výkon algoritmu zhoršuje. Za povšimnutí každopádně stojí velmi nízký počet alokací v programu, což je šetrné k paměti.

Poslední (a nejlepší) varianta řešení distribuce výpočtu rekonstrukčního algoritmu je vnitřní smyčka vnější funkce. Z důvodů, které jsme zmínili výše, byla tato varianta jediná, která při spuštění na plných 24 logických procesorech byla schopná pojmout plný výkon procesoru, díky čemuž poskytuje nejvyšší zrychlení téměř 800\%. Opět ale stojí za zmínku, že přestože má s vyšším počtem jader stále vyšší výkon, neznamená to, že můžeme paralelizovat \uv{donekonečna}: již na rozdílu mezi časem naměřeným pro 16 a 24 jader je vidět, že s vyšším počtem jader je každé další jádro pro výpočet méně přínosné.

\section{Alternativy k použití distribuovaného počítání}\label{alternative_accelerators}

V raných fázích této fáze bylo uvažováno nad akcelerací zmiňovaného algoritmu na platfomě FPGA. Tento přístup se ale ukázal být nepraktický z několika důvodů, ze kterých zmiňme dva nejmarkantnější:

\begin{itemize}
    \item \textbf{Paměťová náročnost.} Jeden z hlavních problémů rekonstručního programu je jeho paměťová náročnost. Již při rekonstrukci 100 částic pracuje program s vysokými desítkami gigabajtů dat, což způsobuje problém, kdy provoz potenciálního FPGA akcelerátoru vyžaduje časté I/O (jakýsi \uv{page miss}) operace, které dále zpomalují výpočet. Tento problém platformy FPGA je těžké odstranit a uvážíme-li, že reálný program je kvůli přítomnosti nehomogenního magnetického pole ještě výrazně paměťově náročnější, alternativa FPGA přestává být atraktivní.
    
    \item \textbf{Prostorová náročnost.} Z podstaty maticového násobení, které je fundamentální pro rekonstrukční algoritmy založené na Kalmanově Filtru, vyplývá, že jeho implementace na hradlových políc by byla velmi náročná. Při násobení matic bychom pro každý výpočet matic $5\times5$ potřebovali 125 násobiček a neexistuje způsob, jak takové množství zmenšit. Když navíc uvážíme, že takových násobení je \textbf{v každé iteraci} Kalmanova Filtru dvanáct, začne být zřejmé, že jakýkoli vysoce specializovaný obvod bude zabírat tolik místa na hradlovém poli, že je jeho implentace nepraktická. Navíc, specializované akcelerátory jako \textit{GPGPU}\footnote{GPGPU -- General Purpose Graphics Processing Unit (obecný grafický akcelerátor)}, poskytují vysoký výkon pro maticové operace (zejména pak násobení a sčítání) a relativně mírnou učící křivku pro jejich programování.
\end{itemize}

Stejně jako algoritmus zvolený pro tuto práci je současné řešení rekonstrukce trajektorií postavené na konceptech celulárního automatu a rozšířeného Kalmanova filtru. Konkrétně řešení navržené pro experiment CBM bylo navrženo pro platformu \textit{Intel Xeon Phi} a využívá paralelního zpracování pomocí vláknové paralelizace a vektorizace SIMD instrukcemi. Tento postup má oproti akceleraci na FPGA velkou výhodu, a tou je jeho \textbf{operační paměť}: Největší zřejmá slabina hradlových polí v této situaci je velmi malé možnosti operační paměti. 

Naopak velmi perspektivní je v tomto odvětví použití grafických procesorů. Z podstaty Kalmanova Filtru vyplývá, že využívá zejména maticového sčítání a násobení, dvě operace, na které jsou moderní grafické procesory velmi dobře vybaveny. V současné době v Evropském Centru pro Jaderný Výzkum (CERN) pracuje experiment ALICE\footnote{ALICE -- A Large Heavy Ion Experiment (velký experiment s těžkými ionty)}, jeden ze čtyř velkých experimentů na Velkém Hadronovém Urychlovači (LHC) a experiment, jehož výpočetní přístup je tomu v CBM velmi blízko (implementace na bázi celulárního automatu a Kalmanova Filtru). Rozdíl oproti CBM je v tom, že výpočty při ALICE jsou akcelerovány pomocí GPGPU \cite{ALICE} a vzhledem k tomu, že v rámci této práce nebylo optimalizačně nijak zasahováno do algoritmu Kalmanova Filtru nabízí otázku, zda by nebylo možné akcelerovat tento výpočet dále pomocí kombinace distrbuovaných výpočtů \textbf{a} urychlení \texttt{KF\_loop} na GPU. 

\chapter{Závěr}

Cílem této práce bylo seznámit se s moderními experimenty na půdě částicové a kvantové fyziky, podrobně prostudovat použité technologie a algoritmy využívané pro rekonstrukci trajektorií v detektorech částic, navrhnout a implementovat akcelerovatelný rekonstrukční algoritmus a tuto akceleraci následně také implementovat. Vzhledem ke rozsáhlosti tématu a matematicko-fyzikálnímu aparátu mimo rámec informačních technologií si tato práce dává také za cíl problematiku rekonstrukce částic v moderních fyzikálních experimentech co nejlépe osvětlit..

V rámci práce jsem se podrobně seznámil s jednotlivými fázemi zpracování dat z detektorů částic, různými technologiemi pro získávání dat z nich a algoritmy pro následnou rekonstrukci trajektorií částic, zejména pak pro detektor CBM při výzkumném středisku FAIR v Darmstadtu. Pro podrobnější studium jsem zvolil dvě výpočetně nejnáročnější (a tedy nejvhodnější pro studium možností jejich akcelerace) části zpracování fyzikálních dat z detektorů: rekonstrukce trajektorií částic a odhad (\textit{fitování}) vlastností takto rekonstruovaných trajektorií. Pro akceleraci těchto dvou fází jsem vybral rekonstrukční algoritmus na bázi celulárních automatů a Kalmanových filtrů a po prostudování stávajících řešení pro detektor CBM a detektory mu konceptuálně podobné jsem se rozhodl podrobně zabývat možnostmi akcelerace algoritmu pomocí distribuovaného počítání. 

Vzhledem k neexistenci vhodných testovacích dat (ať už reálných, či simulovaných) bylo třeba vhodná testovací data vytvořit. Z tohoto důvodu byly v průběhu práce pro různé účely postupně navrženy a implementovány čtyři různé verze simulačního modelu chování částic uvnitř detektorů se silným magnetickým polem. Více modelů bylo implementováno z důvodu, že samotný sledovací algoritmus potřebuje používat simulační model pro výpočet predikčního kroku rozšířeného Kalmanova filtru, a bylo tedy potřeba simulační model navrhnout tak, aby potřebám algoritmu co nejlépe odpovídal. Z tohoto důvodu jsem v různých fázích studia rekonstrukčního algoritmu průběžně podobu modelu upravoval tak, aby nejenže vhodně simuloval dráhy částic v detektoru, ale aby umožňoval i \uv{částečnou} simulaci (tj. aby bylo možné simulovat pouze zadaný úsek dráhy částice). Navzdory snaze sjednotit fyzikální modely tak, aby simulace trajektorií částic i propagační model Kalmanova Filtru fungovaly na stejném matematickém modelu, bylo od této snahy zdůvodu nestability q/p modelu při simulaci ostrých úhlů upuštěno a pro potřeby práce jsou používány dva.

V rámci studia rekonstrukčních algoritmů byl algoritmus implementován v jazyce Julia, nejdříve jako čistě sekvenční program a následně ve čtyřech verzích paralelizován pomocí distribuovaného počítání. Pro dataset nedokonalých částic s chybovým detektorem program oproti sekvenčnímu vykonává výpočet na 24 logických jádrech s téměř osmisetprocentním zrychlením. Jako další pokračování práce se nabízí implementace matice procesního šumu pro výpočet mnohonásobného rozpadu, prozkoumání možností dalšího zrychlení akcelerací distribuovaného rekonstrukčního algoritmu a Kalmanova Filtru za pomocí GPGPU, srovnání výkonu rekonstrukčního kódu v Jazyce Julia vůči například C/C++, popřípadě další zdokonalení vlastního algoritmu, ať už zahrnutím nehomogenního magnetického pole, dělení detektorových stanic podle počtu zásahů nebo implementace rekonstrukce nově vzniknuvších nebo rozbitých částic (s chybějícím zásahem na jedné detektorové stanici).