Monte Carlo 4D integral with variable limit

Question

So, i have this wonderfull function f(th,k,p,n,W) that i need to integrate over 4 variables, but one of them, the variable k, goes from 0 to p (then i integrate over p). I've tried doing this integration using nquad, but it is taking forever... so i was wondering if i can do this integral using monte carlo, but i don't know how to do this with a variable limit. Does any one know how can i do this ?

import numpy as np
import sympy as sym
from scipy.special import gamma, gammaln
from scipy.integrate import nquad
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
import glob

lanbda_bar = 0.205
x0 = 0.01
sigma0= 23/0.389
lanbda_QCD = 0.226 #GeV
Q0 = 0.6 #Gev
rapidez= [0.5, 1.0, 1.5, 2.0, 2.4, 3.0, 3.4]
W=5020 #GeV

csv_files = glob.glob('path/*.csv')

list_n=[] 
list_Pb=[]
for i in csv_files:

    df = pd.read_csv(i)
    n = df['N'].to_numpy()
    PB = df['PN'].to_numpy()

    list_n.append(n) 
    list_Pb.append(PB)


def u(W): return 0.24/(0.13 + 0.32*((W/1000)**0.115))

def y(n,W): return (1/2)*np.log((np.sqrt(((np.cosh(n))**2)+((u(W))**2)) + np.sinh(n))/(np.sqrt(((np.cosh(n))**2)+((u(W))**2)) - np.sinh(n)))

def J(n,W): return np.cosh(n)/(np.sqrt((np.cosh(n))**2 + (u(W))**2))

def mod_p_2(p,k,th): return 0.25*(p**2 + k**2 + 2*p*k*np.cos(th))
def mod_m_2(p,k,th): return 0.25*(p**2 + k**2 - 2*p*k*np.cos(th))

def Q2s_1(n,W): return (Q0**2)*(x0*(W/Q0)*np.exp(y(n,W)))**lanbda_bar
def Q2s_2(n,W): return (208**(1/3))*(Q0**2)*(x0*(W/Q0)*np.exp(-y(n,W)))**lanbda_bar

def phi1(p,k,th,n,W): return ((3*sigma0)/(4*(np.pi**2)))*(mod_p_2(p,k,th)/Q2s_1(n,W))*np.exp(-(mod_p_2(p,k,th)/Q2s_1(n,W)))*(1/0.52 if (12*np.pi)/(27*np.log(Q2s_1(n,W)/lanbda_QCD**2))>0.52 else (27*np.log(Q2s_1(n,W)/lanbda_QCD**2))/(12*np.pi))
def phi2(p,k,th,n,W): return ((3*sigma0)/(4*(np.pi**2)))*(mod_m_2(p,k,th)/Q2s_2(n,W))*np.exp(-(mod_m_2(p,k,th)/Q2s_2(n,W)))*(1/0.52 if (12*np.pi)/(27*np.log(Q2s_2(n,W)/lanbda_QCD**2))>0.52 else (27*np.log(Q2s_2(n,W)/lanbda_QCD**2))/(12*np.pi))
def f(th,k,p,n,W): return (3/8)*(2*np.pi)*J(n,W)*(k/p)*phi1(p,k,th,n,W)*phi2(p,k,th,n,W)*0.3

range_p = [0, np.inf]
range_th = [0, 2*np.pi]
range_n=[[-1/2,1/2],[-1,1],[-1.5,1.5],[-2,2],[-2.4,2.4],[-3,3],[-3.4,3.4]]
def range_k(p,n,W):
    return [0, p]


#----------Here's where i use nquad, i want to change this to monte carlo-----------
list_n_medio=[]
for range_eta,rap in zip(range_n,rapidez):
    result, erro = nquad(f, [range_th,range_k,range_p,range_eta],args=(W,))
    print('|η|<%.1f: OK!'%rap)
    list_n_medio.append(result)
#-----------------------------------------------------------------------------------