Within the widely used EvtGen framework, we have added a new event generator model for $B\to K^* \ell^+ \ell^-$ with improved standard model (SM) decay amplitudes and possible BSM physics contributions, which are implemented in the operator product expansion in terms of Wilson coefficients. This event generator can then be used to investigate the experimental sensitivity to the most general BSM signal resulting from dimension-six operators. We describe the advantages and potential of the newly developed `Sibidanov Physics Generator' in improving experimental sensitivity of searches for lepton non-universal BSM physics and clarifying signatures. The new generator can properly simulate BSM scenarios, interference between SM and BSM amplitudes, and correlations between different BSM observables as well as acceptance bias. We show that exploiting such correlations substantially improves experimental sensitivity. As a demonstration of the utility of the MC generator, we examine the prospects for improved measurements of lepton non-universality in angular distributions for $B \to K^* \ell \ell$ decays from the expected 50 ab$^{-1}$ data set of the Belle II experiment, using a four-dimensional unbinned maximum likelihood fit. We describe promising experimental signatures and correlations between observables. The use of lepton-universality violating $\Delta$-observables significantly reduce uncertainties in the SM expectations due to QCD and resonance effects, are ideally suited for Belle II with the large data sets expected in the next decade. Our simulation studies also show that Belle II should have excellent sensitivity to BSM physics in the Wilson coefficients $C_7$ and $C_7'$, which appears at low $q^2$ in the di-electron channel.