Added examples that will appear in the publication.

Examples/EW/EWHiggsToTwoGluonsOneLoop.m

@@ -0,0 +1,159 @@
+(* ::Package:: *)
+
+(* :Title: EWHiggsToTwoGluonsOneLoop										*)
+
+(*
+	This software is covered by the GNU General Public License 3.
+	Copyright (C) 1990-2016 Rolf Mertig
+	Copyright (C) 1997-2016 Frederik Orellana
+	Copyright (C) 2014-2016 Vladyslav Shtabovenko
+*)
+
+(* :Summary:  Computation of the total decay rate for the decay of a 
+				Higgs into 2 gluons via a top-loop.						*)
+
+(* ------------------------------------------------------------------------ *)
+
+
+
+(* ::Section:: *)
+(*Load FeynCalc and FeynArts*)
+
+
+If[ $FrontEnd === Null,
+		$FeynCalcStartupMessages = False;
+		Print["Computation of the total decay rate for the decay of a Higgs into 2 gluons via a top-loop."];
+];
+If[$Notebooks === False, $FeynCalcStartupMessages = False];
+$LoadAddOns={"FeynHelpers"};
+$LoadFeynArts= True;
+<<FeynCalc`
+$FAVerbose = 0;
+
+
+(* ::Section:: *)
+(*Decay of Higgs into two gluon via a top quark loop*)
+
+
+(* ::Subsection:: *)
+(*Generate Feynman diagrams*)
+
+
+(* ::Text:: *)
+(*Here we consider only the dominant contribution from the top quark mass. However, it is trivial to include also loops from other quark flavors.*)
+
+
+diagsHiggsToGG = InsertFields[CreateTopologies[1,1 -> 2],{S[1]} -> {V[5],V[5]},
+InsertionLevel -> {Particles},Model->"SMQCD"];
+Paint[DiagramExtract[diagsHiggsToGG,{3,6}],ColumnsXRows->{2,1},
+Numbering -> None,SheetHeader->None,ImageSize->{512,256}];
+
+
+(* ::Subsection:: *)
+(*Kinematics*)
+
+
+FCClearScalarProducts[];
+ScalarProduct[k1,k1]=0;
+ScalarProduct[k2,k2]=0;
+ScalarProduct[pH,pH]=SMP["m_H"]^2;
+ScalarProduct[k1,k2]=(SMP["m_H"]^2)/2;
+
+
+(* ::Subsection:: *)
+(*Evaluation of the amplitudes*)
+
+
+(* ::Text:: *)
+(*This is the sum of both amplitudes*)
+
+
+ampHiggsToTwoGluons=FCFAConvert[CreateFeynAmp[DiagramExtract[diagsHiggsToGG,
+{3,6}],PreFactor->-1],IncomingMomenta->{pH},OutgoingMomenta->{k1,k2},
+LoopMomenta->{q},List->False,TransversePolarizationVectors->{k1,k2},
+ChangeDimension->D,DropSumOver->True,SMP->True,
+UndoChiralSplittings->True]//Contract//FCTraceFactor//SUNSimplify
+
+
+(* ::Text:: *)
+(*Tensor reduction is straight-forward, after which we end up with coefficient functions*)
+
+
+ampHiggsToTwoGluons2=TID[ampHiggsToTwoGluons,q,ToPaVe->True]//Collect2[#,B0,C0]&
+
+
+(* ::Text:: *)
+(*The coefficient functions are evaluated with Package-X*)
+
+
+ampHiggsToTwoGluons3=PaXEvaluate[ampHiggsToTwoGluons2,
+PaXImplicitPrefactor->1/(2Pi)^D]//Simplify
+
+
+(* ::Text:: *)
+(*As the result is finite, it is safe to switch from D to 4 dimensions.*)
+
+
+ampHiggsToTwoGluons4=ChangeDimension[ampHiggsToTwoGluons3,4];
+
+
+(* ::Text:: *)
+(*We square the amplitude and sum over the polarizations of the gluons*)
+
+
+ampSq=ampHiggsToTwoGluons4 ComplexConjugate[ampHiggsToTwoGluons4]//
+DoPolarizationSums[#,k1,k2]&//DoPolarizationSums[#,k2,k1]&//Simplify//
+SUNSimplify[#,SUNNToCACF->False]&
+
+
+(* ::Text:: *)
+(*Multiplying the result by the phase space factor we obtain the total decay rate*)
+
+
+$Assumptions={SMP["m_H"]>0,SMP["m_t"]>0};
+decayRateTotal=(1/2*1/(16 Pi SMP["m_H"])*(ampSq/.SUNN->3))//
+ReplaceAll[#,{SMP["e"]^2->4 Pi SMP["alpha_fs"],SMP["g_s"]^4->
+16 Pi^2 SMP["alpha_s"]^2}]&//Simplify
+
+
+(* ::Text:: *)
+(*For convenience, we can eliminate the fine-structure constant and sine of the Weinberg angle in favor of Fermi's constant*)
+
+
+decayRateTotal2=Simplify[decayRateTotal/.{ SMP["m_t"]-> SMP["m_H"]/(2 Sqrt[\[Tau]]),
+SMP["alpha_fs"]->Sqrt[2]/Pi SMP["m_W"]^2 SMP["sin_W"]^2 SMP["G_F"]}]
+
+
+(* ::Text:: *)
+(*To compare to the literature (arXiv:hep-ph/9504378) we extract the function A^2 (tau)*)
+
+
+aSq=decayRateTotal2/(SMP["alpha_s"]^2 SMP["G_F"]SMP["m_H"]^3/(36 Sqrt[2]Pi^3) 9/16)
+
+
+(* ::Text:: *)
+(*The plots show that our results coincide with those of Spira et al.*)
+
+
+aSqLit[x_]:=Piecewise[{{(2(x+(x-1)ArcSin[Sqrt[x]]^2)/x^2)^2,0<x<=1},
+{(2(x+(x-1)(-1/4(Log[(1+\[Sqrt](1-1/x))/(1-\[Sqrt](1-1/x))]- I Pi)^2))/x^2)^2,x>1}}]
+
+
+aSqLitBT=(2(\[Tau]+(\[Tau]-1)ArcSin[Sqrt[\[Tau]]]^2)/\[Tau]^2)^2;
+aSqLitAT=(2(\[Tau]+(\[Tau]-1)(-1/4(Log[(1+\[Sqrt](1-1/\[Tau]))/(1-\[Sqrt](1-1/\[Tau]))]- I Pi)^2))/\[Tau]^2)^2;
+
+
+plot1=Plot[{Re[aSq],Re[aSqLit[\[Tau]]]},{\[Tau],0,2},PlotStyle->{{Dashed,Red},{DotDashed,Blue}},
+PlotLegends->{Re[Subscript[A, FeynCalc]^2],Re[Subscript[A, Literature]^2]},AxesLabel->{\[Tau],Re[A^2]}]
+plot2=Plot[{Im[aSq],Im[aSqLit[\[Tau]]]},{\[Tau],0,2},PlotStyle->{{Dashed,Cyan},{DotDashed,
+Black}},PlotLegends->{Im[Subscript[A, FeynCalc]^2],Im[Subscript[A, Literature]^2]},AxesLabel->{\[Tau],Im[A^2]}]
+
+
+(* ::Subsection:: *)
+(*Check with the literature*)
+
+
+aSqres=(-4*\[Tau] + (-1 + \[Tau])*Log[1 + 2*(-1 + Sqrt[(-1 + \[Tau])/\[Tau]])*\[Tau]]^2)^2/(4*\[Tau]^4);
+Print["Check with the literature: ",
+			If[Simplify[(aSqres-aSq)]===0, 
+			"CORRECT.", "!!! WRONG !!!"]];

Examples/QCD/NRQCDVertexMatchingInTheTwoQuarkSectorOneLoop.m

@@ -0,0 +1,221 @@
+(* ::Package:: *)
+
+(* :Title: NRQCDVertexMatchingInTheTwoQuarkSectorOneLoop                                        *)
+
+(*
+	This software is covered by the GNU General Public License 3.
+	Copyright (C) 1990-2016 Rolf Mertig
+	Copyright (C) 1997-2016 Frederik Orellana
+	Copyright (C) 2014-2016 Vladyslav Shtabovenko
+*)
+
+(* :Summary:  Computation of the QCD side of the QCD/NRQCD matching in the
+			   two quark sector by expanding quark-gluon vertex in the relative
+			   momentum                                                     *)
+
+(* ------------------------------------------------------------------------ *)
+
+
+(* ::Section:: *)
+(*QCD side of the matching to NRQCD (quark-gluon vertex only) at  1-loop *)
+
+
+(* ::Text:: *)
+(*This example uses a custom QCD in background field formalism model created with FeynRules. Please evaluate the file FeynCalc/Examples/FeynRules/QCDBGF/GenerateModelQCDBGF.m before running it for the first time*)
+
+
+(* ::Subsection:: *)
+(*Load FeynCalc, FeynArts and FeynHelpers*)
+
+
+If[ $FrontEnd === Null,
+		$FeynCalcStartupMessages = False;
+		Print["Computation of the QCD side of the QCD/NRQCD matching in the 
+two quark sector by expanding quark-gluon vertex in the relative momentum"];
+];
+$LoadAddOns={"FeynHelpers"};
+$LoadFeynArts = True;
+<<FeynCalc`
+$FAVerbose=0;
+
+
+(* ::Subsection:: *)
+(*Generate Feynman diagrams*)
+
+
+topVertex = CreateTopologies[1,2 -> 1,ExcludeTopologies -> {WFCorrections}];
+diagsVertex = InsertFields[topVertex, {F[3, {1}],V[50]} ->
+		{F[3, {1}]}, InsertionLevel -> {Classes},		
+ Model -> FileNameJoin[{"QCDBGF","QCDBGF"}],GenericModel->FileNameJoin[{"QCDBGF","QCDBGF"}]];
+Paint[diagsVertex, ColumnsXRows -> {2, 1},ImageSize->{512,256}, Numbering -> None,SheetHeader->False];
+
+
+(* ::Subsection:: *)
+(*On-shell kinematics*)
+
+
+FCClearScalarProducts[];
+ScalarProduct[p1,p1]=m^2;
+ScalarProduct[p2,p2]=m^2;
+ScalarProduct[q,q]=q2;
+ScalarProduct[p1,p2]=-1/2FCI@ SPD[q,q]+m^2;
+
+
+(* ::Subsection:: *)
+(*Evaluation*)
+
+
+(* ::Text:: *)
+(*Since we are going to evaluate both diagrams one after another, it is convenient to define helper functions that will do the job.*)
+
+
+(* ::Text:: *)
+(*This function obtains the amplitude for further evaluation with FeynCalc*)
+
+
+ampFunc0[diag_]:=FCFAConvert[FCPrepareFAAmp[CreateFeynAmp[diag, 
+Truncated -> False,PreFactor->-1]],IncomingMomenta->{p1,q},
+OutgoingMomenta->{p2},LoopMomenta->{l},UndoChiralSplittings->True,
+DropSumOver->True,List->False,FinalSubstitutions->{SMP["m_u"]->m,
+q->p2-p1,Pair[Momentum[Polarization[___],___],___]:>1},
+ChangeDimension->D,SMP->True]//Contract//SUNSimplify//
+ReplaceAll[#,SUNTF[{x__},__]:>SUNT[x]]&//SUNSimplify[#,Explicit->True]&//
+ReplaceAll[#,SMP["g_s"]^3->4Pi SMP["alpha_s"]SMP["g_s"]]&
+
+
+(* ::Text:: *)
+(*Tensor reduction and simplification of the Dirac algebra*)
+
+
+ampFunc1[expr_]:=expr//TID[#,l,UsePaVeBasis->True,ToPaVe->True]&//DiracSimplify//
+Collect2[#,Spinor,Factoring->Simplify]&;
+
+
+(* ::Text:: *)
+(*Gordon decomposition*)
+
+
+ampFunc2[expr_]:=(expr/.{FCI[(FVD[p1,i_]+
+FVD[p2,i_])SpinorUBarD[p2,m_].SpinorUBarD[p1,m_]]:>
+FCI[SpinorUBarD[p2,m].(2m GAD[i]-I DiracSigma[GAD[i],
+GSD[p2-p1]]).SpinorUD[p1,m]]})//DotSimplify//Collect2[#,LorentzIndex]&;
+
+
+(* ::Text:: *)
+(*Expansion in q^2 up to first order*)
+
+
+ampFunc3[expr_]:=PaXEvaluateUVIRSplit[expr,PaXSeries->{{q2,0,1}},PaXAnalytic->True,
+PaXImplicitPrefactor->(2Pi)^(-D)]//FCHideEpsilon//Collect2[#,LorentzIndex]&;
+
+
+(* ::Text:: *)
+(*Extracts F_1(q^2)*)
+
+
+vectorPart[expr_,fo_:1]:=((expr//SelectFree2[#,DiracSigma]&)/
+FCI[I SMP["g_s"]SUNT[Glu2]SpinorUBarD[p2,m].GAD[Lor1].SpinorUD[p1,m]])//
+Collect2[#,{SMP["Delta_UV"],SMP["Delta_IR"],q2},Factoring->FullSimplify,FCFactorOut->fo]&;
+
+
+(* ::Text:: *)
+(*Extracts F_2(q^2)*)
+
+
+scalarPart[expr_,fo_:1]:=((expr//SelectNotFree2[#,DiracSigma]&)/
+FCI[I SMP["g_s"]SUNT[Glu2](I/(2m))FCI[SpinorUBarD[p2,m]. DiracSigma[GAD[Lor1],
+GSD[p2-p1]].SpinorUD[p1,m]]])//Collect2[#,{SMP["Delta_UV"],SMP["Delta_IR"],q2},
+Factoring->FullSimplify,FCFactorOut->fo]&;
+
+
+(* ::Subsection:: *)
+(*QCD abelian vertex*)
+
+
+abelianVertex=ampFunc0[DiagramExtract[diagsVertex,1]]
+
+
+abelianVertex1=ampFunc1[abelianVertex]
+
+
+abelianVertex2=ampFunc2[abelianVertex1]
+
+
+abelianVertex3=ampFunc3[abelianVertex2]
+
+
+(* ::Text:: *)
+(*This is  F_1^(V)*)
+
+
+resF1V=vectorPart[abelianVertex3,SMP["alpha_s"](CF-CA/2)1/Pi]
+
+
+(* ::Text:: *)
+(*This is  F_2^(V)*)
+
+
+resF2V=scalarPart[abelianVertex3,SMP["alpha_s"](CF-CA/2)1/Pi]
+
+
+(* ::Subsection:: *)
+(*QCD non-abelian vertex*)
+
+
+nonAbelianVertex=ampFunc0[DiagramExtract[diagsVertex,2]];
+
+
+nonAbelianVertex1=ampFunc1[nonAbelianVertex];
+
+
+nonAbelianVertex2=ampFunc2[nonAbelianVertex1]
+
+
+nonAbelianVertex3=ampFunc3[nonAbelianVertex2];
+
+
+(* ::Text:: *)
+(*This is  F_1^(g)*)
+
+
+resF1G=vectorPart[nonAbelianVertex3,SMP["alpha_s"]CA/(8Pi)]
+
+
+(* ::Text:: *)
+(*This is  F_2^(g)*)
+
+
+resF2G=scalarPart[nonAbelianVertex3,SMP["alpha_s"]CA/(8Pi)]
+
+
+(* ::Subsection:: *)
+(*Check with the literature*)
+
+
+(* ::Text:: *)
+(*Here we collect the results for F_1^(V) and F_2^(V) (abelian on-shell quark-gluon vertex) as well as  F_1^(g) and F_2^(g) (non-abelian on-shell quark-gluon vertex)  from Manohar (arXiv:hep-ph/9701294). The explicit values for these form-factors are given in Eqs. 29-30 and Eqs. 33-34 of the above-mentioned preprint.*)
+
+
+F1V=((SMP["alpha_s"]/Pi)(CF-1/2 CA)(1/(2EpsilonUV)+1/EpsilonIR+1-3/2Log[m/ScaleMu]+
+q2/m^2(-1/(3EpsilonIR)-1/8+1/3Log[m/ScaleMu])))
+
+
+F2V=(SMP["alpha_s"]/Pi)(CF-1/2 CA)(1/2+q2/(12m^2))
+
+
+F1G=(SMP["alpha_s"]/(8Pi))CA(2/EpsilonUV+4/EpsilonIR+4-6Log[m/ScaleMu]+
+q2/m^2(-3/EpsilonIR-1+3Log[m/ScaleMu]))
+
+
+F2G=(SMP["alpha_s"]/(8Pi))CA(4/EpsilonIR+6-4Log[m/ScaleMu]+
+q2/m^2(4/EpsilonIR+1-4Log[m/ScaleMu]))
+
+
+diff=Map[Simplify[(PowerExpand[FCShowEpsilon[#[[1]]]/.
+{ScaleMu^2->ScaleMu^2 E^EulerGamma/(4Pi),1/EpsilonIR->2/EpsilonIR,
+1/EpsilonUV->2/EpsilonUV}])-#[[2]],Assumptions->{ScaleMu>0,m>0}]&,{{resF1V,F1V},
+{resF2V,F2V},{resF1G,F1G},{resF2G,F2G}}]
+
+
+Print["Check with Manohar, arXiv:hep-ph/9701294: ",
+			If[diff==={0,0,0,0}, "CORRECT.", "!!! WRONG !!!"]];