Let us define the real functions $f_a(t)=\begin{cases} 1 & \text{ if } t \in \left[-a;a\right] \\ 0 & \text{ else} \end{cases}$ and $g_{a;b}(t)=\begin{cases} \cos\left(2\pi bt\right) & \text{ if } t \in \left[-a;a\right] \\ 0 & \text{ else} \end{cases}$.
Since $g_{a;b}\in L^1(\mathbb{R})$, we can compute the Fourier Transform of $g_{a;b}$, and we obtain $\mathcal{F}(g)=2a.\operatorname{sinc}(2\pi a(x-b) ) + 2a.\operatorname{sinc}(2\pi a(x+b))$.
We can notice that this transformation is also the convolution product between a distribution and a function in $\mathcal{C}^\infty(\mathbb{R})$: $2a.\operatorname{sinc}(2\pi a(x-b) ) + 2a.\operatorname{sinc}(2\pi a(x+b))= <\delta_b+\delta_{-b};y\mapsto2a.\operatorname{sinc}(2\pi a(x-y) )>_y= (\delta_b+\delta_{-b}) \ast 2a.\operatorname{sinc}(2\pi a\ \cdot) = (\delta_b+\delta_{-b}) \ast \mathcal(F)(f_{a;b}) $
We also, know that the Dirac distribution is linked to the regular distibution associated to the cosine function by $\mathcal{F}T_{\cos(2\pi b\ \cdot)}=\delta_b+\delta_{-b}$.
Therefore, I wander if there is a proper way to obtain the Fourier Transform of $g_{a;b}$ by using the distribution theory?
Maybe there is a theorem for a the product of a function in $L^1(\mathbb{R})$ (I mean $f_{a;b}$) and a distribution. I know one between a Schwartz's function and a tempered distribution, but $f_{a;b}\notin\mathcal{S}(\mathbb{R})$.
One cannot have a result for the product of a distribution with an $L^1$ function. For example if $T$ is the distribution induced by the $L^1$ function $\mathbf 1_A/\sqrt{|x|}$, where $\mathbf 1_A$ is an indicator of a bounded open set containing 0, then note that what is essentially $T^2$, $(\mathbf 1_A/\sqrt{|x|}) T$ is not a distribution (if one tries to follow the usual definition of a product for distributions).
In your particular case, writing $g:=g_{a,b}$ for legibility reasons, we have the expression for $g$, $$g=g_1g_2, \quad g_1(t) = \cos(2π bt) , \quad g_2(t) = \mathbf1_{t\in[-a,a]}, $$which is a product of $g_2\in L^2\subset \mathcal S'$ with $g_1\in C^\infty$ and hence as distributions, $$ \hat g = \hat {g_1}*\hat{g_2} = (\delta_b+\delta_{-b}) * \operatorname{sinc}(2π ax) $$ which gives your result. With a little bit more care, the above line is because: \begin{align} \langle \hat g , \phi\rangle &= \langle g_1g_2 , \hat {\phi}\rangle \quad\text{definition of F.T.}\\ &= \langle g_2 , g_1\hat \phi\rangle \quad\text{definition of product by $C^\infty$}\\ &= \langle g_2 , (\overset{\vee}{g_1} * \phi)^\wedge \rangle \quad\text{property of F.T.}\\ &= \langle \hat {g_2} , \overset{\vee}{g_1} * \phi \rangle \quad\text{definition of F.T.} \\ &= \langle \hat{g_2} * \hat{g_1},\phi\rangle \quad\text{definition of a convolution of a distribution} \end{align} For the last line, see for example this MSE link.