  We present stationary solutions of geometrically thick discs (or tori)
endowed with a self-consistent toroidal magnetic field distribution surrounding
a non-rotating black hole in an analytical, static, spherically-symmetric
$f(R)$-gravity background. These $f(R)$-gravity models introduce a Yukawa-like
modification to the Newtonian potential, encoded in a single parameter $\delta$
which controls the strength of the modified potential and whose specific values
affect the disc configurations when compared to the general relativistic case.
Our models span different magnetic field strengths, from purely hydrodynamical
discs to highly magnetized tori. The characteristics of the solutions are
identified by analyzing the central density, mass, geometrical size, angular
size, and the black hole metric deviations from the Schwarzschild space-time.
In the general relativistic limit ($\delta=0$) our models reproduce previous
results for a Schwarzschild black hole. For small values of the $\delta$
parameter, corresponding to $\sim 10\%$ deviations from general relativity, we
find variations of $\sim 2 \%$ in the event horizon size, a $\sim 5\%$ shift in
the location of the inner edge and center of the disc, while the outer edge
increases by $\sim 10\%$. Our analysis for $|\delta|>0.1$, however, reveals
notable changes in the black hole space-time solution which have a major impact
in the morphological and thermodynamical properties of the discs. The
comparison with general relativity is further investigated by computing the
size of the photon ring produced by a source located at infinity. This allows
us to place constraints on the parameters of the $f(R)$-gravity model based on
the Event Horizon Telescope observations of the size of the light ring in M87
and SgrA$^*$.
