(* ::Package:: *)

(*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] *)
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}]]