(* Definitions *)
dm[mu_]:=DiracMatrix[mu,Dimension->D]
ds[p_]:=DiracSlash[p]
mt[mu_,nu_]:=MetricTensor[mu,nu]
fv[p_,mu_]:=FourVector[p,mu]
sp[p_,q_]:=ScalarProduct[p,q]
SetOptions[OneLoop,Dimension->4];


(* Vacuum Polarization *)
cst1=  e^2;
num1:= fv[-2 p + k,mu] fv[-2 p + k,nu];
amp1:=cst1 num1 \
FeynAmpDenominator[PropagatorDenominator[p,m], \
PropagatorDenominator[p-k,m]];
res1=(-I / Pi^2) OneLoop[p,amp1] ;
res1= PaVeReduce[res1];

cst2= - 2 e^2;
num2:= mt[mu,nu];
amp2:=cst2 num2 \
FeynAmpDenominator[PropagatorDenominator[p,m]];
res2=(-I / Pi^2) OneLoop[p,amp2] ;
res2= PaVeReduce[res2] ;
res=res1 + res2

aux1=Pair[LorentzIndex[mu],Momentum[k]]*Pair[LorentzIndex[nu], Momentum[k]];
aux2=Pair[LorentzIndex[mu], LorentzIndex[nu]];
PI1= Coefficient[res,aux1]/I;
PI2= Coefficient[res,aux2]/I;
CheckTransverse:=Simplify[PI1 - (-PI2/sp[k,k])]

Div={C0[x1_,x2_,x3_,x4_,x5_,x6_]->0,B0[k2_,m02_,m12_]-> 2 div I/(16 Pi^2), 
     A0[m2_]-> m^2 2 div I/(16 Pi^2)}
DivPI =Expand[PI1 /. Div]
deltaZ3 = - Coefficient[DivPI ,div] /. sp[k,k]->0 


(* Scalar Self Energy *)
cst1= e^2;
cst2= - 2 e^2 (1/2) ;
numself1 = fv[2 k +p,mu] (- mt[mu,nu] ) fv[2 k +p,nu];
numself2 =  (- mt[mu,nu] ) mt[mu,nu];

ampself1= cst1 numself1 \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[p+k,m]];
resself1=(-I / Pi^2) OneLoop[p,ampself1] ;
resself1= PaVeReduce[resself1];

ampself2= cst2 numself2 \
FeynAmpDenominator[PropagatorDenominator[p,0]];
resself2=(-I / Pi^2) OneLoop[p,ampself2] ;
resself2= PaVeReduce[resself2];

Sigma = 1/(-I) ( resself1 + resself2 )

DivSigma = Expand[Sigma /. Div] ;
DivSigma = DivSigma /. sp[k,k]->k2

deltaZ2 = D[DivSigma,k2]


(* Vertex *)
cstvertex1 = e^2 ;
numvertex1 = fv[-2 q1 -p,a] (-mt[a,b]) fv[2 q2-p,b] fv[q2-q1-2 p,mu];

ampvertex1 = cstvertex1 numvertex1 \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[p+q1,m],PropagatorDenominator[q2-p,m]];
resvertex1=(-I / Pi^2) OneLoop[p,ampvertex1] ;
resvertex1= PaVeReduce[resvertex1];

cstvertex2 = 2 e^2 (-1) ;
cstvertex3 = 2 e^2 (-1) ;
numvertex2 = mt[mu,a] (-mt[a,b]) fv[-2 q1 +p,b] ;
numvertex3 = mt[mu,a] (-mt[a,b]) fv[2 q2 -p,b] ;

ampvertex2 = cstvertex2 numvertex2 \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[q1-p,m]] ;
resvertex2=(-I / Pi^2) OneLoop[p,ampvertex2] ;
resvertex2= PaVeReduce[resvertex2];

ampvertex3 = cstvertex3 numvertex3 \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[q2-p,m]] ;
resvertex3=(-I / Pi^2) OneLoop[p,ampvertex3] ;
resvertex3= PaVeReduce[resvertex3];
resvertex = ( resvertex1 + resvertex2 + resvertex3 )

DivVertex = Expand[resvertex /. Div];
DivVertex = DivVertex /. sp[q1, q2] -> (k2 - sp[q1, q1] - sp[q2, q2])/2
DivVertex = Simplify[DivVertex] 
aux =  Pair[LorentzIndex[mu], Momentum[q2]];
  
deltaZ1 =  -I Coefficient[DivVertex, aux] 


(* Check Ward Identies *)
Print["CheckWardIdentities:","  deltaZ1-deltaZ2 = ",deltaZ1-deltaZ2]


(* Beta Function for QED *)
deltaZe= -1/2 deltaZ3 ;
beta= 1/2 e^2 D[deltaZe,e];
b1QED= beta/e^3 16 Pi^2 ;

Print["Final Result:"," b1ScalarQED = ",b1QED]



