Fixed a bug in DiracTrick that leads to incorrect

results when sufficiently nested dots are present in the input
(issue #176)

FeynCalc/Dirac/DiracSimplify.m

@@ -166,8 +166,8 @@
 
 			(* 	First of all we need to extract all Dirac structures in the input. *)
 			ex = FCDiracIsolate[ex,FCI->True,Head->dsHead, DotSimplify->True, DiracGammaCombine->OptionValue[DiracGammaCombine],
-				DiracSigmaExplicit->OptionValue[DiracSigmaExplicit], LorentzIndex->True, CartesianIndex->True,
-				ToDiracGamma67->optToDiracGamma67, DiracChain->OptionValue[DiracChain], Split->dsHeadAll];
+				DiracSigmaExplicit->OptionValue[DiracSigmaExplicit], LorentzIndex->All, CartesianIndex->True,
+				ToDiracGamma67->optToDiracGamma67, DiracChain->OptionValue[DiracChain], Split->dsHeadAll, "ExpandNestedDOTs"->True];
 
 			If[	!FreeQ[ex,DiracTrace] && !OptionValue[DiracTrace],
 				ex = ex /. dsHead[zz_]/; !FreeQ[zz,DiracTrace] :> zz

FeynCalc/Dirac/DiracTrace.m

@@ -259,7 +259,10 @@
 
 		time=AbsoluteTime[];
 		FCPrint[1,"DiracTrace: diracTraceEvaluate: Applying DiracTrick.", FCDoControl->diTrVerbose];
-		tmp = DiracTrick[tmp, FCI -> True, InsideDiracTrace->True, FCDiracIsolate->False];
+
+		(* 	If there are no nested DOTs, then FCDiracIsolate is not needed, otherwise one should use it *)
+			tmp = DiracTrick[tmp, FCI -> True, InsideDiracTrace->True, FCDiracIsolate->FeynCalc`Package`containsNestedDOTsQ[tmp]];
+
 		FCPrint[1,"DiracTrace: diracTraceEvaluate: DiracTrick done, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->diTrVerbose];
 		FCPrint[3,"DiracTrace: diracTraceEvaluate: After DiracTrick: ", tmp, FCDoControl->diTrVerbose];
 

FeynCalc/Dirac/DiracTrick.m

@@ -163,9 +163,10 @@ create replacement rules (standard)	and substitute the results back.
 			FCPrint[1, "DiracTrick: Extracting Dirac objects.", FCDoControl->diTrVerbose];
 			(* 	First of all we need to extract all the Dirac structures in the input. *)
 			If[	TrueQ[Head[optFCDiracIsolate]===List],
-				ex = FCDiracIsolate[ex,FCI->True,Head->dsHead, optFCDiracIsolate],
+				ex = FCDiracIsolate[ex,FCI->True,Head->dsHead, optFCDiracIsolate, "ExpandNestedDOTs"->True],
 				ex = FCDiracIsolate[ex,FCI->True,Head->dsHead, DotSimplify->True, DiracGammaCombine->OptionValue[DiracGammaCombine],
-					FCJoinDOTs -> optJoin, LorentzIndex->True, ToDiracGamma67-> OptionValue[ToDiracGamma67], DiracChain->OptionValue[DiracChain]];
+					FCJoinDOTs -> optJoin, LorentzIndex->True, ToDiracGamma67-> OptionValue[ToDiracGamma67], DiracChain->OptionValue[DiracChain],
+					"ExpandNestedDOTs"->True];
 			];
 
 			{freePart,dsPart} = FCSplit[ex,{dsHead}];
@@ -534,6 +535,13 @@ create replacement rules (standard)	and substitute the results back.
 
 		res = res/. DOT -> holdDOT;
 
+		(*	There should be no nested dots *)
+		If[	FeynCalc`Package`containsNestedDOTsQ[res,holdDOT],
+			Message[DiracTrick::failmsg,"Unexpanded nested DOTs detected!"];
+			Abort[]
+		];
+
+
 		If[	gamma5Present,
 			res = res /. holdDOT -> commonGamma5Properties /. commonGamma5Properties -> holdDOT;
 			gamma5Present = !FreeQ2[res,{DiracGamma[5],DiracGamma[6],DiracGamma[7]}];
@@ -689,7 +697,7 @@ create replacement rules (standard)	and substitute the results back.
 		If[	insideDiracTrace && res=!=0,
 			time=AbsoluteTime[];
 			FCPrint[2, "DiracTrick: diracTrickEval: Applying diracTraceSimplify again ", FCDoControl->diTrVerbose];
-			res = FCDiracIsolate[res, DotSimplify->False, FCI->True, DiracGammaCombine->False, FCJoinDOTs->optJoin, Head->dsHead];
+			res = FCDiracIsolate[res, DotSimplify->False, FCI->True, DiracGammaCombine->False, FCJoinDOTs->optJoin, Head->dsHead,"ExpandNestedDOTs"->True];
 			res = res  /. {dsHead[DiracGamma[5]] :> 0, dsHead[DiracGamma[6|7]] :> 1/2 } /. DOT-> holdDOT/.
 			dsHead[holdDOT[x__]] :>  diracTraceSimplify[x];
 			res = res /. diracTraceSimplify -> DOT /. dsHead-> Identity /. holdDOT -> DOT;

FeynCalc/Dirac/FCDiracIsolate.m

@@ -25,6 +25,9 @@
 "FCDiracIsolate failed to isolate Dirac structures in `1`!";
 
 Begin["`Package`"]
+
+containsNestedDOTsQ;
+
 End[]
 
 Begin["`FCDiracIsolate`Private`"]
@@ -55,6 +58,7 @@
 	IsolateNames		-> KK,
 	LorentzIndex		-> False,
 	Polarization		-> False,
+	"ExpandNestedDOTs"	-> False,
 	Spinor				-> True,
 	Split				-> True,
 	TimeConstrained		-> 3,
@@ -67,6 +71,12 @@
 			Cases[SelectNotFree[expr, heads],  l: (_LorentzIndex| _CartesianIndex) :> l[[1]] ,Infinity]]]
 ];
 
