def pvTvz(self,vT,vz,R,z,gl=True,ngl=_DEFAULTNGL2,nsigma=4.):
        """
        NAME:

           pvTvz

        PURPOSE:

           calculate the marginalized (vT,vz) probability at this location (NOT normalized by the density)

        INPUT:

           vT - tangential velocity (can be Quantity)

           vz - vertical velocity (can be Quantity)

           R - radius (can be Quantity)

           z - height (can be Quantity)

           gl - use Gauss-Legendre integration (True, currently the only option)

           ngl - order of Gauss-Legendre integration
           
           nsigma - sets integration limits to [-1,+1]*nsigma*sigma_R(R) for integration over vR (default: 4)

        OUTPUT:

           p(vT,vz,R,z)

        HISTORY:

           2012-12-22 - Written - Bovy (IAS)           
           2018-01-12 - Added Gauss-Legendre integration prefactor nsigma/2 - Trick (MPA)

        """
        sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        if gl:
            if ngl % 2 == 1:
                raise ValueError("ngl must be even")
            #Use Gauss-Legendre integration for all
            if ngl == _DEFAULTNGL:
                glx, glw= self._glxdef, self._glwdef
                glx12, glw12= self._glxdef12, self._glwdef12
            elif ngl == _DEFAULTNGL2:
                glx, glw= self._glxdef2, self._glwdef2
                glx12, glw12= self._glxdef, self._glwdef
            else:
                glx, glw= numpy.polynomial.legendre.leggauss(ngl)
                glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2)
            #Evaluate everywhere
            if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                    actionAngle.actionAngleAdiabaticGrid)):
                vRgl= nsigma*sigmaR1/2.*(glx+1.)
                vRglw= glw
                vRfac= nsigma*sigmaR1  #2 x integration over [0,nsigma*sigmaR1]
            else:
                vRgl= nsigma*sigmaR1/2.*(glx12+1.)
                vRgl= list(vRgl)
                vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.))
                vRgl= numpy.array(vRgl)
                vRglw= glw12
                vRglw= list(vRglw)
                vRglw.extend(glw12)
                vRglw= numpy.array(vRglw)
                vRfac = 0.5*nsigma*sigmaR1   #integration over [-nsigma*sigmaR1,0] and [0,nsigma*sigmaR1]
            #evaluate
            logqeval= self(R+numpy.zeros(ngl),
                           vRgl,
                           vT+numpy.zeros(ngl),
                           z+numpy.zeros(ngl),
                           vz+numpy.zeros(ngl),
                           log=True,use_physical=False)
            return numpy.sum(numpy.exp(logqeval)*vRglw*vRfac)