Trouble computing integral with scipy quad

Question

I'm trying to compute the following definite integral:

Integral I want to compute:

where rho_ch is

and

a = 3.66 * 10^(-15) m (in meters)

b = 0.54 * 10^(-15) m

and rho_0 = 1.23 * 10^(-35) C/m^3.

When I compute it, I get something of the order of 10^(-61), but I believe it should be something closer to 1.

I take the upper limit of the integral to 10^(-10) because the integral should have converged by then. My guess is that the problem has do with overflow in the rho_ch function, although I tried a fix for that, but it didn't work.

Does anyone have any idea of what I'm doing wrong?

import numpy as np
from scipy.constants import *
from scipy import integrate

Z_t = 20
a = 1.07*40**(1/3)
b = 0.54
rho_0 = 0.077*e
X_0 = np.array([rho_0, a, b])*10**(-15)
E = 250*10**(6)*e
p_sq = (E/c)**2-(m_e*c)**2
beta_sq = 1-(1/(E/(m_e*c**2))**2)
theta = np.pi/6

def q(theta):
    return 2*np.sqrt(p_sq*np.sin(theta/2)**2)

def F_quad(theta, X):
    return 4*np.pi*hbar/(Z_t*e*q(theta))*integrate.quad(integrand2, 0, 10**(-10), args=(theta, X))[0]

integrand2 = lambda r, theta, X: r*X[0]/(1+np.exp((r-X[1])/X[2]))*np.sin(q(theta)*r/hbar) if (r-X[1])/X[2] <= 600 else 0