(* ::Package:: *)

Quit


(* ::Subsection::Closed:: *)
(*Setup*)


(* Find a set with the smallest number of parameters for which u or f vanish when they are set to 0 *)
Options[FindSetS] = {PrintMessage -> True};
FindSetS[u_,f_,ti_List,opts:OptionsPattern[]] := Module[{i,testsets},
i=1;
While[i<=Length[ti],
  testsets = Subsets[ti,{i}];
  Table[
    If[TrueQ[
      FullSimplify[(u /. Table[testsets[[j,k]]->0,{k,1,i}]) == 0] || 
      FullSimplify[(f /. Table[testsets[[j,k]]->0,{k,1,i}]) == 0]
      ],
      Return[testsets[[j]],Module];
    ],
  {j,1,Length[testsets]}];
  i++;
];
If[TrueQ[OptionValue[PrintMessage]],Print["Singularities have been decomposed."]]; Return[Null]
]


Options[SectorDecompose] = {MaxSteps -> 200};

SectorDecompose[fac_,{u_,uexp_},{f_,fexp_},nprops_,nloops_,opts:OptionsPattern[]] := Module[{steps,decomposed,sectors,sec,setS,subst,unew,fnew,facnew},
steps = 0;
(* list of decomposed sectors *)
decomposed = {};
(* create primary sectors and add them to list of sectors to be decomposed *)
sectors = Table[
  {i,Simplify[{fac,u/ToExpression["x"<>ToString[i]]^nloops,f/ToExpression["x"<>ToString[i]]^(nloops+1)} /. Table[
      ToExpression["x"<>ToString[j]] -> ToExpression["x"<>ToString[i]]*ToExpression["t"<>ToString[j]],
      {j,Delete[Table[k,{k,1,nprops}],{i}]}]
    ]},
  {i,1,nprops}];
(* while sectors is not empty, decompose the first element *)
While[Length[sectors]>0 && steps<OptionValue[MaxSteps],
sec = sectors[[1]];
sectors = Delete[sectors,{1}];
(* determine the set S for the sector under consideration *)
(* We only decompose U or F when they appear in the denominator, otherwise they cause no singularities *)
setS = FindSetS[If[TrueQ[uexp<0/.eps->0],sec[[2,2]],1],If[TrueQ[fexp<0/.eps->0],sec[[2,3]],1],
  Delete[Table[ToExpression["t"<>ToString[i]],{i,1,nprops}],{ToExpression[StringTake[ToString[sec[[1]]],1]]}],PrintMessage -> False];
(* If no set S exists add to list of decomposed sectors, otherwise decompose *)
If[Head[setS]=!=List,AppendTo[decomposed,sec],
  (* decomposition of sec into Length[setS] sectors which are added to list of remaining sectors *)
  Table[
     (* substitution of variables *)
     subst = Table[setS[[j]]->setS[[i]]*If[j==i,1,setS[[j]]],{j,1,Length[setS]}];
     unew = Factor[sec[[2,2]] /. subst];
     fnew = Factor[sec[[2,3]] /. subst];
     facnew = (sec[[2,1]]/.subst)*setS[[i]]^(Length[setS]-1)(*Jacobian*);
     (* If unew or fnew factorizes pull out the factor of ti *)
     If[TrueQ[FullSimplify[unew == 0 /. {setS[[i]]->0}]],
     facnew*=setS[[i]]^uexp; unew = FullSimplify[unew/setS[[i]]];];
     If[TrueQ[FullSimplify[fnew == 0 /. {setS[[i]]->0}]],
     facnew*=setS[[i]]^fexp; fnew = FullSimplify[fnew/setS[[i]]];];
     (* add ith new sector *)
     AppendTo[sectors,
     {ToExpression[ToString[sec[[1]]]<>StringDrop[ToString[setS[[i]]],{1}]],
     {Simplify[facnew,Assumptions->And@@(Table[ToExpression["t"<>ToString[j]]>0,{j,1,nprops}])], unew, fnew}}],
  {i,1,Length[setS]}];
];
steps++;
];
If[Length[sectors]>0,Print["Decomposition incomplete. Remaining:\n",sectors]];
decomposed
]


(* Integration only works with single poles, integration goes up to finite order *)
Options[IntegrateSectors] = {PrintStatus -> False};

