Double integral with non constant limits in quadpy?

I need to integrate a function defined in 2Dims (z and radius r), for which I don't have an expression.

I can just query the function at any position (z,r) and get the returned value.

I have the integration range across z: [-z_range, z_range], which I partition into N_z points as:

z_is = -z_range + (np.arange(N_z) + 0.5) * (2.*z_range/N_z)

For each value in z_is, the range to integrate across r is: [0, r_thresh_at_this_z].

r_thresh_at_this_z is obtained from the value of z called z_i as:

def get_r_thresh(z_i):
    return expression_of_z_i # returns a positive float.

So the range on the radial integral is dependent on the value on z.

I have the function f as:

def f(r,z):
    return interpolator_for_f(r,z)

I want to use the quadpy package in the most efficient way as it has been created to be able to be used in this way.

I was thinking to use a for loop to loop through the z_is and to perform a gaussian quadrature across [0, r_thres_for_that_z], for each value in z_is.

I could use:

results = np.zeros((N_z))
errs = np.zeros((N_z))
for i in range(N_z):
   def f(r):
      return interpolator_for_f(r,z_is[i])
   results[i], errs[i] = quadpy.quad(f, [0, r_thresh_at_this_z])

But I feel a for loop is not the most efficient way to use quadpy.

Can you tell me what I am missing in doing this integral only with fast numpy arrays, so no for-loops?

[I have read about the shapes of the input x to the function f. In my case d = 1 because it's a line integral across r. n = N_z because I want to perform N_z such line integrals which I will then add up to obtain 1 single scalar, the result of the (whole) double integral.

p = 1000 because say I want 1000 integration points across r, for each value of z.

So I will need to sample the function at N_z * 1000 points.

Function f shall return an array shaped (N_z, 1000)

Is the identification of these parameters helpful in any way?]

Thank you!