First I need to introduce some background which you may already know, and which I touched upon a couple of time in previous answers, so as to make my answer useful to more people. The cross-section for a process such as $pp\to hX$, where $p$ stands for proton, $h$ for some hadron whose momentum is measured, and $X$ for any other particle ignored in the measurement reads very schematically
$$\frac{d\sigma}{d^3p} = \sum_i \sum_j \sum_k \int dx_i f_i(x_i)\int dx_jf_j(x_j)\int dz_kF_k(z_k)|\mathcal{M}(ij\to kX)|^2.$$
There $i$, $j$ and $k$ stand for a species of parton, i.e. quark or gluon. $f_i(x_i)$ is the probability to find in a proton a parton $i$ with a fraction $x_i$ of the proton momentum and it is called for that reason a parton distribution function (PDF). $F_k(z_k)$ is the probability that a parton $k$ will hadronise into a particle $h$ carrying a fraction $z_k$ of the parton momentum and it is therefore called a parton fragmentation function. Finally $\mathcal{M}(ij\to kX)$ is the amplitude for the given partonic process.
I am now in a position to answer your question. The computation of $|\mathcal{M}(ij\to kX)|^2$ will lead to divergences, which come in two flavours:
the ultraviolet (UV) ones, which appear because of large momentum in the loops of the Feynman diagrams representing the amplitude;
the infrared (IR) ones, which appear because (i) either a virtual or a real particle can reach a zero momentum, or because (ii) a massless particle radiates another massless particle.
The UV divergences are cured by introducing the renormalisation scale $\mu_R$, of which the coupling constant $\alpha_S$ becomes a function of. The IR divergences in case (i) cancel out (as predicted by the Kinoshita-Lee-Nauenberg theorem) but not in case (ii): they are cured by introducing the factorisation scale $\mu_F$, of which the parton distribution and fragmentation functions will become a function of. So we end up with a new finite expression
$$\begin{align}
\frac{d\sigma}{d^3p} = \sum_i \sum_j \sum_k \int dx_i f_i(x_i, \mu_F)&\int dx_jf_j(x_j, \mu_F)\int dz_k F_k(z_k,\mu_F)\\
&\times\frac{d\sigma_{ij\to kX}}{d^3p_k}(x_i,x_j,z_k,\mu_F,\mu_R,\alpha_S(\mu_R)).
\end{align}$$
At this point, it is very important to understand that $\mu_F$ and $\mu_R$ are spurious parameters and that physical observable should ideally not depend on them. This would be true if we could sum the entire perturbation series, which is not only practically impossible but theoretically unsound. But at least, the more terms of the series we compute, the less the observable depends on those scales. The worst case is to stay at the leading-order because then the observables are monotonous functions of each of those scales. But even at next-leading order, there is a leftover dependence. Not even next-next-leading order entirely cures the problem. In the face of this issue, there are a couple of traditional tricks:
fix $\mu_F$ and $\mu_R$ all at a supposedly physically meaningful value, for example the one in your question for a particle of transverse momentum $p_T$;
find the values of $\mu_F$ and $\mu_R$ which minimise the observable.
The idea behind the last point is that such a minimum is the point where the observable will vary the least with those scales, and since we know theoretically that they should not vary with them, this is a good as it gets.
Eventually, there is no ideal solution but at least the following three prescriptions shall be strictly followed:
- use as high an order of perturbation theory as possible (and avoid leading-order like the plague);
- make sure that all the computation you use or compare with each other use the same prescription to fix the factorisation and renormalisation scale;
- if you use as scale factors some typical energy of your process ($\sqrt{S}$, or the formula in your question), at least study how the theoretical prediction vary about that point when you vary the scale factors.
The last point will give you an idea of the theoretical uncertainty. It is often enough to compute for $\mu$, $2\mu$ and $\mu/2$ where $\mu$ would be the supposedly physical scale.
The best for the last two points are the theoretical computation which keep $\mu_R$ and $\mu_F$ in the formulae without giving them a special value, leaving it to the user of those formula to make the choice.