(* ::Package:: *)

(*add properties to the scalar product*)
ScalarProduct[xx_,yy_+zz_] := ScalarProduct[xx,yy]+ScalarProduct[xx,zz];
ScalarProduct[xx_,aa_*yy_/;NumberQ[aa] ] := aa*ScalarProduct[xx,yy];
ScalarProduct[xx_,0] := 0;

log = 0;(*if 0 no log terms in the integrals, if 1 ok*)
l4 = 1;(*if 0 no \[CapitalLambda]^4 terms in the integrals, if 1 ok*)

SimplifyAmpList::usage="Simplify the result return by FeynArts (a FeynAmpList) and compute the integrals for tree level 
or one loop diagrams involving only scalars (TL simplify the result if #2==1)"
(*Begin["Private`"]*)
SimplifyAmpList:=(Module[{k,i,result,gen,result2},
SimplifyAmpList::mt1l="More than one loop, not implemented";result=0;
(*run over the FeynAmps*)
For[k=1,k<=Length[#],k++,
result2=0;
If[Length[#[[k,2]]]==0,gen = #[[k,3]];,
(*integrate if one loop*)
If[Length[#[[k,2]]]==1,gen=LoopIntegrate[#[[k,3]],#[[k,2,1]]];,Message[SimplifyAmpList::mt1l]](*end if*)
](*end if*);
(*make the different replacement du to the different particles in the loop*)
For[i=1,i<=Length[#[[k,4,2]]],i++,result2 = result2+(gen)/.MakeReplaceList[#[[k,4,1]],#[[k,4,2,i]]];
];(*end for*)
If[#2==1,
result2 = TLSimplify[result2/.M$FACouplings];,result2 = Simplify[result2/.M$FACouplings];];
result = result+result2;
];(*end for*)
Simplify[result]
](*end module*))&(*end SimpifyAmpList*)


LoopIntegrate::usage="takes a generic one loop amplitude FeynAmp[[3]] and the momentum in the loop FeynAmp[[2,1]] 
and return the result of the integration"
LoopIntegrate:=(Module[{n,res,temp,dP,D},
LoopIntegrate::tlnp="Too large number of propagtor(not implemented)";
(*number of propagators<=2*)
If[Length[#[[4]]]<=2,
(*number of propagator=2*)
If[Length[#[[4]]]==2,
res=0;
temp = (#[[5]]*#[[6]]);
ScalarProduct[xx_,x*yy_] := x*ScalarProduct[xx,yy];
dP = Simplify[(Coefficient[#[[4,1,1]]^2*(1-x)/2+#[[4,2,1]]^2*x,#2,1]/2)];
D = Simplify[(ScalarProduct[Coefficient[#[[4,1,1]],#2,0],Coefficient[#[[4,1,1]],#2,0]]+#[[4,1,2]]^2)*(1-x)
+(ScalarProduct[Coefficient[#[[4,2,1]],#2,0],Coefficient[#[[4,2,1]],#2,0]]+#[[4,2,2]]^2)*x-ScalarProduct[dP,dP]];
(*Check change of variable : Print[(ScalarProduct[#[[4,1,1]],#[[4,1,1]]]+#[[4,1,2]]^2)*(1-x)+(ScalarProduct[#[[4,2,1]],#[[4,2,1]]]+#[[4,2,2]]^2)*x];
Print[Simplify[(((ScalarProduct[#[[4,1,1]],#[[4,1,1]]]+#[[4,1,2]]^2)*(1-x)+(ScalarProduct[#[[4,2,1]],#[[4,2,1]]]+#[[4,2,2]]^2)*x)/.{#2->#2-dP})-D]];*)
temp = temp/.{#2->#2-dP};temp=ExpandAll[temp];
For[n=1,n<=Length[temp],n++,
Switch[temp[[n]],
_*ScalarProduct[#2,#2]*ScalarProduct[#2,#2],
(*Print[temp[[n]]];*)res = res+(l4*\[CapitalLambda]^4/2-2*Integrate[((D)/.ScalarProduct[xx_,yy_]->-ScalarProduct[xx,yy]),{x,0,1}]*\[CapitalLambda]^2)*temp[[n]]/ScalarProduct[#2,#2]/ScalarProduct[#2,#2];
(*Print[(l4*\[CapitalLambda]^4/2-2Integrate[D,{x,0,1}]*\[CapitalLambda]^2)*temp[[n]]/ScalarProduct[#2,#2]/ScalarProduct[#2,#2]];*),
_*ScalarProduct[#2,#2]*ScalarProduct[#2,_],
(*Print[temp[[n]]];*)res = res-(2*\[CapitalLambda]^2/d)*Integrate[(temp[[n]]/ScalarProduct[#2,#2])/.ScalarProduct[#2,xx_]->ScalarProduct[dP,xx],{x,0,1}];
(*Print[(2*\[CapitalLambda]^2/d)*Integrate[(temp[[n]]/ScalarProduct[#2,#2])/.ScalarProduct[#2,xx_]->ScalarProduct[dP,xx],{x,0,1}]];*),
_*ScalarProduct[#2,#2],
res = res-(\[CapitalLambda]^2)*Integrate[(temp[[n]]/ScalarProduct[#2,#2]),{x,0,1}];(*-sign, wick*)
(*Print[temp[[n]]];Print[-(\[CapitalLambda]^2)*Integrate[(temp[[n]]/ScalarProduct[#2,#2]),{x,0,1}]];*),
_*ScalarProduct[#2,_]*ScalarProduct[#2,bof_],
res = res+(-\[CapitalLambda]^2/d)*Integrate[(temp[[n]])/.{ScalarProduct[#2,xx_]*ScalarProduct[#2,yy_]->ScalarProduct[yy,xx]},{x,0,1}];
(*Print[temp[[n]]];Print[(-\[CapitalLambda]^2/d)*Integrate[(temp[[n]])/.ScalarProduct[#2,xx_]*ScalarProduct[#2,yy_]->ScalarProduct[yy,xx],{x,0,1}]];*),
_*ScalarProduct[#2,_]*ScalarProduct[#2,_],
res = res+(-\[CapitalLambda]^2/d)*Integrate[(temp[[n]])/.{ScalarProduct[#2,xx_]->(ScalarProduct[xx,xx])^(1/2)},{x,0,1}];
](*end switch*)
];(*end for*)
res=(\[ImaginaryI]*\[Pi]^2*res),
(*number of propagator=1*)
res=0;For[n=1,n<=Length[#[[5]]],n++,
Switch[#[[5,n]],(*on the kinematic factor of the numerator*)
_*ScalarProduct[#2,#2],
res = res+#[[5,n]]/ScalarProduct[#2,#2]*\[ImaginaryI] \[Pi]^2(l4*(\[CapitalLambda]^4)/2-#[[4,1,2]]^2*\[CapitalLambda]^2+#[[4,1,2]]^4*log*Log[1+\[CapitalLambda]^2/#[[4,1,2]]^2]);,
_*ScalarProduct[#2,_],
res =res;,
_,
res=res+#[[5,n]](-\[ImaginaryI]) \[Pi]^2(\[CapitalLambda]^2-#[[4,1,2]]^2*log*Log[1+\[CapitalLambda]^2/#[[4,1,2]]^2])];](*end for*)];(*end if*),
Message[LoopIntegrate::tlnp]
];(*end if*)
res*#[[1]]*#[[2]]*#[[3]]](*end module*))&(*end LoopIntegrate*)

MakeReplaceList::usage="takes 2 lists of the same length, l1 and l2 and return a list of replacement L[[i]]=l1[[i]]->l2[[i]]"
MakeReplaceList:=(Module[{i,Rep},MakeReplaceList::dll="The lists are not of the same length";
If[Length[#]==Length[#2],
Rep = {};
For[i=1,i<=Length[#],i++,Rep = Append[Rep,#[[i]]-> #2[[i]]];];Rep,
Message[MakeReplaceList::dll]](*end if*)
](*end module*))&(*end MakeReplaceList*)

(*End[]
EndPackage[]*)


TLSimplify::usage="return the expression after simplification using the tree level relation 
for the masses and the parameters m0, rmu, rmd, rms in the isospin limit"
TLSimplify:=(Module[{res},
m0=Sqrt[2/3(mk^2-mp^2)(1-2Sqrt[2]/Tan[2T])];
mup=mp^2/r;
md=mp^2/r;
ms = (2mk^2-mp^2)/r;
res = Simplify[#/.{MK0 ->mk, MKm->mk,Mpi0->mp,Mpim->mp,Meta->Sqrt[1/3(4mk^2-mp^2+2Sqrt[2](mk^2-mp^2)Tan[T])],
Metap->Sqrt[1/3(4mk^2-mp^2-2Sqrt[2](mk^2-mp^2)Cot[T])]}];
Clear[m0,mup,md,ms];
res
](*end module*))&(*end TLSimplify*)


MomSimplify::usage="return the expression after simplification using momentum conservation and p^2=m^2 for the process #2(FeynAmpList[[0,8,2]])"
MomSimplify:=(Module[{temp},
temp = #/.MomConservation[Join[-#2[[1]],#2[[2]]]];temp = Simplify[temp/.MomConservation[Join[-#2[[1]],#2[[2]]]]];(*Print[MomConservation[Join[-#2[[1]],#2[[2]]]]];*)
temp = Simplify[temp/.MakeMomentumRL[Join[-#2[[1]],#2[[2]]]]];temp
](*end Module*)
)&(*end MomSimplify*)

MakeMomentumRL::usage="For a process, make the replacement list, p^2->m^2"
MakeMomentumRL:=(Module[{list,k},
list = {};
For[k=2,k<=Length[#],k++,list = Append[list,ScalarProduct[#[[k,2]],#[[k,2]]]->#[[k,3]]^2];
];(*end for*)
list
](*end module*)
)&(*end momentumRL*)

MomConservation::usage="p1->sum of outgoing momentum-sum of ingoing momentum+p1 and p1^2=m1^2"
MomConservation:=(Module[{p1,k},
p1=0;
For[k=2,k<=Length[#],k++,p1 = p1 +#[[k,2]];];
If[Length[#]>=3,
{-#[[1,2]]->p1,ScalarProduct[#[[Length[#],2]],#[[Length[#]-1,2]]]->Simplify[#[[1,3]]^2/2-ScalarProduct[p1,p1]/2+ScalarProduct[#[[Length[#],2]],#[[Length[#]-1,2]]]]},
-#[[1,2]]->p1
](*end if*)

](*end module*)
)&(*end MomConservation*)
