plotting – 1D Integral always evaluates to 0 and Replace All

I am a bit puzzled by Mathematica’s behaviour. I am trying to evaluate some 1D integrals which depend on parameters phi1,phi2 numerically (functions f1 and f2), however NIntegrate always outputs 0 whereas plotting the integrand with Plot shows that the integral is non-zero.

Here is the code in question:

t = 1;
lambda = 0.3;
hbarOmega = 4;
q = 0.5;
alpha = 4 lambda^2 t (Sin(q/2))/q;

A(k_) := 2 t (Sin(q/2)/q) Sin(k) (1/(hbarOmega + 4 t Sin(q/2) Sin(k)) + 
 1/(hbarOmega - 4 t Sin(q/2) Sin(k)));

B(k_) := 2 t (Sin(q/2)/q) Sin(k) (1/(hbarOmega + 4 t Sin(q/2) Sin(k)) - 
 1/(hbarOmega - 4 t Sin(q/2) Sin(k)));

e1(k_) := 4 t lambda^2  (Sin(q/2)/q) B(k + q/2) Sin(k + q/2);
e2(k_) := -4 t lambda^2  (Sin(q/2)/q) B(k - q/2) Sin(k - q/2);
en(k_) := -2 t Cos(k) + e1(k - q) + e2(k + q);
F(k_, phi1_, phi2_) := - alpha (A(k) phi1 + Sin(k) phi2);
mu(phi1_, phi2_) := - 2 alpha Re(phi1 Conjugate(phi2));

hMf(k_, phi1_, phi2_) := {{en(k + q/2),F(k, phi1, phi2)}, {Conjugate(F(k, phi1, phi2)), en(k - 
q/2)}};


u(k_, phi1_, phi2_) := Transpose(Eigenvectors(hMf(k, phi1, phi2)));
enMf(k_, phi1_, phi2_) := ConjugateTranspose(u(k, phi1, phi2)).hMf(k, phi1, phi2).u(k, phi1, 
phi2);

o(k_, phi1_, phi2_) := ConjugateTranspose(u(k, phi1, phi2)).{{0, 0}, {1, 0}}.u(k, phi1,phi2);


f1(phi1_?NumericQ,phi2_?NumericQ) := (1/(2 Pi)) NIntegrate(
Sin(k)*(o(k, phi1, phi2)((1, 1)) Boole(
     enMf(k, phi1, phi2)((1, 1)) <= mu(phi1, phi2)) + 
   o(k, phi1, phi2)((2, 2)) Boole(
     enMf(k, phi1, phi2)((2, 2)) <= mu(phi1, phi2))), {k, -Pi, Pi});

f2(phi1_?NumericQ,phi2_?NumericQ) := (1/(2 Pi)) NIntegrate(A(k) (o(k, phi1, phi2)((1, 1)) 
Boole(enMf(k, phi1, phi2)((1, 1)) <= mu(phi1, phi2)) + 
   o(k, phi1, phi2)((2, 2)) Boole(
     enMf(k, phi1, phi2)((2, 2)) <= mu(phi1, phi2))), {k, -Pi,Pi});

f1(1, 1)

Plot(Sin(k)*(o(k, phi1, phi2)((1, 1)) Boole(
   enMf(k, phi1, phi2)((1, 1)) <= mu(phi1, phi2)) + 
 o(k, phi1, phi2)((2, 2)) Boole(
   enMf(k, phi1, phi2)((2, 2)) <= mu(phi1, phi2))) /. {phi1 -> 1, phi2 -> 1}, {k, -Pi, Pi})

Plot(Sin(k)*(o(k, 1, 1)((1, 1)) Boole(enMf(k, 1, 1)((1, 1)) <= mu(1, 1)) + o(k, 1, 1)((2, 2)) 
Boole(enMf(k, 1, 1)((2, 2)) <= mu(1, 1))), {k, -Pi, Pi})

So for example, when I evaluate the integral given by f1(1,1), I get an output of 0 alongside with a warning that Mathematica couldn’t find all eigenvectors of a 2×2 matrix which is bizarre, because when I plot the eigenvalues separately with:

Plot({enMf(k, 1, 1)((1, 1)), enMf(k, 1, 1)((2, 2))}, {k, -Pi, Pi})

the warning is not present.

However when I plot the integrand inside f1 with the same parameters phi1=1,phi2=1 with:

Plot(Sin(k)*(o(k, 1, 1)((1, 1)) Boole(
  enMf(k, 1, 1)((1, 1)) <= mu(1, 1)) + 
o(k, 1, 1)((2, 2)) Boole(
  enMf(k, 1, 1)((2, 2)) <= mu(1, 1))), {k, -Pi, Pi})

it shows that the integrand is always positive and that the integral should be non-zero.

Moreover, there is another oddity: replacing the parameters phi1,phi2 within the integrand with /. for the same values phi1=1,phi2=1 like so:

Plot(Sin(k)*(o(k, phi1, phi2)((1, 1)) Boole(
   enMf(k, phi1, phi2)((1, 1)) <= mu(phi1, phi2)) + 
 o(k, phi1, phi2)((2, 2)) Boole(
   enMf(k, phi1, phi2)((2, 2)) <= mu(phi1, phi2))) /. {phi1 -> 1,phi2 -> 1}, {k, -Pi, Pi})

gives a different plot compared to just substituting the values phi1=1,phi2=1 directly:

Plot(Sin(k)*(o(k, 1, 1)((1, 1)) Boole(
  enMf(k, 1, 1)((1, 1)) <= mu(1, 1)) + 
o(k, 1, 1)((2, 2)) Boole(
  enMf(k, 1, 1)((2, 2)) <= mu(1, 1))), {k, -Pi, Pi})

which I find very bizarre and assume has to do with delayed evaluation. Any suggestions would be appreciated.