+containsNestedDOTsQ[expr_, dot_: DOT]:=
+Block[{chk,holdDOT2,holdDOT3,null1,null2},
+	chk = Cases[(expr /. dot->holdDOT2) + null1 + null2, holdDOT2[__],Infinity];
+	!FreeQ[holdDOT3@@@chk,holdDOT2]
+];
+
 holdDOT[]=1;
 
 FCDiracIsolate[a_ == b_, opts:OptionsPattern[]] :=
@@ -79,7 +89,7 @@
 	Block[{	res, null1, null2, ex,tmp, head, selectionList,
 			time, fcdiVerbose, headsList, headsOrig, optTimeConstrained,
 			optHead, headR, allHeads, allHeadsEval, headNoMatrix, collectList,
-			optSplit},
+			optSplit, aux},
 
 		If [OptionValue[FCVerbose]===False,
 			fcdiVerbose=$VeryVerbose,
@@ -172,9 +182,21 @@
 			time=AbsoluteTime[];
 			FCPrint[1, "FCDiracIsolate: Applying DotSimplify.", FCDoControl->fcdiVerbose];
 			tmp = FCSplit[ex, headsList, Expanding->OptionValue[Expanding]];
-			ex = tmp[[1]]+ DotSimplify[tmp[[2]],Expanding->False,FCI->True, FCJoinDOTs->OptionValue[FCJoinDOTs]];
+			aux = DotSimplify[tmp[[2]],Expanding->False,FCI->True, FCJoinDOTs->OptionValue[FCJoinDOTs]];
 			FCPrint[1, "FCDiracIsolate: Done applying DotSimplify, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcdiVerbose];
-			FCPrint[3, "FCDiracIsolate: After DotSimplify: ", ex, FCDoControl->fcdiVerbose]
+			FCPrint[3, "FCDiracIsolate: After DotSimplify: ", aux, FCDoControl->fcdiVerbose];
+
+			If[	OptionValue["ExpandNestedDOTs"],
+				If[	containsNestedDOTsQ[aux,DOT],
+					FCPrint[1, "FCDiracIsolate: Nested DOTs detected. Rerunning DotSimplify with Expanding set to True.", FCDoControl->fcdiVerbose];
+					time=AbsoluteTime[];
+					aux = DotSimplify[aux,Expanding->True,FCI->True, FCJoinDOTs->OptionValue[FCJoinDOTs]];
+					FCPrint[1, "FCDiracIsolate: Done applying DotSimplify, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcdiVerbose];
+					FCPrint[3, "FCDiracIsolate: After DotSimplify: ", aux, FCDoControl->fcdiVerbose];
+				];
+			];
+
+			ex = tmp[[1]]+ aux;
 		];
 
 		If[	OptionValue[FCTraceExpand],
@@ -203,6 +225,7 @@
 				ex + null1 + null2] /. {null1 | null2 -> 0} /. head[1] -> 1)
 		];
 		FCPrint[1, "FCDiracIsolate: Done isolating heads, timing: ", N[AbsoluteTime[] - time, 4], FCDoControl->fcdiVerbose];
