Off[General::spell];
Off[General::spell1];

================================
Load FeynCalc:

<< HighEnergyPhysics`fc`

FeynCalc 4.1.0.3b     Evaluate ?FeynCalc for help or visit www.feyncalc.org

?FeynCalc

================================
Input Raw Amps From FeynArts:

rawamps=(

FeynAmpList[Model -> {SM}, Generation1 -> True, 
    Generation2 -> True, Generation3 -> True, 
    ElectronHCoupling -> True, LightFHCoupling -> True, 
    QuarkMixing -> False, UnitaryGauge -> False, 
    RemoveEmptyTops -> True, ProcessName -> V1F1V1F1, 
    Process -> 
     {{V[1], p1, 0}, {F[1], p2, ME}} -> 
      {{V[1], k1, 0}, {F[1], k2, ME}}][FeynAmp[GraphName[V1F1V\
      1F1, T1, I1, N1], EL^2*
     Conjugate[PolarizationVector[k1, λ4]]*
     LeptonSpinor[k2, ME] . DiracMatrix[λ2] . 
      (ME + DiracSlash[-k1 + p2]) . DiracMatrix[λ4] . 
      LeptonSpinor[p2, ME]*PolarizationVector[p1, λ2]*
     PropagatorDenominator[k1 - p2, ME]], 
   FeynAmp[GraphName[V1F1V1F1, T2, I1, N2], 
    EL^2*Conjugate[PolarizationVector[k1, λ4]]*
     LeptonSpinor[k2, ME] . DiracMatrix[λ4] . 
      (ME + DiracSlash[p1 + p2]) . DiracMatrix[λ2] . 
      LeptonSpinor[p2, ME]*PolarizationVector[p1, λ2]*
     PropagatorDenominator[p1 + p2, ME]]]     
);

amps={rawamps[[1,2]],rawamps[[2,2]]}

================================
Conugate Amps

camps=ComplexConjugate[amps]

================================
Square Amps

msquared1=Outer[Times,amps,camps];
Dimensions[msquared1]

================================
Set ME=0 for simplicity:

ME=0;

================================
Contract Polarization Vectors: 
e[mu] e*[nu] -> g[mu,nu]

Clear[ContractBosons];
ContractBosons={

  Pair[LorentzIndex[ComplexIndex[li1_]], 
          Momentum[Polarization[p_, _]]]* 
  Pair[LorentzIndex[             li2_] , 
          Momentum[Polarization[p_, _]]]
   :> MetricTensor[ComplexIndex[li1],li2]
  
 };

================================
Evaluate Square Amps

msquared1

tmp1= msquared1  //.ContractBosons

tmp2 =  tmp1  // Contract

tmp3 =  tmp2  // Map[FermionSpinSum, #, {2}] &

Be careful not to turn on DiracTraceEvaluate  until the FermionSpinSum is done.

(* Turn  DiracTraceEvaluate Off *)
SetOptions[DiracTrace, 
          DiracTraceEvaluate -> 
            True];
SetOptions[DiracTrace, 
        DiracTraceEvaluate -> False];

(* Turn  DiracTraceEvaluate On *)
SetOptions[DiracTrace, 
        DiracTraceEvaluate -> False]; 
      SetOptions[DiracTrace, DiracTraceEvaluate -> True];

msquared2 = tmp3 // PropagatorDenominatorExplicit

================================
Examine Mandelstam Variables

SetMandelstam[s,t,u,p1,p2,-k2,-k1,q,0,ME,ME];

( {ScalarProduct[p1+p2,p1+p2]
  ,ScalarProduct[p1-k2,p1-k2] 
  ,ScalarProduct[p1-k1,p1-k1]}
      //ExpandScalarProduct 
      //Simplify
   )

================================
Evaluate Square Amps

q /: q^2 =-Q^2;

msquared2=msquared2 //ExpandScalarProduct //Simplify

msquared3= 
 TrickMandelstam[ msquared2, {s,t,u,-Q^2+2 ME^2}]

msquared4= msquared3 /(8 EL^4)  /.{ME->0}

msquared5= Plus @@ Plus @@ msquared4 //Simplify