your gaussian function seems to have a factor of 2 in, try this instead:
def gaussian(wz, r, I):
"""
wz: Beam width at z=0;
r : Radial coordinates
I : Intensity distribution
"""
Fin = I*np.exp(-(r/wz)**2)
return Fin
Are you sure that 2 supposed to be there?