+		FCPrint[3, "FCDiracIsolate: After isolating heads: ", ex, FCDoControl->fcdiVerbose];
 
 		If[ Together[(res /. restHead|head -> Identity)-ex] =!= 0,
 			Message[FCDiracIsolate::fail, ex];

Tests/Dirac/DiracTrace.test

@@ -256,5 +256,15 @@ FCDiracIsolate -> False]", "4 Pair[LorentzIndex[mu], LorentzIndex[nu]]"},
 {"fcstDiracTrace-ID53", "DiracTrace[DiracSigma[GA[mu], GA[nu]].GA[rho, si],
 DiracTraceEvaluate -> True, FCDiracIsolate -> False]", "4 I (Pair[LorentzIndex[mu], LorentzIndex[si]] Pair[LorentzIndex[nu],
 	LorentzIndex[rho]] - Pair[LorentzIndex[mu], LorentzIndex[rho]] Pair[LorentzIndex[nu],
-	LorentzIndex[si]])"}
+	LorentzIndex[si]])"},
+{"fcstDiracTrace-ID60", "FCSetDiracGammaScheme[\"NDR\"];\n DiracTrace[GS[p1 + p2 + p3] . GA[i1, i2, i3] . GA[5],
+DiracTraceEvaluate -> True, FCE -> True]", "(4*I)*(LC[i1, i2, i3][p1] + LC[i1, i2, i3][p2] + LC[i1, i2, i3][p3])"},
+{"fcstDiracTrace-ID61", "FCSetDiracGammaScheme[\"BMHV\"];\n DiracTrace[(GA[5] . GAD[\[Mu]] . GAD[lor1] . GAD[lorb] . GAD[lor2] .
+	GAD[lorz] . GAD[lor3] . GAD[\[Sigma]] -
+	GAD[\[Sigma]] . GAD[lor1] . GAD[lorb] . GAD[lor2] . GAD[lorz] .
+	GAD[lor3] . GA[5] . GAD[\[Mu]]) . (GSD[x] . GAD[\[Nu]] .
+	GAD[lor5] . GAD[lorz] . GAD[lor4]*FVD[x, lor1]*FVD[x, lor3]*
+	MTD[lor2, lor5]*MTD[lor4, lora]*MTD[lora, \[Rho]]*
+	MTD[lorb, \[Nu]]), DiracTraceEvaluate -> True, FCE -> True]",
+	"-16 I (56 - 30 D + 3 D^2) SPD[x, x] LC[\[Mu], \[Rho], \[Sigma]][x]"}
 });

Tests/Dirac/DiracTrick.test

@@ -1761,8 +1761,9 @@ InsideDiracTrace\[Rule]True]//FCE",
 {"fcstDiracTrickInsideDiracTrace-ID29",
 "DiracTrick[(m+GS[p]).(-fa GA[5].GA[nu]+f1 GA[nu]).(m+GS[pp]).(fa \
 GA[mu].GA[5]+f1 GA[mu]),InsideDiracTrace\[Rule]True]//FCE",
-"(m + GS[p]) . (fa*GA[nu] . GA[5] + f1*GA[nu]) . (m + GS[pp]) . \
-(fa*GA[mu] . GA[5] + f1*GA[mu])"},
+"f1^2*GS[p] . GA[nu] . GS[pp] . GA[mu] + fa^2*GS[p] . GA[nu] . GS[pp] . GA[mu] -
+f1*fa*GA[nu] . GS[pp] . GA[mu] . GS[p] . GA[5] + f1*fa*GS[p] .
+GA[nu] . GS[pp] . GA[mu] . GA[5] + f1^2*m^2*MT[mu, nu] - fa^2*m^2*MT[mu, nu]"},
 {"fcstDiracTrickInsideDiracTrace-ID30",
 "DiracTrick[GS[p].GA[nu, 5, mu, 5], InsideDiracTrace -> True]",
 "0"},

Tests/Dirac/FCDiracIsolate.test

