Let me start by simplifying your inputs. I will show only inputs, not outputs, except for one at the very end.

1. This is an abstract tensor computation. You only need to load xTensor:

   <<xAct`xTensor`
   

2. Configure signs as soon as possible:

   $PrePrint = ScreenDollarIndices;
   $RiemannSign = -1;
   

3. Construct your background manifold:

   DefConstantSymbol[dim]
   DefManifold[M, dim, IndexRange[a, e]]
   DefMetric[-1, g[-a, -b], lc]
   

4. Declare the new derivative. You don't need to declare the torsion separately. You can get the torsion of cd with Torsion[cd].

   DefCovD[cd[-a], Torsion -> True, {"|", "\[ScriptCapitalD]"}]
   DefTensor[f[-a], M]
   

5. Now set S. You don't need to declare it as a tensor separately. This definition tells us all we need to know about it:

   IndexSet[S[a_, -b_, -c_], g[a, -b] f[-c] + g[a, -c] f[-b] - g[-b, -c] f[a]]
   

Now that we have all objects, let us get the curvature tensors.

6. Riemann. I can imagine you are interested in this difference, where I already expand the Riemann tensors into derivatives of Christoffels:

   Riemanncd[-a, -b, -c, d] - Riemannlc[-a, -b, -c, d] // RiemannToChristoffel
   

   Then you need to impose your rule relating the two Christoffels:

   % /. Christoffelcd[a_, -b_, -c_] :> Christoffellc[a, -b, -c] + S[a, -b, -c] // Expand
   

   And simplify a bit:

   % // ContractMetric // Simplification
   

7. Ricci. We just need to contract two indices:

   % delta[-d, b] // ContractMetric // Simplification
   

8. Ricci scalar. We need to contract with the metric. I will use the background metric in both cases. Not sure this is what you want/need:

   % g[a, c] // ContractMetric // Simplification
   

   In this I'll also expand the Christoffels of the background metric into derivatives of that metric, and finally simplify:

   % // ChristoffelToMetric // Simplification
   

   This final result is simple enough that can be copied here. Recall this is the difference of Ricci scalars:

   Out[]= (-1 + dim) ((-2 + dim) f[-a] f[a] + PD[-a][f[a]] + f[a] g[c, d] PD[-a][g[-c, -d]] + g[-a, -c] PD[c][f[a]])
   

A previous answer I wrote fully missed the point. It helped to see the code in Wolfram Community.

The problem is that in the call

SetOptions[ToCanonical, UseMetricOnVBundle -> g];

you have a metric on the RHS of the rule, while this expects a vbundle, in your case TangentM. This seems to break the canonicalization algorithm.