$LoadFeynArts=True;
<<FeynCalc`

$FAVerbose=0;

SetOptions[InitializeModel,ModelEdit:>(M$ClassesDescription=M$ClassesDescription/.MZ->MZc)];
InitializeModel[{SM,UnitarySM},GenericModel->{Lorentz,UnitaryLorentz}];

NoElectronHCoupling=ExcludeFieldPoints->{FieldPoint[0][-F[2,{1}],F[1,{1}],S[3]],FieldPoint[0][-F[2,{1}],F[2,{1}],S[1]],FieldPoint[0][-F[2,{1}],F[2,{1}],S[2]]};

part=InsertFields[CreateTopologies[0,2->3],{F[2,{1}],-F[2,{1}]}->{F[3,{3}],-F[3,{3}],S[1]},Restrictions->NoElectronHCoupling,InsertionLevel->{Classes},Model->{SM,UnitarySM},GenericModel->{Lorentz,UnitaryLorentz}];
Paint[part,PaintLevel->{Classes}];

listdiag=FCFAConvert[CreateFeynAmp[part],SMP->True,ChangeDimension->4,IncomingMomenta->{p1,p2},OutgoingMomenta->{k1,k2,k3},DropSumOver->True,List->False,UndoChiralSplittings->True];
diag1=listdiag[[5]]
diag2=listdiag[[4]]
diag3=listdiag[[3]]
diag4=listdiag[[2]]
diag5=listdiag[[1]]

ClearScalarProducts;

SP[p1,p1]=me^2;
SP[p2,p2]=me^2;
SP[p1,p2]=(s-2 me^2)/2;

SP[k1,k1]=mt^2;
SP[k2,k2]=mt^2;
SP[k3,k3]=mh^2;
SP[k1,k2]=(s3-2 mt^2)/2;
SP[k1,k3]=(s2-mt^2-mh^2)/2;
SP[k2,k3]=(s1-mt^2-mh^2)/2;

diag1=diag1/.{Lor1->μ,Lor2->ν,Lor3->α,Lor4->β};
diag1C=ComplexConjugate[diag1]/.{μ->μlin,ν->νlin,α->αlin,β->βlin};

M12=diag1 diag1C//Contract

res=((M12//FermionSpinSum)/.DiracTrace->Tr)//Contract;

FreeQ2[res,{DiracTrace,LorentzIndex}]

res