IntegrateSectors[list_List,{uexp_,fexp_},ti_List,reps_List,opts:OptionsPattern[]] := Module[{n,res,tiint,exps,poles,Int,NInt},
n=1; res = 0;
While[n<=Length[list],
  tiint = Complement[ti,{ToExpression["t"<>StringTake[ToString[list[[n,1]]],{1}]]}];
  exps = Exponent[list[[n,2,1]],tiint];
  Table[If[(exps[[i]]/.eps->0)<-1,Print["Higher poles not implemented!"];Abort[]],{i,1,Length[tiint]}];
  poles = {}; Table[If[(exps[[i]]/.eps->0)==-1,AppendTo[poles,tiint[[i]]]],{i,1,Length[tiint]}];
  (* Add integral and then process with subtractions *)
  res += Int[list[[n,2,1]]*list[[n,2,2]]^uexp*list[[n,2,3]]^fexp /. reps,tiint,Length[poles]];
  While[Length[poles]>0,
  res = res /. {Int[expr_,vars_,epspow_] :> 
     Int[expr - (Limit[expr*poles[[1]]^(-Exponent[list[[n,2,1]],poles[[1]]]),poles[[1]]->0,Assumptions->And@@(Table[ti[[i]]>0,{i,1,Length[ti]}])])*poles[[1]]^Exponent[list[[n,2,1]],poles[[1]]],vars,epspow-1] +
     1/(1+Exponent[list[[n,2,1]],poles[[1]]])*Int[Limit[expr*poles[[1]]^(-Exponent[list[[n,2,1]],poles[[1]]]),poles[[1]]->0,Assumptions->And@@(Table[ti[[i]]>0,{i,1,Length[ti]}])],Complement[vars,{poles[[1]]}],epspow]};
  poles = Delete[poles,{1}];
  ];
  res = res /. {Int[expr_,vars_,epspow_] :> NInt[CoefficientList[Series[expr,{eps,0,epspow}],eps],vars]};
  res = res /. {NInt[coeffs_,vars_] :> Sum[
           NInt[Simplify[coeffs[[i]]],Sequence@@(Table[{vars[[j]],0,1},{j,1,Length[vars]}])]*eps^(i-1),
        {i,1,Length[coeffs]}]};
  res = Expand[res /. {NInt[expr_] -> expr} /. NInt -> NIntegrate];
  
  If[TrueQ[OptionValue[PrintStatus]],Print[DateString[],"     Sector ",n," of ",Length[list]," complete."]];
  n++;
];
res
];


(* ::Subsection::Closed:: *)
(*Two-loop vertex: Examples for decomposed sectors*)


(* Consider the planar 2-loop vertex with 2 on-shell lines with indices equal to one at Euclidean momentum q^2 = -1 *)
FeynmanParametrization[{{l1+p1,0,1},{l1+l2+p1,0,1},{l1+l2-p2,0,1},{l1-p2,0,1},{l1,0,1},{l2,0,1}},{l1,l2},{p1,p2},Replacements->{p1^2->0,p2^2->0,p1*p2->q2/2},PrintMQJ->True,GraphPolynomials->True]


(* Fhat = Int[fac*U^(3*eps)*F^(-2-2*eps)*Delta[1-Sum[xi]],{xi,0,Infinity}] with *)


fac = 1;
U = (x2+x3) (x1+x4+x5)+(x1+x2+x3+x4+x5) x6;
F = -q2 (x1 x3 x4 + x1 x2 (x3+x4) + x2 x3 (x4+x5) + (x1+x2) (x3+x4) x6) /. {q2->-1};


(* Example 1 *)


(* Take the 6th primary sector as an example *)
(* Fhat6 = Int[fac6*U6^(3*eps)*F6^(-2-2*eps),{ti,0,1}] with *) 
fac6 = Simplify[fac /. Table[ToExpression["x"<>ToString[i]]->x6*ToExpression["t"<>ToString[i]],{i,1,5}]]
U6 = Simplify[U/x6^2 /. Table[ToExpression["x"<>ToString[i]]->x6*ToExpression["t"<>ToString[i]],{i,1,5}]]
F6 = Simplify[F/x6^3 /. Table[ToExpression["x"<>ToString[i]]->x6*ToExpression["t"<>ToString[i]],{i,1,5}]]


FindSetS[U6,F6,{t1,t2,t3,t4,t5}]


(* Since the set S has 2 elements we must split the {t1,t2} square into two sectors *)
(* Consider the sector 61 as an example *)
Factor[U6 /. {t2->t1*t2}]
Factor[F6 /. {t2->t1*t2}]


(* We see that t1 can be pulled out of F6. Thus (taking into account the Jacobian factor t1)*)
(* Fhat61 = Int[fac61*U61^(3*eps)*F61^(-2-2*eps),{ti,0,1}] with *) 
fac61 = Simplify[fac6*t1^(1-2-2*eps) /. {t2->t1*t2}]
U61 = Simplify[U6 /. {t2->t1*t2}]
F61 = Simplify[F6/t1 /. {t2->t1*t2}]


FindSetS[U61,F61,{t1,t2,t3,t4,t5}]


(* Since the set S has 2 elements we must split the {t3,t4} square into two sectors *)
(* Consider the sector 613 as an example *)
Factor[U61 /. {t4->t3*t4}]
Factor[F61 /. {t4->t3*t4}]


