(* ::Package:: *)

(* Obtain Feynman parameter integral *)
FP = FeynmanParametrization[{{l1,0,a1},{l1+p1,0,a2},{l1+p1+p2,0,a3},{l1+p3,0,a4}},{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}]


(* Jacobian of variable transform *)
With[{xi={u1*v1,u1*(1-v1),u2*v2,u2*(1-v2)},yzj={u1,u2,v1,v2}},
Table[D[xi[[i]],yzj[[j]]],{i,1,4},{j,1,4}]]
Det[%]


(* Variable transform *)
FP = Simplify[u1*u2*FP /. {x[1]->u1*v1,x[2]->u1*(1-v1),x[3]->u2*v2,x[4]->u2*(1-v2)},Assumptions->u1>0&&u2>0&&0<v1<1&&0<v2<1]
(* Delta function *)
FP = Simplify[Simplify[FP /. {u2->1-u1},Assumptions->0<u1<1&&0<v1<1&&0<v2<1&&a1>0&&a2>0&&a3>0&&a4>0]*(
   u1^(2-a1-a2-a3-a4-eps)*(1-u1)^(2-a1-a2-a3-a4-eps)*(-(t (-1+v1) (-1+v2)+s v1 v2))^(2-a1-a2-a3-a4-eps))/(-(1-u1) u1 (t (-1+v1) (-1+v2)+s v1 v2))^(2-a1-a2-a3-a4-eps) /. {
  (u1 v1)^(-1+a1)->u1^(-1+a1)*v1^(-1+a1),(u1-u1 v1)^a2->u1^a2*(1-v1)^a2,((-1+u1) (-1+v2))^a4->(1-u1)^a4*(1-v2)^a4,(v2-u1 v2)^(-1+a3)->(1-u1)^(-1+a3)*v2^(-1+a3)}]
(* u1 integration *)
FP = Integrate[FP,{u1,0,1},GenerateConditions->False]


(* Apply Mellin-Barnes representation to the denominator *)
MB = Simplify[FP*(Gamma[-2+a1+a2+a3+a4+eps+z]Gamma[-z]/(Gamma[-2+a1+a2+a3+a4+eps]))*x^z*(1-v1)^z (1-v2)^z*(-s)^(2-a1-a2-a3-a4-eps)*v1^(2-a1-a2-a3-a4-eps-z)*v2^(2-a1-a2-a3-a4-eps-z)/((-s v1 v2+t (-1+v1+v2-v1 v2))^(2-a1-a2-a3-a4-eps))]
(* Perform v1,2 integration *)
MB = Integrate[MB,{v1,0,1},{v2,0,1},GenerateConditions->False]


(* F1111 *)
MB1111 = MB /. {a1->1,a2->1,a3->1,a4->1}


(* The poles at -1 and -1-eps pinch the contour *)
(* Split in integration along contour with Re[z]=1/2 and -residue at -1-eps  *)


(* residue contribution *)
MB1111res = FullSimplify[-Residue[MB1111,{z,-1-eps}]]
MB1111res = FullSimplify[Series[MB1111res*Exp[eps*EulerGamma]*(-s)^(2+eps)*x,{eps,0,1}]]


(* evaluate integral along straight line by closing the contour to the right *)
MB1111int = FullSimplify[-Residue[Normal[Series[MB1111,{eps,0,1}]],{z,n},Assumptions->n\[Element]Integers&&n>=0],Assumptions->n\[Element]Integers&&n>=0]
MB1111int = Apart[MB1111int/.{n->n-1},n]/.{n->n+1}
MB1111int = Series[Sum[MB1111int,{n,0,Infinity}]*Exp[eps*EulerGamma]*(-s)^(2+eps)*x,{eps,0,1}]


(* Result up to O[eps], prefactor is Exp[-eps*EulerGamma]*(-s)^-eps/(s*t)*)
FullSimplify[MB1111res + MB1111int]


(*Expansion in t/s*)
(*1/t*)
-Residue[MB1111,{z,-1-eps}]
(*t^0*)
-(Residue[MB1111,{z,-eps}] + Residue[MB1111,{z,0}])
(*t*)
-(Residue[MB1111,{z,1-eps}] + Residue[MB1111,{z,1}])


(*Expansion in s/t*)
(*1/s*)
Residue[MB1111,{z,-1}]
(*s^0*)
(Residue[MB1111,{z,-2-eps}] + Residue[MB1111,{z,-2}])
(*s*)
(Residue[MB1111,{z,-3-eps}] + Residue[MB1111,{z,-3}])