Print[" ** W -> L n ** "]
MW=.
ML=.
EL=.
PxL=.
PyL=.
PzL=.
EN=.
PxN=.
PyN=.
PzN=.
PtN=.
QLN2 = ML^2 +2EL*EN - 2(PxL*PxN+PyL*PyN + PzL*PzN);
Print["  MW2 = ",QLN2];

MET=.
ϕN=.

EN = Sqrt[MET^2+PzN^2];

QLN2=FullSimplify[QLN2];
Print["  MW2 = ",QLN2]

Eq1 = QLN2-MW^2;
solucio = FullSimplify[Solve[Eq1 ==0,PzN]];
PzN1=PzN/.solucio[[1]];
PzN2=PzN/.solucio[[2]];
Print["  PzN1 = ",PzN1];
Print["  PzN2 = ",PzN2];

a=Apart[(EL-PzL) (EL+PzL)];
b = -(-(ML^2)+(MW^2)+2 PxL PxN+2 PyL PyN) PzL;
Print["a = ",a]
Print["b = ",b]
Print["-b/a = ",-b/a]
Print["PzN1+PzN2 = ",FullSimplify[PzN1+PzN2]]
Print["PzN1-PzN2 = ",FullSimplify[PzN1-PzN2]]
Discrim=FullSimplify[(PzN1-PzN2)*a]^2;
Print["Sqrt_Term = ",Discrim]
Print["c = ",FullSimplify[(Discrim-b^2)/(-4a)]]

PxN = MET*Cos[ϕN];
PyN = MET*Sin[ϕN];
(*solDelta = Together[FullSimplify[Solve[Discrim==0,MET]]];*)
Print["Sqrt_Term = ",Discrim]
solDelta = FullSimplify[Solve[Discrim==0,MET]];
solDeltaN1=MET/.solDelta[[1]];
solDeltaN2=MET/.solDelta[[2]];
Print["solDeltaN1 = ",solDeltaN1]
Print["solDeltaN2 = ",solDeltaN2]