@@ -311,6 +311,22 @@ c dAll[dch[
 	q2].Spinor[-Momentum[p4], m, 1]]]"},
 {"fcstFCDiracIsolate-ID35", "FCDiracIsolate[GA[i] == x, Head -> dch, FCE -> True]", "dch[GA[i]] == x"},
 {"fcstFCDiracIsolate-ID36", "FCDiracIsolate[{GA[i], GS[p]}, Head -> dch, FCE -> True]",
-"{dch[GA[i]], dch[GS[p]]}"}
+"{dch[GA[i]], dch[GS[p]]}"},
+{"fcstFCDiracIsolate-ID37", "FCDiracIsolate[(GA[5] . GAD[\[Mu]] . GAD[lor1] . GAD[lorb] .
+	GAD[lor2] . GAD[lorz] . GAD[lor3] . GAD[\[Sigma]] -
+	GAD[\[Sigma]] . GAD[lor1] . GAD[lorb] . GAD[lor2] . GAD[lorz] .
+	GAD[lor3] . GA[5] . GAD[\[Mu]]) . (GSD[x] . GAD[\[Nu]] .
+	GAD[lor5] . GAD[lorz] . GAD[lor4]*FVD[x, lor1]*FVD[x, lor3]*
+	MTD[lor2, lor5]*MTD[lor4, lora]*MTD[lora, \[Rho]]*
+	MTD[lorb, \[Nu]]), DotSimplify -> True, Head -> dsHead,
+FCE -> True,\"ExpandNestedDOTs\"->True]","dsHead[GA[5] . GAD[\[Mu]] . GAD[lor1] . GAD[lorb] . GAD[lor2] .
+	GAD[lorz] . GAD[lor3] . GAD[\[Sigma]] . GSD[x] . GAD[\[Nu]] .
+	GAD[lor5] . GAD[lorz] . GAD[lor4]] FVD[x, lor1] FVD[x, lor3] MTD[
+lor2, lor5] MTD[lor4, lora] MTD[lora, \[Rho]] MTD[lorb, \[Nu]] -
+dsHead[GAD[\[Sigma]] . GAD[lor1] . GAD[lorb] . GAD[lor2] .
+	GAD[lorz] . GAD[lor3] . GA[5] . GAD[\[Mu]] . GSD[x] .
+	GAD[\[Nu]] . GAD[lor5] . GAD[lorz] . GAD[lor4]] FVD[x, lor1] FVD[
+x, lor3] MTD[lor2, lor5] MTD[lor4, lora] MTD[lora, \[Rho]] MTD[
+lorb, \[Nu]]"}
 });
 

Tests/Dirac/Tr.test

@@ -314,7 +314,9 @@ ExpandScalarProduct // Expand // FCE; \n FCClearScalarProducts[]; \n tmp",
 	ExpandScalarProduct // Expand // FCE; \n FCClearScalarProducts[]; \n tmp",
 	"-((4 SP[l, p]^2)/SP[l, o]) - (4 eps SP[l, p]^2)/SP[l, o] - (4 SP[l, p] SP[o, p])/SP[l, o] - (4 eps SP[l, p] SP[o, p])/SP[l, o]"},
 {"fcstTr-ID32", "Tr[FV[k1 + k2, mu] (GA[al] FV[k1, al] + m).(GA[mu] f1 - 1/2 I/m 1/2 I (GA[mu, nu] - GA[nu, mu]) FV[k1 - k2, nu] f2).(GS[k1] + m)]//FCE",
-"8 f1 m (SP[k1, k1] + SP[k1, k2])"},
+"(4*(2*f1*m^2*SP[k1, k1] - f2*SP[k1, k1]*SP[k1, k1 - k2] + \
+2*f1*m^2*SP[k1, k2] - f2*SP[k1, k1 - k2]*SP[k1, k2] + f2*SP[k1, \
+k1]*SP[k1, k1 + k2] - f2*SP[k1, k2]*SP[k1, k1 + k2]))/m"},
 {"fcstTr-ID33", "Tr[1/32 GA[m].(1 - GA[5]).GS[k].(1 + GA[5]).GA[n].(1 - GA[5]).GS[p].(1 + GA[5])] // FCE",
 "FV[k, n] FV[p, m] + FV[k, m] FV[p, n] - MT[m, n] SP[k, p] +I LC[m, n][k, p]"},
 {"fcstTr-ID34", "Tr[GA[5].(GA[5] + GS[p]).GA[5].GS[q]]//FCE", "-4 SP[p, q]"},