xff[g, mx, m[Psi]] = Quiet[Table[Re[x /. FindRoot[x - Log[(0.038*2*1.22*10^19*m[Psi]*sigv[x, g, mx, m[Psi]])/
     Sqrt[107*x]], {x, 25}]] /. {g -> matt[[i, 1]],
 mx -> matt[[i, 2]], m[Psi] -> matt[[i, 3]]}, {i, 1, 10000000}]]

replace1 = {me -> 0.0, mμ -> 0.0, mτ -> 0.0, md -> 0.0,
   ms -> 0.0, mb -> 0.0, mu -> 0.0, mc -> 0.0, mt -> 173.21,
   mW -> 80.385, mZ -> 91.188, vev -> 246, mh -> 125};
replace = {me -> 0.000511, mμ -> 0.1057, mτ -> 1.777,
   md -> 0.0048, ms -> 0.095, mb -> 4.18, mu -> 0.0023, mc -> 1.275,
   mt -> 173.21, mW -> 80.385, mZ -> 91.188, vev -> 246, mh -> 125};

rulel = {nc -> 1, gfv -> 0, gfa -> -0.707 g};
ruleU = {nc -> 3, gfv -> 0.884 g, gfa -> 0};
ruleD = {nc -> 3, gfv -> 0.707 g, gfa -> 0};
ruleN = {nc -> 1, gfv -> 0.354 g, gfa -> -0.354 g};
rule1 = {gψa -> 1.7536718949856565` g, gψv -> 0};
rule2 = {gψa -> 1.7536718949856565` g, gψv -> 0};

n1 = n4 = 1/Sqrt[2] // N;
θ = 1;
g1 = 5/8 (3 n1 + n4) g;
g2 = 1/8 (3 n1 + n4) g;

g11 = g1*Cos[θ]^2 + g2*Sin[θ]^2;
g22 = g2*Cos[θ]^2 + g1*Sin[θ]^2;
g12 = g21 = Sin[2 θ]/2 (g1 - g2);

geff = 2 + 2 (1 + Δ)^(3/2) Exp[-x (Δ)];

(* Total Decay width of X: *)
Γtt1[g_, mx_] = ((3*nc*mx)/(24*π) (1)^(1/2) (gfa^2 (1) + gfv^2 (1)) /. rulel) +
                ((3 nc*mx)/(24*π) (1)^(1/2) (gfa^2 (1) + gfv^2 (1)) /. ruleN) +
                ((3*nc*mx)/(24*π) (1)^(1/2) (gfa^2 (1) + gfv^2 (1)) /. ruleU) +
                ((3*nc*mx)/(24*π) (1)^(1/2) (gfa^2 (1) + gfv^2 (1)) /. ruleD);

Γfd = 0.1828052084804079` g^2 mx;

Δ = 1;

(* Thermal cross section via integration formula: *)

sigma[x_, g_, mx_, mψ_, mf_] =
   3/(12*π*s ((s - mx^2)^2 + mx^2*(Γfd)^2)) ((1 - 4*mf^2/s)/(1 - 4*mψ^2/s))^(1/2)*
   (gfa^2 g11^2 (4 mψ^2 (mf^2 (7 - 6 s/mx^2 + 3 s^2/mx^4) - s) +
   s (s - 4 mf^2)) + gfv^2 g11^2 (s + 2 mf^2) (s - 4 mψ^2)); /. mf -> 0

sig[x_, g_, mx_, mψ_] = (sigma[x, g, mx, mψ, me] /. rulel /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, mμ] /. rulel /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, mτ] /. rulel /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, md] /. ruleD /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, md] /. ruleD /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, ms] /. ruleD /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, mu] /. ruleU /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, mc] /. ruleU /. replace /. rule1) +
                        (sigma[x, g, mx, mψ, mt] /. ruleU /. replace /. rule1);

sig1[x_, g_, mx_, mψ_] = 4/(geff)^2 (1 + 2 g12^2/g11^2 (1 + Δ)^(3/2) Exp[-x (Δ)] +
  g22^2/g11^2 (1 + Δ)^3 Exp[-2 x (Δ)]) sig[x, g, mx, mψ];

sigv[x_?NumericQ, g_?NumericQ, mx_?NumericQ, mψ_?NumericQ] :=
  x/(8 mψ^5*BesselK[2, x]*BesselK[2, x]) (NIntegrate[
     sig1[x, g, mx, mψ]*s^(1/2)*BesselK[1, s^(1/2)*x/mψ], {s, 4 mψ^2, Infinity},
     Method -> {Automatic, "SymbolicProcessing" -> 0},
     PrecisionGoal -> 2]);

(* CALCULATION OF XF: *)

matt = Tuples[{Range[0.005, 0.7, 0.02], Range[100, 5100, 500],
    Range[100, 2600, 400]}];

Clear[xff]
xff[g, mx, mψ] =
   Quiet[Table[
     Re[x /. FindRoot[
         x - Log[(0.038*2*1.22*10^19*mψ*sigv[x, g, mx, mψ])/
            Sqrt[107*x]], {x, 25}]] /. {g -> matt[[i, 1]],
       mx -> matt[[i, 2]], mψ -> matt[[i, 3]]}, {i, 1,
      2695}]]; // AbsoluteTiming

(* {3567.33,Null} *)