Let' s calculate one-loop counterterms of minimal second order operator using methods described in Calculating one-loop counterterms. The operator is: \[ D_i{}^j\,=\,\delta_i{}^j\,\Box + W_i{}^j \]
The input quantities needed for the algorithm described in Calculating one-loop counterterms are: \begin{eqnarray*} &&K^{\mu\nu}{}_\alpha{}^\beta \,=\, g^{\mu\nu} \delta_\alpha^\beta, \quad S^\mu{}_\alpha{}^\beta \,=\, 0\,, \quad W_\alpha{}^\beta \,=\, W_\alpha{}^\beta\,, \\&& (Kn)^{-1}{}_{\alpha}{}^\beta \,=\, \delta_{\alpha}^\beta\,, \quad F_{\mu\nu}{}_\alpha{}^\beta \,=\, F_{\mu\nu}{}_\alpha{}^\beta\,. \end{eqnarray*}
The following code produces the result:
//setup symmetries of Riemann tensor
addSymmetries 'R_abcd', -[[0, 1]].p, [[0, 2], [1, 3]].p
setSymmetric 'R_ab'
//(Kn)^(-1)
def iK = 'iK_a^b = d^b_a'.t
def K = 'K^lm_a^b = d^b_a*g^{lm}'.t
def S = 'S^lab = 0'.t
def W = 'W_a^b = W_a^b'.t
def F = 'F_lmab = F_lmab'.t
def div = oneloopdiv2(iK, K, S, W, F)
def counterterms = EliminateDueSymmetries >> div.counterterms
println counterterms
> counterterms = (1/30)*R**2+(1/12)*F_{bm}^{c}_{d}*F^{bmd}_{c}
+(1/2)*W^{a_{5}}_{a}*W^{a}_{a_{5}}
+(1/6)*R*W^{a_{5}}_{a_{5}}+(1/15)*R_{ad}*R^{ad}
Multiplying the produced result by $1\left/16\pi(d-4)\right.$ and integrating over the space-time volume gives:
\[ \Gamma^{(1)}_{\infty} = \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{1}{30} R^2 +\frac{1}{12} F_{\nu\beta }{}^{\epsilon}{}_{\rho} F^{\nu\beta\rho}{}_{\epsilon}+\frac{1}{15} R_{\delta \nu } R^{\delta\nu }+\frac{1}{2} W^{\alpha }{}_{\rho} W^{\rho}_{\alpha}+\frac{1}{6} R W^{\beta }{}_{\beta}\right), \]
where $F_{\mu\nu\alpha\beta}$ is a curvature tensor with respect to the principal bundle, $R_{\mu\nu}$ is a Ricci tensor and $R$ is a Riemann scalar curvature.