(* ::Package:: *)

(* ::Subsection:: *)
(*Formula for general Feynman parametrization*)


(*Implement general Feynman parametrization from 1502.06595*)

(* props = {{q1,m1,nu1}, ... , {qn,mn,nun}} *)
(* li = list of loop momenta *)
(* pi = list of external momenta *)

x::usage = "Symbol for Feynman parameters.";
Dimension::usage = "Dimensionality of the loop integrals.";
eps::usage = "Dimensional regulator, d = 4 - 2*eps.";
Replacements::usage = "Allows to set replacement rules for scalar products.";
PrintMQJ::usage = "Option to display the quantities M, Q and J.";
GraphPolynomials::usage = "Option to display the graph polynomials.";

Options[FeynmanParametrization] = {
Dimension -> 4 - 2*eps, 
Replacements -> {},
PrintMQJ -> False,
GraphPolynomials -> False};

FeynmanParametrization[props_List,li_List,pi_List,opts:OptionsPattern[]] := Module[
{Nloops,Next,Nprops,nu,Nnu,dim,reps,alpha,beta,M,Q,J,U,F,FeynmanIntegrand},

Nloops = Length[li];
Next = Length[pi];
Nprops = Length[props];
nu = Table[props[[j,3]],{j,1,Nprops}];
Nnu = Sum[nu[[j]],{j,1,Nprops}];

dim = OptionValue[Dimension];
reps = OptionValue[Replacements];

alpha = Table[Coefficient[props[[j,1]],li[[l]]],{j,1,Nprops},{l,1,Nloops}];
beta = Table[-Coefficient[props[[j,1]],pi[[e]]],{j,1,Nprops},{e,1,Next}];

M = Table[Sum[x[j]*alpha[[j,l1]]*alpha[[j,l2]],{j,1,Nprops}],{l1,1,Nloops},{l2,1,Nloops}];
Q = Table[Sum[x[j]*alpha[[j,l]]*beta[[j,e]]*pi[[e]],{j,1,Nprops},{e,1,Next}],{l,1,Nloops}];
J = Sum[x[j]*(Sum[beta[[j,e]]*pi[[e]],{e,1,Next}]^2-props[[j,2]]^2),{j,1,Nprops}]//.reps;

If[TrueQ[OptionValue[PrintMQJ]],
  Print["M = ",M];
  Print["Q = ",Q];
  Print["J = ",J];
];

U = Expand[Det[M]];
F = Simplify[Expand[U*(Dot[Q,Inverse[M],Q] - J)]//.reps];

If[TrueQ[OptionValue[GraphPolynomials]],
  Print["U = ",U];
  Print["F = ",F];
];

FeynmanIntegrand = (-1)^Nnu * Gamma[Nnu-Nloops*dim/2]/Product[Gamma[nu[[j]]],{j,1,Nprops}] * 
  Product[x[j]^(nu[[j]]-1),{j,1,Nprops}] * U^Collect[Nnu - (Nloops+1)*dim/2,eps] * F^Collect[Nloops*dim/2 - Nnu,eps];

FeynmanIntegrand
];


(* ::Subsubsection::Closed:: *)
(*Example 1*)


(* Fully massive 2-point function *)
FeynmanParametrization[{{l1+q/2,m,1},{l1-q/2,m,1}},{l1},{q},Replacements->{q^2->s},GraphPolynomials->True]
FullSimplify[%/.{x[1]->x1,x[2]->1-x1}]


(* ::Subsubsection::Closed:: *)
(*Example 2*)


(* 2-loop massless propagator integral *)
FeynmanParametrization[{{l1,0,1},{l1-q,0,1},{l2,0,1},{l2-q,0,1},{l1-l2,0,1}},{l1,l2},{q},Replacements->{q^2->s},GraphPolynomials->True]
% /. {eps->0,x[i_]:>ToExpression["x"<>ToString[i]]}
(* Choose DiracDelta[1-x5] *)
Integrate[%/.x5->1,{x1,0,Infinity},{x2,0,Infinity},{x3,0,Infinity},{x4,0,Infinity}]


(* ::Subsubsection::Closed:: *)
(*Exercise*)


(* 2-loop tadpole with one mass *)
FeynmanParametrization[{{l1,0,a1},{l2,0,a2},{l1+l2,m,a3}},{l1,l2},{},GraphPolynomials->True] /. {x[i_]:>ToExpression["x"<>ToString[i]]}
(* Choose DiracDelta[1-x3] to simplify the factor (x1+x2) in the F polynomial *)
Simplify[%/.{x3->1},Assumptions->m>0]
Normal[Integrate[%,{x2,0,Infinity},Assumptions->x1>0]]
Normal[Integrate[%,{x1,0,Infinity}]]


(* 2-loop tadpole with two identical masses *)
FeynmanParametrization[{{l1,m,a1},{l2,m,a2},{l1+l2,0,a3}},{l1,l2},{},GraphPolynomials->True] /. {x[i_]:>ToExpression["x"<>ToString[i]]}
(* Choose DiracDelta[1-x1-x2] to simplify the factor (x1+x2) in the F polynomial *)
Simplify[%/.{x2->1-x1},Assumptions->m>0]
Normal[Integrate[%,{x3,0,Infinity},Assumptions->0<x1<1]]
Normal[Integrate[%,{x1,0,1}]]


(* ::Subsection::Closed:: *)
(*Manual specification of procedure*)


(* Input is a scalar integral MyInt[const,{{prop1,i1},...,{propN,iN}},{l1,...,lL},{p1,...,pP},{}] *)

(* When no more loop integrals are left write out integrand explicitly *)
MyInt[num_,props_List,{},pi_List,xi_List] := MyInt[Simplify[
    num*Product[props[[i,1]]^(-props[[i,2]]),{i,1,Length[props]}],
  Assumptions->And@@Table[0<xi[[i]]<1 && 0<ToExpression[ToString[xi[[i]]]<>"bar"]<1,{i,1,Length[xi]}]],
  xi];

(* denote scalar producs by sp[a,b] *)
Attributes[sp]={Orderless};
sp[a_Plus,b_] := Sum[sp[a[[i]],b],{i,1,Length[a]}];
sp[a_,b_Plus] := Sum[sp[a,b[[i]]],{i,1,Length[b]}];
sp[Times[n_,a__],b_] := n*sp[Times[a],b] /; NumericQ[n];
sp[a_,Times[n_,b__]] := n*sp[a,Times[b]] /; NumericQ[n];

(* combine propagators n1 and n2 with usual Feynman parametrization: *)
CombinePropagators[MyInt[num_,props_List,li_List,pi_List,xi_List],{n1_,n2_}] := Module[{x,xbar,in1,in2,propn1n2},
x = ToExpression["x"<>ToString[Length[xi]+1]];
xbar = ToExpression[ToString[x]<>"bar"];
in1 = props[[n1,2]]; in2 = props[[n2,2]];
propn1n2 = {FullSimplify[x*props[[n1,1]] + (1-x)*props[[n2,1]]]/.{1-x->xbar,x-1->-xbar},in1+in2};

MyInt[num*Gamma[in1+in2]/(Gamma[in1]Gamma[in2])*x^(in1-1)xbar^(in2-1),
Delete[ReplacePart[props,n1->propn1n2],{n2}],li,pi,Flatten[Append[xi,{x}]]]
];

(* Perform integral over loop momentum li[[ln]] taking into account only the propagator n *)
(* This assumes that all instances of the considered loop momenta have been combined into one propagator *)
IntegrateLoop[MyInt[num_,props_List,li_List,pi_List,xi_List],{ln_,n_}] := Module[{prop,in,lint,c0,moms},
prop = props[[n,1]]; in = props[[n,2]]; lint = li[[ln]];
c0 = Coefficient[prop,sp[lint,lint]]; If[TrueQ[c0==0],Print["Propagator must be quadratic in loop momentum!"];Abort[]];
(* normalize propagator,rescale numerator in output: *) prop /= -c0; 
(* complete the square *)
moms = Union[li[[;;ln-1]],li[[ln+1;;]],pi];
prop = Expand[prop/.{lint -> lint + Sum[Coefficient[prop,sp[lint,moms[[i]]]]*moms[[i]]/2,{i,1,Length[moms]}]}];

prop = prop + sp[lint,lint] //. 
  Table[ sp[Times[a___,xi[[i]],b___],c_] -> xi[[i]]*sp[Times[a,b],c], {i,1,Length[xi]}] //. 
  Table[ sp[c_,Times[a___,xi[[i]],b___]] -> xi[[i]]*sp[c,Times[a,b]], {i,1,Length[xi]}] //. 
  Table[ sp[Times[a___,ToExpression[ToString[xi[[i]]]<>"bar"],b___],c_] -> ToExpression[ToString[xi[[i]]]<>"bar"]*sp[Times[a,b],c], {i,1,Length[xi]}] //. 
  Table[ sp[c_,Times[a___,ToExpression[ToString[xi[[i]]]<>"bar"],b___]] -> ToExpression[ToString[xi[[i]]]<>"bar"]*sp[c,Times[a,b]], {i,1,Length[xi]}];
  
prop = FullSimplify[prop] /. 
  Table[1-xi[[i]] -> ToExpression[ToString[xi[[i]]]<>"bar"],{i,1,Length[xi]}] /. 
  Table[xi[[i]]-1 -> -ToExpression[ToString[xi[[i]]]<>"bar"],{i,1,Length[xi]}];

MyInt[Gamma[in-2+eps]/Gamma[in]*num/(-c0)^in,ReplacePart[props,n->{prop,in-2+eps}],Delete[li,{ln}],pi,xi]
];

(* Normalizes the propagator n s.t. the coefficient of sp[li[[ln]],li[[ln]]] is -1 *)
NormalizePropagator[MyInt[num_,props_List,li_List,pi_List,xi_List],{ln_,n_}] := Module[{prop,in,lint,c0},
prop = props[[n,1]]; in = props[[n,2]]; lint = li[[ln]];
c0 = -Coefficient[prop,sp[lint,lint]]; If[TrueQ[c0==0],Print["Propagator must be quadratic in loop momentum!"];Abort[]];

MyInt[If[TrueQ[Head[c0]==Times],num*Product[c0[[i]]^(-in),{i,1,Length[c0]}],num*c0^(-in)],ReplacePart[props,n->{Expand[prop/(c0)],in}],li,pi,xi]
]; 


int10111 = MyInt[1,{{-sp[l1,l1]-2sp[l1,q],1},{-sp[l1,l1]+mH2,1},{-sp[l2,l2]+mH2,1},{-sp[l1-l2,l1-l2]+mH2,1}},{l1,l2},{q},{}];
int10111 = CombinePropagators[int10111,{4,3}];
int10111 = IntegrateLoop[int10111,{2,3}];
int10111 = NormalizePropagator[int10111,{1,3}];
int10111 = CombinePropagators[int10111,{3,2}];
int10111 = CombinePropagators[int10111,{2,1}];
int10111 = IntegrateLoop[int10111,{1,1}]


int10111 = MyInt[1,{{-sp[l1,l1]-2sp[l1,q],1},{-sp[l2,l2]-2sp[l2,q],0},{-sp[l1,l1]+mH2,1},{-sp[l2,l2]+mH2,1},{-sp[l1-l2,l1-l2]+mH2,1}},{l1,l2},{q},{}];
int10111 = CombinePropagators[int10111,{5,4}];
int10111 = IntegrateLoop[int10111,{2,4}];
int10111 = NormalizePropagator[int10111,{1,4}];
int10111 = CombinePropagators[int10111,{4,3}];
int10111 = CombinePropagators[int10111,{3,1}];
int10111 = IntegrateLoop[int10111,{1,2}]


intni = MyInt[1,{{-sp[l1,l1]-2sp[l1,q],n1},{-sp[l1,l1]+mH2,n3},{-sp[l2,l2]+mH2,n4},{-sp[l1-l2,l1-l2]+mH2,n5}},{l1,l2},{q},{}];
intni = CombinePropagators[intni,{4,3}];
intni = IntegrateLoop[intni,{2,3}];
intni = NormalizePropagator[intni,{1,3}];
intni = CombinePropagators[intni,{3,2}];
intni = CombinePropagators[intni,{2,1}];
intni = IntegrateLoop[intni,{1,1}]