(********************* Program QEDVertex.m ***********************)
(*
 Version campatible with FeynCalc 10.1.0

 Date: 10/12/2024
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)
 
(* First input FeynCalc *)
(* Uncomment below if you want to call from this program. If open a new
   mathematica notebook and load FeynCalc from there you should not load
   it again
*)
(*
<< FeynCalc`
*)

(* Tell FeynCalc to reduce the result to scalar functions *)

SetOptions[PaVeReduce,A0ToB0->True]
$LimitTo4 = True;

(* Useful Function *)

TakeDTo4 = Function[exp, expaux1 = exp /. D -> 4 - eps; 
  expaux2 = Normal[Series[expaux1, {eps, 0, 1}]];
  c0 = Coefficient[expaux2, eps, 0]; c1 = Coefficient[expaux2, eps, 1];
  c1div = c1 /. PaVe[0, {z1_}, {z2_, z3_}] -> 2/eps;
  expaux3 = c0 + eps c1div // Simplify;
  Simplify[expaux3 /. eps -> 0]]


(* Now write the numerator of the Feynman diagram. We define the 
   constant

       C=  alpha/(4 pi) 

   The kinematics is: q = p1 -p2 and the internal momenta is k.

*)

num:=SpinorUBarD[p1,m].GAD[ro].(GSD[p1]-GSD[k]+m).GAD[mu].(GSD[p2]-GSD[k]+m).GAD[ro].SpinorUD[p2,m]


  amp:=C num FAD[{k,lbd},{k-p1,m},{k-p2,m}]

(* Define the on-shell kinematics *)
  
onshell={ScalarProduct[p1,p1]->m^2,ScalarProduct[p2,p2]->m^2, \
         ScalarProduct[p1,p2]->m^2-q2/2}

(* Define the divergent part of the relevant PV functions*)

div={PaVe[0,{a_},{b_,c_}]->Div}

res1=(-I / Pi^2) OneLoop[k,amp] 

res=res1 /. onshell

auxV1= res /.onshell
auxV2= PaVeReduce[auxV1]
auxV3= PaVeReduce[auxV2] /. div
divV=Simplify[Div*Coefficient[auxV3,Div]]

(* Check  that the divergencies do not cancel *)

testdiv:=Simplify[divV]

ans1=res;
ola=Select[Variables[ans1],(Head[#]===StandardMatrixElement)&]
Print["StandardMatrixElements\n",ola]
var=Select[Variables[ans1],(Head[#]===StandardMatrixElement)&]
  Set @@ {var, {ME[1],ME[2],ME[3],ME[4]}}

(* Extract the different Matrix Elements 

Mathematica writes the result in terms of 4 non independent
Standard Matrix Elements. To have a simpler result we substitute
these elements by simpler expressions (ME[1],ME[2],ME[3],ME[4]). 

PR=GA[6]
PL=GA[7]		 
			 
 {StandardMatrixElement[u[p1, m1] . PR. u[p2, m2]],

  StandardMatrixElement[u[p1, m1] . PL. u[p2, m2]], 
 
  StandardMatrixElement[u[p1, m1] . ga[mu] . PR . u[p2, m2]],

  StandardMatrixElement[u[p1, m1] . ga[mu] . PL . u[p2, m2]]}

*)

(*

ME[5] = StandardMatrixElement[u[p1, m1] .  u[p2, m2]]}

ME[6] = StandardMatrixElement[u[p1, m1] . GA[mu].  u[p2, m2]]}  

*)

ans2=PaVeReduce[PaVeReduce[ans1]]/. {ME[1]->ME[5]-ME[2],ME[3]->ME[6]-ME[4]}//FCE//Simplify

  
CE5=Coefficient[ans2, ME[5]] 
CE6=Coefficient[ans2, ME[6]] 
CE51=Coefficient[CE5,FV[p1,mu]]
CE52=Coefficient[CE5,FV[p2,mu]]
  
ans3=CE51 (FV[p1,mu]+FV[p2,mu]) ME[5] +  CE6 ME[6]

test1:=Simplify[CE51-CE52]
test2:=Simplify[ans2-ans3]

ans4= ans3 /. {(FV[p1,mu]+FV[p2,mu]) ME[5] -> 2 m ME[6] -2m ME[7]}

ans5=TakeDTo4[ans4]
    
CGamma:=Coefficient[ans5,ME[6]]
CSigmaAux:=Coefficient[ans5,ME[7]]

test3:=Simplify[ans5-CGamma ME[6] -CSigmaAux ME[7]]

  
F2:=CSigmaAux /. lbd->0//Simplify 

delZ1aux= - CGamma /. q2->0 //Simplify

delZ1:= delZ1aux /. lbd->0//Simplify

F1:=CGamma + delZ1 /. lbd->0 //Simplify

Print["delZ1=",delZ1,"\n"]
Print["F1=",F1,"\n"]
Print["F2=",F2,"\n"]    
Print["F1[0]=",F1/.q2->0,"\n"]  
F2zero=F2/.q2->0;
Print["F2[0]=",F2zero,"\n"]  
F2zeroFinal=F2zero /. {PaVe[0,{m^2},{0,m^2}]-> PaVe[0,{0},{m^2,m^2}]+2,C-> \[Alpha]/(4 Pi)}//Simplify
Print["F2[0]=",F2zeroFinal,"\n"]  


(***************** End of Program QEDVertex.m ********************)
