This analysis aims to conduct a comprehensive study of the decay kinematics of the semileptonic decay \signal, with $D_s^{*-} \to D_s^- \gamma$ and $D_s^- \to K^+ K^- \pi^-$, using data collected by the LHCb experiment in Run 2. A first measurement of the form factors describing the $B_s^0$ meson semileptonic decay is provided, performing a four-dimensional binned fit in the space given by the variables describing the decay kinematics, namely $q^2$, $\cos\theta_{\ell}$, $\cos\theta_d$ and $\chi$. Taking into account the detector acceptance, as well as the reconstruction efficiencies and the resolution effects, the full differential distribution is obtained; then, a fit to this distribution is performed using different parameterisations for the $B_s^0 \rightarrow D_s^*$ transition form factors. Furthermore, the unfolded distributions are compared with the theoretical predictions and the Belle-II experiment results. Finally, a model-independent approach is tested and its compatibility with the model-dependent results is studied.