(* We see that t3 can be pulled out of F61. Thus (taking into account the Jacobian factor t3)*)
(* Fhat613 = Int[fac613*U613^(3*eps)*F613^(-2-2*eps),{ti,0,1}] with *) 
fac613 = Simplify[fac61*t3^(1-2-2*eps) /. {t4->t3*t4}]
U613 = Simplify[U61 /. {t4->t3*t4}]
F613 = Simplify[F61/t3 /. {t4->t3*t4}]


(* Due to the constant term in F613 all the IR singularities in this sector have been decomposed *)
FindSetS[U613,F613,{t1,t2,t3,t4,t5}]


(* Since the set S has 3 elements we must split the {t1,t3,t5} cube into three sectors *)
(* Consider the sector 6131 as an example *)
Factor[U613 /. {t3->t1*t3,t5->t1*t5}]
Factor[F613 /. {t3->t1*t3,t5->t1*t5}]


(* We see that t1 can be pulled out of U631. Thus (taking into account the Jacobian factor t1^2)*)
(* Fhat6131 = Int[fac6131*U6131^(3*eps)*F6131^(-2-2*eps),{ti,0,1}] with *) 
fac6131 = Simplify[fac613*t1^(2+3*eps) /. {t3->t1*t3,t5->t1*t5},Assumptions->And@@(Table[ToExpression["t"<>ToString[i]]>0,{i,1,5}])]
U6131 = Simplify[U613/t1 /. {t3->t1*t3,t5->t1*t5}]
F6131 = Simplify[F613 /. {t3->t1*t3,t5->t1*t5}]


(* At this point both graph polynomials contain constant terms implying that the decomposition is complete in this sector *)
FindSetS[U6131,F6131,{t1,t2,t3,t4,t5}]


(* Example 2 *)


(* Take the 1st primary sector as an example *)
(* Fhat1 = Int[fac1*U1^(3*eps)*F1^(-2-2*eps),{ti,0,1}] with *) 
fac1 = Simplify[fac /. Table[ToExpression["x"<>ToString[i]]->x1*ToExpression["t"<>ToString[i]],{i,2,6}],Assumptions->And@@(Table[ToExpression["t"<>ToString[i]]>0,{i,2,6}])]
U1 = Simplify[U/x1^2 /. Table[ToExpression["x"<>ToString[i]]->x1*ToExpression["t"<>ToString[i]],{i,2,6}]]
F1 = Simplify[F/x1^3 /. Table[ToExpression["x"<>ToString[i]]->x1*ToExpression["t"<>ToString[i]],{i,2,6}]]


FindSetS[U1,F1,{t2,t3,t4,t5,t6}]


(* Since the set S has 2 elements we must split the {t3,t4} square into two sectors *)
(* Consider the sector 13 as an example *)
Factor[U1 /. {t4->t3*t4}]
Factor[F1 /. {t4->t3*t4}]


(* We see that t3 can be pulled out of F1. Thus (taking into account the Jacobian factor t3)*)
(* Fhat13 = Int[fac13*U13^(3*eps)*F13^(-2-2*eps),{ti,0,1}] with *) 
fac13 = Simplify[fac1*t3^(1-2-2*eps) /. {t4->t3*t4},Assumptions->And@@(Table[ToExpression["t"<>ToString[i]]>0,{i,2,6}])]
U13 = Simplify[U1 /. {t4->t3*t4}]
F13 = Simplify[F1/t3 /. {t4->t3*t4}]


FindSetS[U13,F13,{t2,t3,t4,t5,t6}]


(* Since the set S has 3 elements we must split the {t2,t3,t6} cube into three sectors *)
(* Consider the sector 132 as an example *)
Factor[U13 /. {t3->t2*t3,t6->t2*t6}]
Factor[F13 /. {t3->t2*t3,t6->t2*t6}]


(* We see that t2 can be pulled out of U13 and F13. Thus (taking into account the Jacobian factor t1^2)*)
(* Fhat132 = Int[fac132*U132^(3*eps)*F6131^(-2-2*eps),{ti,0,1}] with *) 
fac132 = Simplify[fac13*t2^(2+3*eps-2-2*eps) /. {t3->t2*t3,t6->t2*t6},Assumptions->And@@(Table[ToExpression["t"<>ToString[i]]>0,{i,2,6}])]
U132 = Simplify[U13/t2 /. {t3->t2*t3,t6->t2*t6}]
F132 = Simplify[F13/t2 /. {t3->t2*t3,t6->t2*t6}]


(* At this point both graph polynomials contain constant terms implying that the decomposition is complete in this sector *)
FindSetS[U132,F132,{t2,t3,t4,t5}]


