term. The second term can thus be evaluated in 4 dimensions with numerical methods. The square bracket in the first term can be evaluated analytically once and for all. One usually defines $$V_{{\ensuremath{\operatorname{sv}}}} (\Phi_B) = \lim_{\epsilon \rightarrow 0} \left[ V (\Phi_B) +
  \int {\mathrm{d}}\Phi_{{\ensuremath{\operatorname{rad}}}} R_s (\Phi_B, \Phi_{{\ensuremath{\operatorname{rad}}}}) \right],$$ and the first term becomes the 4-dimensional expression $$\int {\mathrm{d}}\Phi_B [B (\Phi_B) + V_{{\ensuremath{\operatorname{sv}}}} (\Phi_B)]  \mathcal{O}
  (\Phi_B),$$ that can be computed numerically.

The development of the subtraction method started since the very early QCD computations, already appearing in the bud in the calculation of the Drell-Yan process of ref. [@Altarelli:1979ub]. A more systematic use of it was made in ref. [[@Ellis:1980wv]]{}, in the context of $e^+ e^-$ annihilation into hadrons. In ref. [[@Kunszt:1989km]]{} the calculation of ref. [[@Ellis:1980wv]]{} was implemented as a parton level generator, such that any given observable could be computed with it without any dedicated analytic work, and was in fact used to compute a number of commonly used IR-safe observables for QCD studies at LEP. Subsequently, the subtraction method implemented in parton level generators was applied also for processes initiated by hadrons [[@Mele:1990bq]]{}, and it became common practice to compute the $R_s$ term by using the collinear and 