This example uses a custom QCD model created with FeynRules. Please evaluate the file FeynCalc/Examples/FeynRules/QCD/GenerateModelQCD.m before running it for the first time.
description = "Renormalization, massless QCD, MS and MSbar, 1-loop";
If[ $FrontEnd === Null,
$FeynCalcStartupMessages = False;
Print[description];
];
If[ $Notebooks === False,
$FeynCalcStartupMessages = False
];
LaunchKernels[4];
$LoadAddOns = {"FeynArts", "FeynHelpers"};
<< FeynCalc`
$FAVerbose = 0;
$ParallelizeFeynCalc = True;
FCCheckVersion[10, 2, 0];
If[ToExpression[StringSplit[$FeynHelpersVersion, "."]][[1]] < 2,
Print["You need at least FeynHelpers 2.0 to run this example."];
Abort[];
]\text{FeynCalc }\;\text{10.2.0 (dev version, 2025-12-22 21:09:03 +01:00, fcd53f9b). For help, use the }\underline{\text{online} \;\text{documentation},}\;\text{ visit the }\underline{\text{forum}}\;\text{ and have a look at the supplied }\underline{\text{examples}.}\;\text{ The PDF-version of the manual can be downloaded }\underline{\text{here}.}
\text{If you use FeynCalc in your research, please evaluate FeynCalcHowToCite[] to learn how to cite this software.}
\text{Please keep in mind that the proper academic attribution of our work is crucial to ensure the future development of this package!}
\text{FeynArts }\;\text{3.12 (27 Mar 2025) patched for use with FeynCalc, for documentation see the }\underline{\text{manual}}\;\text{ or visit }\underline{\text{www}.\text{feynarts}.\text{de}.}
\text{If you use FeynArts in your research, please cite}
\text{ $\bullet $ T. Hahn, Comput. Phys. Commun., 140, 418-431, 2001, arXiv:hep-ph/0012260}
\text{FeynHelpers }\;\text{2.0.0 (2025-12-22 19:07:44 +01:00, c92fb9f5). For help, use the }\underline{\text{online} \;\text{documentation},}\;\text{ visit the }\underline{\text{forum}}\;\text{ and have a look at the supplied }\underline{\text{examples}.}\;\text{The PDF-version of the manual can be downloaded }\underline{\text{here}.}
\text{ If you use FeynHelpers in your research, please evaluate FeynHelpersHowToCite[] to learn how to cite this work.}
We keep scaleless B0 functions, since otherwise the UV part would not come out right.
$KeepLogDivergentScalelessIntegrals = True;FAPatch[PatchModelsOnly -> True];
(*Successfully patched FeynArts.*)Nicer typesetting
FCAttachTypesettingRule[mu, "\[Mu]"];
FCAttachTypesettingRule[nu, "\[Nu]"];
FCAttachTypesettingRule[rho, "\[Rho]"];
FCAttachTypesettingRule[si, "\[Sigma]"];params = {InsertionLevel -> {Particles}, Model -> FileNameJoin[{"QCD", "QCD"}],
GenericModel -> FileNameJoin[{"QCD", "QCD"}], ExcludeParticles -> {F[3 | 4, {2 | 3}], F[4, {1}]}};
top[i_, j_] := CreateTopologies[1, i -> j,
ExcludeTopologies -> {Tadpoles, WFCorrections, WFCorrectionCTs}];
topTriangle[i_, j_] := CreateTopologies[1, i -> j,
ExcludeTopologies -> {Tadpoles, WFCorrections, WFCorrectionCTs, SelfEnergies}];
topCT[i_, j_] := CreateCTTopologies[1, i -> j,
ExcludeTopologies -> {Tadpoles, WFCorrections, WFCorrectionCTs}];
topTriangleCT[i_, j_] := CreateCTTopologies[1, i -> j,
ExcludeTopologies -> {Tadpoles, WFCorrections, WFCorrectionCTs, SelfEnergyCTs}];
{diagQuarkSE, diagQuarkSECT} = InsertFields[#, {F[3, {1}]} -> {F[3, {1}]},
Sequence @@ params] & /@ {top[1, 1], topCT[1, 1]};
{diagGluonSE, diagGluonSECT} = InsertFields[#, {V[5]} -> {V[5]},
Sequence @@ params] & /@ {top[1, 1], topCT[1, 1]};
{diagGhostSE, diagGhostSECT} = InsertFields[#, {U[5]} -> {U[5]},
Sequence @@ params] & /@ {top[1, 1], topCT[1, 1]};
{diagQuarkGluonVertex, diagQuarkGluonVertexCT} = InsertFields[#,
{F[3, {1}], V[5]} -> {F[3, {1}]}, Sequence @@ params] & /@ {topTriangle[2, 1], topTriangleCT[2, 1]};diag1[0] = diagQuarkSE[[0]][Sequence @@ diagQuarkSE,
Sequence @@ diagQuarkSECT];
diag2[0] = diagGluonSE[[0]][Sequence @@ diagGluonSE,
Sequence @@ diagGluonSECT];
diag3[0] = diagGhostSE[[0]][Sequence @@ diagGhostSE,
Sequence @@ diagGhostSECT];
diag4[0] = diagQuarkGluonVertex[[0]][Sequence @@ diagQuarkGluonVertex,
Sequence @@ diagQuarkGluonVertexCT];Paint[diag1[0], ColumnsXRows -> {2, 1}, SheetHeader -> None,
Numbering -> Simple, ImageSize -> {512, 256}];
Paint[diag2[0], ColumnsXRows -> {4, 1}, SheetHeader -> None,
Numbering -> Simple, ImageSize -> {512, 128}];
Paint[diag3[0], ColumnsXRows -> {2, 1}, SheetHeader -> None,
Numbering -> Simple, ImageSize -> {512, 256}];
Paint[diag4[0], ColumnsXRows -> {3, 1}, SheetHeader -> None,
Numbering -> Simple, ImageSize -> {512, 256}];
The 1/(2Pi)^D prefactor is implicit.
Quark self-energy including the counter-term
ampQuarkSE[0] = FCFAConvert[CreateFeynAmp[diag1[0], Truncated -> True,
GaugeRules -> {}, PreFactor -> 1],
IncomingMomenta -> {p}, OutgoingMomenta -> {p},
LorentzIndexNames -> {mu, nu}, DropSumOver -> True,
LoopMomenta -> {l}, UndoChiralSplittings -> True,
ChangeDimension -> D, SMP -> True,
FinalSubstitutions -> {Zm -> SMP["Z_m"], Zpsi -> SMP["Z_psi"],
SMP["m_u"] -> 0}];Gluon self-energy including the counter-term
ampGluonSE[0] = FCFAConvert[CreateFeynAmp[diag2[0], Truncated -> True,
GaugeRules -> {}, PreFactor -> 1],
IncomingMomenta -> {p}, OutgoingMomenta -> {p},
LorentzIndexNames -> {mu, nu, rho, si}, DropSumOver -> True,
LoopMomenta -> {l}, UndoChiralSplittings -> True,
ChangeDimension -> D, SMP -> True,
FinalSubstitutions -> {ZA -> SMP["Z_A"], Zxi -> SMP["Z_xi"],
SMP["m_u"] -> 0}];Ghost self-energy including the counter-term
ampGhostSE[0] = FCFAConvert[CreateFeynAmp[diag3[0], Truncated -> True,
GaugeRules -> {}, PreFactor -> 1],
IncomingMomenta -> {p}, OutgoingMomenta -> {p},
LorentzIndexNames -> {mu, nu}, DropSumOver -> True,
LoopMomenta -> {l}, UndoChiralSplittings -> True,
ChangeDimension -> D, SMP -> True,
FinalSubstitutions -> {Zu -> SMP["Z_u"]}];Quark-gluon vertex including the counter-term
ampQGlVertex[0] = FCFAConvert[CreateFeynAmp[diag4[0], Truncated -> True,
GaugeRules -> {}, PreFactor -> 1],
IncomingMomenta -> {p1, k}, OutgoingMomenta -> {p2},
LorentzIndexNames -> {mu, nu, rho}, DropSumOver -> True, LoopMomenta -> {l},
UndoChiralSplittings -> True, ChangeDimension -> D,
SMP -> True, FinalSubstitutions ->
{ZA -> SMP["Z_A"], Zg -> SMP["Z_g"], Zpsi -> SMP["Z_psi"],
SMP["m_u"] -> 0}];Tensor reduction allows us to express the quark self-energy in tems of the Passarino-Veltman coefficient functions.
ampQuarkSE[1] = ampQuarkSE[0] // SUNSimplify[#, FCParallelize -> True] & //
DiracSimplify[#, FCParallelize -> True] & // TID[#, l, UsePaVeBasis -> True,
ToPaVe -> True, FCParallelize -> True] & // Total;Discard all the finite pieces of the 1-loop amplitude.
ampQuarkSEDiv[0] = ampQuarkSE[1] // PaVeUVPart[#, Prefactor -> 1/(2 Pi)^D] &;ampQuarkSEDiv[1] = FCReplaceD[ampQuarkSEDiv[0], D -> 4 - 2 Epsilon] //
Series[#, {Epsilon, 0, 0}] & // Normal // FCHideEpsilon // Simplify\frac{i \delta _{\text{Col1}\;\text{Col2}} \gamma \cdot p \left(C_F (\Delta +\gamma +\log (\pi )) \xi _{\text{G}} g_s^2+16 \pi ^2 \left(Z_{\psi }-1\right)\right)}{16 \pi ^2}
ampQuarkSEDiv[2] = ampQuarkSEDiv[1] // ReplaceRepeated[#, {
SMP["Z_m"] -> 1 + alpha SMP["d_m"],
SMP["Z_psi"] -> 1 + alpha SMP["d_psi"]}] & // Series[#, {alpha, 0, 1}] & //
Normal // ReplaceAll[#, alpha -> 1] & // SelectNotFree2[#, SMP["Delta"], SMP["d_m"],
SMP["d_psi"]] &\frac{i \Delta C_F \xi _{\text{G}} g_s^2 \delta _{\text{Col1}\;\text{Col2}} \gamma \cdot p}{16 \pi ^2}+i \delta _{\psi } \delta _{\text{Col1}\;\text{Col2}} \gamma \cdot p
ampQuarkSEDiv[3] = ampQuarkSEDiv[2] // SUNSimplify //
Collect2[#, DiracGamma, Factoring -> Simplify] &\frac{i \left(C_A-2 C_F\right) \delta _{\text{Col1}\;\text{Col2}} \gamma \cdot p \left(32 \pi ^2 C_A \delta _{\psi }+\Delta \left(C_A^2-1\right) \xi _{\text{G}} g_s^2\right)}{32 \pi ^2}
sol[1] = Solve[ampQuarkSEDiv[3] == 0, SMP["d_psi"]] //
Flatten // ReplaceAll[#, Rule[a_, b_] :> Rule[a, SUNSimplify[b]]] & //
ReplaceAll[#, SMP["g_s"]^2 -> 4 Pi SMP["alpha_s"]] &;
solMS1 = sol[1] /. {
SMP["d_psi"] -> SMP["d_psi^MS"], SMP["Delta"] -> 1/Epsilon
}
solMSbar1 = sol[1] /. {
SMP["d_psi"] -> SMP["d_psi^MSbar"]
}\left\{\delta _{\psi }^{\text{MS}}\to -\frac{C_F \xi _{\text{G}} \alpha _s}{4 \pi \varepsilon }\right\}
\left\{\delta _{\psi }^{\overset{---}{\text{MS}}}\to -\frac{\Delta C_F \xi _{\text{G}} \alpha _s}{4 \pi }\right\}
Tensor reduction allows us to express the gluon self-energy in tems of the Passarino-Veltman coefficient functions.
ampGluonSE[1] = ({ampGluonSE[0][[1]], Nf ampGluonSE[0][[2]],
Sequence @@ ampGluonSE[0][[3 ;;]]}) // SUNSimplify[#, FCParallelize -> True] & //
DiracSimplify[#, FCParallelize -> True] &;ampGluonSE[2] = TID[ampGluonSE[1], l, UsePaVeBasis -> True, ToPaVe -> True, FCParallelize -> True];Discard all the finite pieces of the 1-loop amplitude
ampGluonSEDiv[0] = ampGluonSE[2] // Total // PaVeUVPart[#, Prefactor -> 1/(2 Pi)^D] &;ampGluonSEDiv[1] = FCReplaceD[ampGluonSEDiv[0], D -> 4 - 2 Epsilon] //
Series[#, {Epsilon, 0, 0}] & //
Normal // FCHideEpsilon // SUNSimplify;ampGluonSEDiv[2] = ampGluonSEDiv[1] // ReplaceRepeated[#, {
SMP["Z_A"] -> 1 + alpha SMP["d_A"],
SMP["Z_xi"] -> 1 + alpha SMP["d_A"]}] & // Series[#, {alpha, 0, 1}] & //
Normal // ReplaceAll[#, alpha -> 1] & // SelectNotFree2[#, SMP["Delta"], SMP["d_A"],
SMP["d_xi"]] &\frac{i \Delta C_A \xi _{\text{G}} g_s^2 p^{\mu } p^{\nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{32 \pi ^2}-\frac{i \Delta p^2 C_A \xi _{\text{G}} g_s^2 g^{\mu \nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{32 \pi ^2}-\frac{13 i \Delta C_A g_s^2 p^{\mu } p^{\nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{96 \pi ^2}+\frac{13 i \Delta p^2 C_A g_s^2 g^{\mu \nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{96 \pi ^2}-i p^2 \delta _A g^{\mu \nu } \delta ^{\text{Glu1}\;\text{Glu2}}+i \delta _A p^{\mu } p^{\nu } \delta ^{\text{Glu1}\;\text{Glu2}}+\frac{i \Delta N_f g_s^2 p^{\mu } p^{\nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{24 \pi ^2}-\frac{i \Delta p^2 N_f g_s^2 g^{\mu \nu } \delta ^{\text{Glu1}\;\text{Glu2}}}{24 \pi ^2}
sol[3] = Solve[ampGluonSEDiv[2] == 0, SMP["d_A"]] // Flatten //
ReplaceAll[#, SMP["g_s"]^2 -> 4 Pi SMP["alpha_s"]] & // Simplify
solMS2 = sol[3] /. {SMP["d_A"] -> SMP["d_A^MS"], SMP["Delta"] -> 1/Epsilon}
solMSbar2 = sol[3] /. {SMP["d_A"] -> SMP["d_A^MSbar"]}\left\{\delta _A\to -\frac{\Delta \alpha _s \left(3 C_A \xi _{\text{G}}-13 C_A+4 N_f\right)}{24 \pi }\right\}
\left\{\delta _A^{\text{MS}}\to -\frac{\alpha _s \left(3 C_A \xi _{\text{G}}-13 C_A+4 N_f\right)}{24 \pi \varepsilon }\right\}
\left\{\delta _A^{\overset{---}{\text{MS}}}\to -\frac{\Delta \alpha _s \left(3 C_A \xi _{\text{G}}-13 C_A+4 N_f\right)}{24 \pi }\right\}
Tensor reduction allows us to express the ghost self-energy in tems of the Passarino-Veltman coefficient functions.
ampGhostSE[1] = ampGhostSE[0] // SUNSimplify[#, FCParallelize -> True] & //
DiracSimplify[#, FCParallelize -> True] &;ampGhostSE[2] = TID[ampGhostSE[1], l, UsePaVeBasis -> True, ToPaVe -> True, FCParallelize -> True];Discard all the finite pieces of the 1-loop amplitude
ampGhostSEDiv[0] = ampGhostSE[2] // Total // PaVeUVPart[#, Prefactor -> 1/(2 Pi)^D] &\frac{i 2^{-D-1} \pi ^{-D} p^2 \delta ^{\text{Glu1}\;\text{Glu2}} \left(\pi ^2 \left(-C_A\right) \xi _{\text{G}} g_s^2+3 \pi ^2 C_A g_s^2-2^{D+3} \pi ^D Z_u+2^{D+1} D \pi ^D Z_u+2^{D+3} \pi ^D-2^{D+1} D \pi ^D\right)}{D-4}
ampGhostSEDiv[1] = FCReplaceD[ampGhostSEDiv[0], D -> 4 - 2 Epsilon] //
Series[#, {Epsilon, 0, 0}] & //
Normal // FCHideEpsilon // SUNSimplify-\frac{1}{64 \pi ^2}i p^2 \delta ^{\text{Glu1}\;\text{Glu2}} \left(-\Delta C_A \xi _{\text{G}} g_s^2+\gamma \left(-C_A\right) \xi _{\text{G}} g_s^2+\log (4 \pi ) C_A \xi _{\text{G}} g_s^2-2 \log (\pi ) C_A \xi _{\text{G}} g_s^2-2 \log (2) C_A \xi _{\text{G}} g_s^2+3 \Delta C_A g_s^2+3 \gamma C_A g_s^2-3 \log (4 \pi ) C_A g_s^2+6 \log (\pi ) C_A g_s^2+6 \log (2) C_A g_s^2-64 \pi ^2 Z_u+64 \pi ^2\right)
ampGhostSEDiv[2] = ampGhostSEDiv[1] // ReplaceRepeated[#, {
SMP["Z_u"] -> 1 + alpha SMP["d_u"]}] & // Series[#, {alpha, 0, 1}] & //
Normal // ReplaceAll[#, alpha -> 1] & //
SelectNotFree2[#, SMP["Delta"], SMP["d_u"]] & // Simplify\frac{i p^2 \delta ^{\text{Glu1}\;\text{Glu2}} \left(\Delta C_A \left(\xi _{\text{G}}-3\right) g_s^2+64 \pi ^2 \delta _u\right)}{64 \pi ^2}
sol[4] = Solve[ampGhostSEDiv[2] == 0, SMP["d_u"]] // Flatten //
ReplaceAll[#, SMP["g_s"]^2 -> 4 Pi SMP["alpha_s"]] & // Simplify
solMS3 = sol[4] /. {SMP["d_u"] -> SMP["d_u^MS"], SMP["Delta"] -> 1/Epsilon}
solMSbar3 = sol[4] /. {SMP["d_u"] -> SMP["d_u^MSbar"]}\left\{\delta _u\to -\frac{\Delta C_A \left(\xi _{\text{G}}-3\right) \alpha _s}{16 \pi }\right\}
\left\{\delta _u^{\text{MS}}\to -\frac{C_A \left(\xi _{\text{G}}-3\right) \alpha _s}{16 \pi \varepsilon }\right\}
\left\{\delta _u^{\overset{---}{\text{MS}}}\to -\frac{\Delta C_A \left(\xi _{\text{G}}-3\right) \alpha _s}{16 \pi }\right\}
Tensor reduction allows us to express the quark-gluon vertex in tems of the Passarino-Veltman coefficient functions.
ampQGlVertex[1] = ampQGlVertex[0] // SUNSimplify[#, FCParallelize -> True] & //
DiracSimplify[#, FCParallelize -> True] &;ampQGlVertex[2] = TID[ampQGlVertex[1], l, UsePaVeBasis -> True, ToPaVe -> True, FCParallelize -> True];Discard all the finite pieces of the 1-loop amplitude
ampQGlVertexDiv[0] = ampQGlVertex[2] // Total // PaVeUVPart[#, Prefactor -> 1/(2 Pi)^D] &;ampQGlVertexDiv[1] = FCReplaceD[ampQGlVertexDiv[0], D -> 4 - 2 Epsilon] //
Series[#, {Epsilon, 0, 0}] & //
Normal // FCHideEpsilon // SUNSimplify;ampQGlVertexDiv[2] = ampQGlVertexDiv[1] // ReplaceRepeated[#, {
SMP["Z_g"] -> 1 + alpha SMP["d_g"],
SMP["Z_A"] -> 1 + alpha SMP["d_A"],
SMP["Z_psi"] -> 1 + alpha SMP["d_psi"]
}] & // Series[#, {alpha, 0, 1}] & //
Normal // ReplaceAll[#, alpha -> 1] & //
SelectNotFree2[#, SMP["Delta"], SMP["d_g"], SMP["d_A"], SMP["d_psi"]] & // Simplify-\frac{i \gamma ^{\mu } g_s \left(C_A-2 C_F\right) T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}} \left(64 \pi ^2 C_A \delta _{\psi }+32 \pi ^2 C_A \delta _A+3 \Delta C_A^2 \xi _{\text{G}} g_s^2+64 \pi ^2 C_A \delta _g+3 \Delta C_A^2 g_s^2-2 \Delta \xi _{\text{G}} g_s^2\right)}{64 \pi ^2}
ampQGlVertexDiv[3] = ampQGlVertexDiv[2] // Collect2[#, SMP] &-\frac{i \Delta \gamma ^{\mu } g_s^3 \left(C_A-2 C_F\right) \left(3 C_A^2 \xi _{\text{G}}+3 C_A^2-2 \xi _{\text{G}}\right) T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}}}{64 \pi ^2}-i C_A \gamma ^{\mu } \delta _{\psi } g_s \left(C_A-2 C_F\right) T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}}-\frac{1}{2} i C_A \delta _A \gamma ^{\mu } g_s \left(C_A-2 C_F\right) T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}}-i C_A \gamma ^{\mu } \delta _g g_s \left(C_A-2 C_F\right) T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}}
ampQGlVertexDiv[4] = (ampQGlVertexDiv[3] /. {
SMP["d_A"] -> SMP["d_A^MS"],
SMP["d_psi"] -> SMP["d_psi^MS"],
SMP["d_g"] -> SMP["d_g"],
SMP["Delta"] -> 1/Epsilon
} /. solMS1 /. solMS2) // ReplaceAll[#, SMP["g_s"]^3 -> 4 Pi SMP["alpha_s"] SMP["g_s"]] & //
Collect2[#, Epsilon] & // SUNSimplify\frac{i \gamma ^{\mu } g_s T_{\text{Col3}\;\text{Col1}}^{\text{Glu2}} \left(-11 C_A \alpha _s+2 N_f \alpha _s-24 \pi \varepsilon \delta _g\right)}{24 \pi \varepsilon }
sol[5] = Solve[ampQGlVertexDiv[4] == 0, SMP["d_g"]] // Flatten // Simplify
solMS4 = sol[5] /. {SMP["d_g"] -> SMP["d_g^MS"]}
solMSbar4 = sol[5] /. {SMP["d_g"] -> SMP["d_g^MSbar"], 1/Epsilon -> SMP["Delta"]}\left\{\delta _g\to \frac{\alpha _s \left(2 N_f-11 C_A\right)}{24 \pi \varepsilon }\right\}
\left\{\delta _g^{\text{MS}}\to \frac{\alpha _s \left(2 N_f-11 C_A\right)}{24 \pi \varepsilon }\right\}
\left\{\delta _g^{\overset{---}{\text{MS}}}\to \frac{\Delta \alpha _s \left(2 N_f-11 C_A\right)}{24 \pi }\right\}
```mathematica knownResult = { SMP[“d_psi^MS”] -> -SMP[“alpha_s”]/(4 Pi) 1/Epsilon CF GaugeXi[“G”], SMP[“d_psi^MSbar”] -> -SMP[“alpha_s”]/(4 Pi) SMP[“Delta”] CF GaugeXi[“G”], SMP[“d_A^MS”] -> SMP[“alpha_s”]/(4 Pi) 1/Epsilon (1/2 CA (13/3 - GaugeXi[“G”]) - 2/3 Nf), SMP[“d_A^MSbar”] -> SMP[“alpha_s”]/(4 Pi) SMP[“Delta”] (1/2 CA (13/3 - GaugeXi[“G”]) - 2/3 Nf), SMP[“d_u^MS”] -> SMP[“alpha_s”]/(4 Pi) CA/Epsilon (3 - GaugeXi[“G”])/4, SMP[“d_u^MSbar”] -> SMP[“alpha_s”]/(4 Pi) CA SMP[“Delta”] (3 - GaugeXi[“G”])/4, SMP[“d_g^MS”] -> ((-11CASMP[“alpha_s”])/(24EpsilonPi) + (NfSMP[“alpha_s”])/(12EpsilonPi)), SMP[“d_g^MSbar”] -> ((-11CASMP[“alpha_s”])/(24 Pi) SMP[“Delta”] + (NfSMP[“alpha_s”])/(12*Pi) SMP[“Delta”]) } // Factor2; FCCompareResults[Join[solMS1, solMSbar1, solMS2, solMSbar2, solMS3, solMSbar3, solMS4,solMSbar4] // Factor2, knownResult, Text -> {“to Muta, Foundations of QCD, Eqs 2.5.131-2.5.147:”, “CORRECT.”, “WRONG!”}, Interrupt -> {Hold[Quit[1]], Automatic}] Print[“Time used:”, Round[N[TimeUsed[], 4], 0.001], ” s.”];
```mathematica
\text{$\backslash $tCompare to Muta, Foundations of QCD, Eqs 2.5.131-2.5.147:} \;\text{CORRECT.}
\text{True}
\text{$\backslash $tCPU Time used: }28.323\text{ s.}