How to solve the roundoff error problem in scipy.integrate.quad on specific integrand?

Question

I'm doing the following numerical calculation.

import numpy as np
from scipy.integrate import quad

def E_k(gap, q, mu):
    return np.sqrt(abs(gap) * abs(gap) + (q * q - mu) * (q * q - mu))

quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
                 0, np.inf, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)

It is aimed to set the value of mu to be in the range of (-30, 1), the value of gap to be in the range of (0, 3), and T to be positive. The quad will give wrong results (it tells me roundoff error detected) when T is small (approximately less than 0.1) and mu is positive or negatively close to zero (approximately -0.0005 < mu < 1).

I have tried to divide the integrating interval at point mu when mu is positive.

quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
                 0, mu, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)[0] + \
quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
                 mu, np.inf, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)[0]

But it seems not to work. I have also increased the limit and adjusted epsabs, epsrel in quad, but the results are still uncorrected.