I think the simplest answer is just to understand that Black and Scholes uses risk neutral probabilities, vis-a-vis "d1" and "d2" in commonly cited derivations. These are very similar to "Z-scores" in statistics, but they are derived with the assumption that a dynamically risk-free portfolio can be constructed -- this allows us to derive a closed form solution for the differential equation. Just as how a standard Z-score can be converted to a probability, the Black and Scholes converts d1 and d2 to a "risk neutral" probabilities through the cumulative distribution function (CDF) of the standard normal ("i.e., Gaussian) variety.
If you want to find the raw probability of an options expiring above or below a given strike, it would be so much easier to just utilize the Z-score formula, as in:
$$Z = \frac{\ln(S / K) - m}{s \sqrt{T-t}}$$$$Z = \frac{\ln(S / K) - m(T-t)}{s \sqrt{T-t}}$$
$$p = \phi(Z)$$
where $S$ is the stock price, $K$ is the arbitrary strike price, $m$ is expected annualized drift, $s$ is the annualized volatility, and $T - t$ is the time till expiration in years, $p$ is the probability that the underlying will expire at or above the strike at time $T - t$, and $\phi$ is the CDF of the Gaussian distribution.