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= (-1) e2
num1:=Tr[dm[mu] . (ds[p] +m) . dm[nu]. ( ds[p] -ds[k] + m)]
amp1:=cst1 num1 \
FeynAmpDenominator[PropagatorDenominator[p,m], \
PropagatorDenominator[p-k,m]]

res1=(-I / Pi^2) OneLoop[p,amp1] 
res1= PaVeReduce[res1]

aux=Pair[LorentzIndex[mu],Momentum[k]]*Pair[LorentzIndex[nu], Momentum[k]] - 
    Pair[LorentzIndex[mu], LorentzIndex[nu]]* Pair[Momentum[k], Momentum[k]]

PI= Coefficient[res1,aux]/I
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[PI /. Div]
deltaZ3= Coefficient[-DivPI ,div] /. sp[k,k]->0 
Print["deltaZ3 = ",deltaZ3]


(* Electron Self Energy *)
cstself= e^2;
numself = dm[mu] . (ds[p] + ds[k] +m ) . dm[nu] (-mt[mu,nu])

ampself= cstself numself \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[p+k,m]];
resself=(-I / Pi^2) OneLoop[p,ampself] ;
resself= PaVeReduce[resself];

Sigma = 1/(-I) resself

DivSigma = Expand[Sigma /. Div] ;
DivSigma = DivSigma /. sp[k,k]->k2
aux = DiracGamma[Momentum[k]]
deltaZ2 = Coefficient[DivSigma,aux]


(* Vertex *)
cstvertex = e^2 ;
numvertex = dm[a].(ds[q1]-ds[p]+m).dm[mu].(-ds[q2] -ds[p]+m).dm[b] (-mt[a,b])

ampvertex = cstvertex numvertex \
FeynAmpDenominator[PropagatorDenominator[p,0], \
PropagatorDenominator[q1-p,m],PropagatorDenominator[-q2-p,m]];
resvertex=(-I / Pi^2) OneLoop[p,ampvertex] ;
resvertex= PaVeReduce[resvertex];

DivVertex = Expand[resvertex /. Div];
DivVertex = Coefficient[DivVertex,div] div //Simplify ;
aux =  DiracGamma[LorentzIndex[mu]];

deltaZ1 =  -I Coefficient[DivVertex, aux] 


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

(* Beta Function for QED *)
deltaZe= -1/2 deltaZ3 /. e2->e^2;
beta= 1/2 e^2 D[deltaZe,e] ;
b1QED= beta/e^3 16 Pi^2; 
Print["Final Result:"," b1QED = ",b1QED]
