  We study the strong rate of convergence of the Euler--Maruyama scheme for a
multidimensional stochastic differential equation (SDE) $$ dX_t = b(X_t) \, dt
+ dL_t, $$ with irregular $\beta$-H\"older drift, $\beta > 0$, driven by a
L\'evy process with exponent $\alpha \in (0, 2]$. For $\alpha \in [2/3, 2]$, we
obtain strong $L_p$ and almost sure convergence rates in the entire range
$\beta > 1 - \alpha/2$, where the SDE is known to be strongly well-posed. This
significantly improves the current state of the art, both in terms of
convergence rate and the range of $\alpha$. Notably, the obtained convergence
rate does not depend on $p$, which is a novelty even in the case of smooth
drifts. As a corollary of the obtained moment-independent error rate, we show
that the Euler--Maruyama scheme for such SDEs converges almost surely and
obtain an explicit convergence rate. Additionally, as a byproduct of our
results, we derive strong $L_p$ convergence rates for approximations of
nonsmooth additive functionals of a L\'evy process. Our technique is based on a
new extension of stochastic sewing arguments and L\^e's quantitative
John-Nirenberg inequality.
