Help with NIntegrate with singularity and variable integration limit

I have a function of the type:

$$F(x,z)= \int_{-\infty}^{\infty} dx' \int^{z}_{-\infty} dz' f(z,z',x,x')\frac{\partial g(z',x')}{\partial z'}$$

where the integration order may be interchanged. I have tried to use NIntegrate in a function of x and z, but with no success. Th end goal would be to have a function $$F(z,x)$$ and use the Plot function to generate a plot vs z for a given x i.e. $$F(z,0)$$.

The code I made is shown below:

(*CONSTANTS*)
n = 10^-12/(1.6*10^-19);
re = 2.817940920000000*10^-15;
e = 1;
rho = 10;
gamma = 1800/0.511;
beta = Sqrt[1 - 1/gamma^2];
kappa = Sqrt[((x - xprime)/rho)^2 + 4 (1 + (x - xprime)/rho) Sin[phi/2]^2];

(*PHIPOS DERIVATION*)
v = (3 (1 - beta^2 - beta^2 (x/rho)))/(beta^2 (1 + x/rho));
eta = -((6 ((z - zprime)/(2 rho)))/(beta^2 (1 + x/rho)));
chi = (3 (4 ((z - zprime)/(2 rho))^2 - beta^2 ((x - xprime)/rho)^2))/(4 beta^2 (1 + (x - xprime)/rho));
m = -(v/3) + (chi/3 + v^2/36) omega^(-1/3) + omega^(1/3);
omega = eta^2/16 - (chi*v)/6 + v^3/216 + Sqrt[(eta^2/16 - (chi*v)/6 + v^3/216)^2 - (chi/3 + v^2/36)^3];

phiPos = (Sqrt[2 m] + Sqrt[-2 (m + v) - (2 eta)/Sqrt[2 m]]);

(*POTENTIAL WITH PHIPOS*)
potential = (e*beta^2 (Cos[phi] - 1/(1 + (x - xprime)/rho)))/(2 rho^2 (kappa - beta (1 + (x - xprime)/rho) Sin[phi]));
potentialPos = ((2 rho)/beta^2) potential /. phi -> phiPos;

(*GAUSSIAN*)
sigmaZ = 25*10^-6;
sigmaX = 25*10^-6;
gaussian = 1/(2 Pi) Exp[-(zprime^2/(2 sigmaZ^2))] Exp[-(xprime^2/(2sigmaX^2))];
dgaussian = D[gaussian, zprime];

(*FINAL SOLUTION*)
solution[x_?NumericQ, z_?NumericQ] :=  NIntegrate[dgaussian*potentialPos,{zprime, -\[Infinity],z}, {xprime, -\[Infinity], \[Infinity]}]

Plot[solution[0, z], {z, -10^-4, 10^-4}]