(* ::Subsection::Closed:: *)
(*Onel-loop triangle: Complete decomposition and numerical integration*)


FeynmanParametrization[{{l1,0,1},{l1+p1,0,1},{l1+p1+p2,0,1}},{l1},{p1,p2},Replacements->{p1^2->0,p2^2->0,p1*p2->s/2},GraphPolynomials->True]//Simplify


secs = SectorDecompose[1,{x1+x2+x3,-1+2*eps},{-s x1 x3,-1-eps},3,1]


secsint = IntegrateSectors[secs,{-1+2*eps,-1-eps},{t1,t2,t3},{s->-1}]
Series[-Gamma[1+eps]*Exp[eps*EulerGamma]*secsint,{eps,0,0}]


<<X`


LoopRefine[LoopIntegrate[1,l1,{l1,0,1},{l1+p1,0,1},{l1+p1+p2,0,1}]/.{p1.p1->0,p2.p2->0,p1.p2->s/2}]
(* At symmetric Euclidean point this gives *)
%/.{s->-1,Mu->1}
%//N


(* ::Subsection::Closed:: *)
(*Onel-loop box: Complete decomposition and numerical integration*)


FeynmanParametrization[{{l1,0,1},{l1+p1,0,1},{l1+p1+p2,0,1},{l1+p3,0,1}},{l1},{p1,p2,p3},Replacements->{p1^2->0,p2^2->0,p3^2->0,p1*p2->s/2,p1*p3->-t/2,p2*p3->(s+t)/2},GraphPolynomials->True]//Simplify


secs = SectorDecompose[1,{x1+x2+x3+x4,2*eps},{-s x1 x3 - t x2 x4,-2-eps},4,1]
Length[secs]


secsint = IntegrateSectors[secs,{2*eps,-2-eps},{t1,t2,t3,t4},{s->-1,t->-1}]
Series[Gamma[2+eps]*Exp[eps*EulerGamma]*secsint,{eps,0,0}]


<<X`


LoopRefine[LoopIntegrate[1,l1,{l1,0,1},{l1+p1,0,1},{l1+p1+p2,0,1},{l1+p3,0,1}]/.{p1.p1->0,p2.p2->0,p3.p3->0,p1.p2->s/2,p1.p3->-t/2,p2.p3->(s+t)/2}]
(* At symmetric Euclidean point this gives *)
%/.{s->-1,t->-1,Mu->1}
%//N


(* ::Subsection::Closed:: *)
(*Two-loop vertex: Complete decomposition and numerical integration*)


secs = SectorDecompose[1,{(x2+x3) (x1+x4+x5)+(x1+x2+x3+x4+x5) x6,3*eps},{(x1 x3 x4 + x1 x2 (x3+x4) + x2 x3 (x4+x5) + (x1+x2) (x3+x4) x6),-2-2eps},6,2]


Length[secs]


(* Sectors with eps^(-4) poles *)
{51432,{(t2 t3)^eps (t1 t2 t3 t4)^(-1-2 eps),(t1+t4) (1+t1 t3+t2 t4)+t6+t1 (1+t2) t3 t6+t2 (1+t3) t4 t6,1+t1+t1 t3+t4+t2 t4+(1+t2) (1+t3) t6}}
{51436,{(t3 t6)^eps (t1 t3 t4 t6)^(-1-2 eps),(1+t1 t3) (1+t1 t2+t4)+(t4+(t1 t2+t4) (t3+t4)) t6,1+t3+t4+t2 (1+t1+t1 t3+(1+t3+t4) t6)}}


(* Coefficient of eps^(-4) pole agrees with result 1/4 from literature, e.g. (6.22) in Smirnov's book *)
Integrate[(t2 t3)^eps (t1 t2 t3 t4)^(-1-2 eps),{t1,0,1},{t2,0,1},{t3,0,1},{t4,0,1},GenerateConditions->False]*Integrate[(1+t1+t1 t3+t4+t2 t4+(1+t2) (1+t3) t6)^(-2)/.{t1->0,t2->0,t3->0,t4->0},{t5,0,1},{t6,0,1}] +
Integrate[(t3 t6)^eps (t1 t3 t4 t6)^(-1-2 eps),{t1,0,1},{t3,0,1},{t4,0,1},{t6,0,1},GenerateConditions->False]*Integrate[(1+t3+t4+t2 (1+t1+t1 t3+(1+t3+t4) t6))^(-2)/.{t1->0,t3->0,t4->0,t6->0},{t2,0,1},{t5,0,1}]


(*Dies at some point*)
(*secsint = IntegrateSectors[secs,{3*eps,-2-2eps},{t1,t2,t3,t4,t5,t6},{},PrintStatus \[Rule] True]
Series[Gamma[2+2 eps]*Exp[2*eps*EulerGamma]*secsint,{eps,0,0}]*)