diff --git a/lenspyx/__init__.py b/lenspyx/__init__.py
index 5eebc5d..cff7c8d 100644
--- a/lenspyx/__init__.py
+++ b/lenspyx/__init__.py
@@ -1,3 +1,3 @@
-from lenspyx.lensing import alm2lenmap, alm2lenmap_spin, synfast, get_geom
+from lenspyx.lensing import alm2lenmap, alm2lenmap_spin, synfast, get_geom, dlm2angles
 from._version import __version__
 
diff --git a/lenspyx/cachers.py b/lenspyx/cachers.py
index 3a1a354..d69620d 100644
--- a/lenspyx/cachers.py
+++ b/lenspyx/cachers.py
@@ -24,7 +24,7 @@ def remove(self, fn):
 class cacher_npy(cacher):
     def __init__(self, lib_dir, verbose=False):
         if not os.path.exists(lib_dir):
-            os.makedirs(lib_dir)
+            os.makedirs(lib_dir, exist_ok=True)
         self.lib_dir = lib_dir
         self.verbose = verbose
 
diff --git a/lenspyx/experimental.py b/lenspyx/experimental.py
new file mode 100644
index 0000000..f931849
--- /dev/null
+++ b/lenspyx/experimental.py
@@ -0,0 +1,90 @@
+import numpy as np
+from ducc0.sht import synthesis_general as syngducc, adjoint_synthesis_general as adjsyngducc
+try:
+    import capsht
+except ImportError:
+    print("capsht not found, you will not be able to use the functions in this module")
+try:
+    from capsht.experimental import synthesis_general_cap, synthesis_general_band, adjoint_synthesis_general_cap, adjoint_synthesis_general_band
+except ImportError:
+    print("synthesis_general_cap or synthesis_general_band not found in capsht.experimental, are you up to date?")
+
+def _epsapo(thtcap, epsilon, lmax, version=1, dl_7=None):
+    #dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2
+    assert version == 1, 'C++ code now only implemented for version 1'
+    if version == 0:
+        if dl_7 is None:
+            dl_7 = 15
+        dl = dl_7 * ((-np.log10(epsilon) / (7 )) ** 1) ** 0.5
+    elif version == 1: 
+        if dl_7 is None:
+            dl_7 = 2*7*np.log(10.)/np.pi
+        dl = dl_7 * (-np.log10(epsilon) / (7. ))
+    else:
+        raise ValueError('version %s not implemented'%version)
+    return np.sqrt(dl / lmax * np.pi / thtcap)
+
+def synthesis_general(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, 
+                      thtcap:float=None, eps_apo:float=None, tht_min:float=None, tht_max:float=None, verbose:bool=False, **kwargs):
+    """Wrapper to capsht synthesis_general function, hiding the choice of eps_apo and SHT algorithm
+    
+
+        See ducc0.sht.synthesis_general for arguments, optional arguments and outputs
+
+        relevant keyword, *mode*, *map*, *mmax*
+    
+    """
+    if tht_min is not None and tht_max is not None: # attempt at synthesis_general_band
+        eps_apo = eps_apo or 1.2 * _epsapo(tht_max-tht_min, epsilon, lmax)
+        thta_p = tht_min - 0.5 * eps_apo * (tht_max - tht_min)
+        thtb_p = tht_max + 0.5 * eps_apo * (tht_max - tht_min)
+        if (thta_p >= 0.) and (thtb_p <= np.pi):
+            if verbose:
+                print('syng type: band %.1f deg %.1f deg' % (thta_p/np.pi*180, thtb_p/np.pi*180))
+            assert 0, 'fix band to new scheme'
+            return synthesis_general_band(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, 
+                                      thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs)
+        if thta_p < 0.: # Can try synthesis_general_cap later on
+            thtcap = tht_max
+            eps_apo = None
+    if thtcap is not None: # attempt at synthesis_general_cap
+        eps_apo = eps_apo or _epsapo(thtcap, epsilon, lmax)
+        epsilon_nufft = kwargs.pop('epsilon_nufft', epsilon)
+        if verbose:
+            print('syng type: sent to cap %.1f epsapo %.2f' % (thtcap/np.pi*180, eps_apo))
+        return synthesis_general_cap(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon_nufft, thtcap=thtcap, eps_apo=eps_apo, **kwargs)
+    if verbose:
+        print('syng type : general')
+    return syngducc(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon,  **kwargs)
+
+def adjoint_synthesis_general(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, 
+                      thtcap:float=None, eps_apo:float=None, tht_min:float=None, tht_max:float=None, verbose:bool=False, **kwargs):
+    """Wrapper to capsht synthesis_general function, hiding the choice of eps_apo
+    
+
+        See ducc0.sht.synthesis_general for arguments, optional arguments and outputs
+
+        relevant keyword, *mode*, *alm*, *mmax*
+
+    
+    """
+    if tht_min is not None and tht_max is not None: # attempt at synthesis_general_band
+        eps_apo = eps_apo or 1.2 * _epsapo(tht_max-tht_min, epsilon, lmax)
+        thta_p = tht_min - 0.5 * eps_apo * (tht_max - tht_min)
+        thtb_p = tht_max + 0.5 * eps_apo * (tht_max - tht_min)
+        if (thta_p >= 0.) and (thtb_p <= np.pi):
+            if verbose:
+                print('adjsyng type: band %.1f deg %.1f deg' % (thta_p/np.pi*180, thtb_p/np.pi*180))
+            return adjoint_synthesis_general_band(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, 
+                                      thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs)
+        if thta_p < 0.: # Can try synthesis_general_cap later on
+            thtcap = tht_max
+            eps_apo = None
+    if thtcap is not None: # attempt at synthesis_general_cap
+        eps_apo = eps_apo or _epsapo(thtcap, epsilon, lmax)
+        if verbose:
+            print('adjsyng type: sent to cap %.1f deg, eps_apo %.2f' % (thtcap/np.pi*180, eps_apo))
+        return adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs)
+    if verbose:
+        print('adjsyng type : general')
+    return syngducc(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon,  **kwargs)
\ No newline at end of file
diff --git a/lenspyx/qest/qest.py b/lenspyx/qest/qest.py
index 70ac64d..760b754 100644
--- a/lenspyx/qest/qest.py
+++ b/lenspyx/qest/qest.py
@@ -14,7 +14,7 @@
 from lenspyx.qest.ivfs import OpFilt
 
 
-def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, verbose=False, get_alm2=None, geometry: Geom or None=None):
+def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, mmax=None, mmax_qlm=None, verbose=False, get_alm2=None, geometry: Geom or None=None):
     """Evaluates a quadratic estimator gradient and curl terms.
 
 
@@ -32,13 +32,14 @@ def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, verbose=False, get_
             glm and clm healpy arrays (gradient and curl terms of the QE estimate)
 
     """
+    assert mmax in [None, lmax_ivf], 'there is a bug with non-trivial mmmax, fix this first'
     qe_list = _get_qes(qe_key, lmax_ivf, cls_weight)
     if geometry is None:
         qe_spin = np.max([qe[0].spin_ou + qe[1].spin_ou for qe in uqe.qe_compress(qe_list)])
         geometry = Geom.get_thingauss_geometry((2 * lmax_ivf + lmax_qlm) // 2 + 1, qe_spin)
-    return _eval_qe(qe_list, get_alm, lmax_qlm, verbose=verbose, get_alm2=get_alm2, geo=geometry)
+    return _eval_qe(qe_list, get_alm, lmax_qlm, mmax_qlm=mmax_qlm, mmax=mmax, verbose=verbose, get_alm2=get_alm2, geo=geometry)
 
-def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, get_alm2=None, mmax_qlm:int or None=None, nthreads=0):
+def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, get_alm2=None, mmax:int or None=None, mmax_qlm:int or None=None, nthreads=0):
     """Evaluation of a QE from its list of leg definitions.
 
         Args:
@@ -105,25 +106,25 @@ def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, ge
             if len(conjugate) > 0:
                 for j in conjugate:
                     print("in-spins conjugate leg and out-spin", qes[j][1].spins_in, qes[j][1].spin_ou)
-        a = q[0](get_alm, geo)
+        a = q[0](get_alm, geo, mmax=mmax)
         if qe_spin:
-            dc += fac1 * a * q[1](get_alm2, geo)
+            dc += fac1 * a * q[1](get_alm2, geo, mmax=mmax)
             for j in conjugate:
-                dc += a.conj() * qes[j][1](get_alm2, geo)
+                dc += a.conj() * qes[j][1](get_alm2, geo, mmax=mmax)
         else: # We must consider the real part only
-            dc += fac1 * (a * q[1](get_alm2, geo)).real
+            dc += fac1 * (a * q[1](get_alm2, geo, mmax=mmax)).real
             for j in conjugate:
-                dc += (a.conj() * qes[j][1](get_alm2, geo)).real
+                dc += (a.conj() * qes[j][1](get_alm2, geo, mmax=mmax)).real
         if symmetrize: # same, swapping alm2 and alm1
-            a = q[0](get_alm2, geo)
+            a = q[0](get_alm2, geo, mmax=mmax)
             if qe_spin:
-                dc += fac1 * a * q[1](get_alm, geo)
+                dc += fac1 * a * q[1](get_alm, geo, mmax=mmax)
                 for j in conjugate:
-                    dc += a.conj() * qes[j][1](get_alm, geo)
+                    dc += a.conj() * qes[j][1](get_alm, geo, mmax=mmax)
             else:
-                dc += fac1 * (a * q[1](get_alm, geo)).real
+                dc += fac1 * (a * q[1](get_alm, geo, mmax=mmax)).real
                 for j in conjugate:
-                    dc += (a.conj() * qes[j][1](get_alm, geo)).real
+                    dc += (a.conj() * qes[j][1](get_alm, geo, mmax=mmax)).real
 
     gclm = geo.adjoint_synthesis(m=dr, spin=qe_spin, lmax=lmax_qlm, mmax=mmax_qlm, nthreads=nthreads)
     if symmetrize:
diff --git a/lenspyx/qest/utils_qe.py b/lenspyx/qest/utils_qe.py
index d51d722..ed5c79b 100644
--- a/lenspyx/qest/utils_qe.py
+++ b/lenspyx/qest/utils_qe.py
@@ -22,7 +22,7 @@ def __mul__(self, other):
     def __add__(self, other):
         assert self.spin_in == other.spin_in and self.spin_ou == other.spin_ou
         lmax = max(self.get_lmax(), other.get_lmax())
-        cl = np.zeros(lmax + 1, dtype=float)
+        cl = np.zeros(lmax + 1, dtype=self.cl.dtype)
         cl[:len(self.cl)] += self.cl
         cl[:len(other.cl)] += other.cl
         return qeleg(self.spin_in, self.spin_ou, cl)
@@ -50,7 +50,7 @@ def __iadd__(self, other_qe: qeleg):
         self.cls.append(np.copy(other_qe.cl))
         return self
 
-    def __call__(self, get_alm: callable, geometry: Geom, nthreads: int = 0):
+    def __call__(self, get_alm: callable, geometry: Geom, mmax=None, nthreads: int = 0):
         """Returns the spin-weighted real-space map of the estimator.
 
             We first build X_lm in the wanted _{si}X_lm _{so}Y_lm and then convert this alm2map_spin conventions.
@@ -59,10 +59,10 @@ def __call__(self, get_alm: callable, geometry: Geom, nthreads: int = 0):
         if nthreads <= 0:
             nthreads = cpu_count()
         lmax = self.get_lmax()
-        mmax = lmax
+        if mmax is None:
+            mmax = lmax
         ncomp, npix = 1 + (self.spin_ou != 0), geometry.npix()
-        alm_size = Alm.getsize(lmax, mmax)
-        gclm = np.zeros((ncomp, alm_size), dtype=complex)
+        gclm = np.zeros((ncomp, Alm.getsize(lmax, mmax)), dtype=complex)
         for i, (si, cl) in enumerate(zip(self.spins_in, self.cls)):
             assert si in [0, -2, 2], str(si) + ' input spin not implemented'
             alms = [get_alm('e'), get_alm('b')] if abs(si) == 2 else [-get_alm('t'), 0]
diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py
index c34de23..9544cb0 100644
--- a/lenspyx/remapping/utils_geom.py
+++ b/lenspyx/remapping/utils_geom.py
@@ -4,7 +4,7 @@
 import ducc0
 from ducc0.misc import GL_thetas, GL_weights
 from ducc0.fft import good_size
-from ducc0.sht.experimental import synthesis, adjoint_synthesis, synthesis_deriv1
+from ducc0.sht import synthesis, adjoint_synthesis, synthesis_deriv1, synthesis_general, adjoint_synthesis_general
 
 def st2mmax(spin, tht, lmax):
     r"""Converts spin, tht and lmax to a maximum effective m, according to libsharp paper polar optimization formula Eqs. 7-8
@@ -49,6 +49,14 @@ def npix(self):
         """
         return int(np.sum(self.nph))
 
+    def mmax(self, spin:int, lmax:int):
+        """Safe value of mmax for a given lmax and spin, given the latitude range of the geometry
+         
+         
+        """
+        return min(int(st2mmax(spin, np.pi * 0.5 - np.min(np.abs(np.pi * 0.5 - self.theta)), lmax)) + 1, lmax)
+
+
     def fsky(self):
         """Fractional area of the sky covered by the pixelization
 
@@ -175,6 +183,7 @@ def adjoint_synthesis(self, m: np.ndarray, spin:int, lmax:int, mmax:int, nthread
                                  nthreads=nthreads, ringstart=self.ofs, alm=alm,  **kwargs)
 
     def alm2map_spin(self, gclm:np.ndarray, spin:int, lmax:int, mmax:int, nthreads:int, zbounds=(-1., 1.), **kwargs):
+        # FIXME: method only here for backwards compatiblity
         # FIXME: method only here for backwards compatiblity
         assert zbounds[0] == -1 and zbounds[1] == 1., zbounds
         return self.synthesis(gclm, spin, lmax, mmax, nthreads, **kwargs)
@@ -342,7 +351,7 @@ def get_tgl_geometry(lmax:int, smax:int, good_size_real=True):
 
         """
         return Geom.get_thingauss_geometry(lmax, smax, good_size_real=good_size_real)
-
+    
 
 class pbounds:
     """Class to regroup simple functions handling sky maps longitude truncation
diff --git a/lenspyx/tests/helper.py b/lenspyx/tests/helper.py
index 5eb80d0..6b38751 100644
--- a/lenspyx/tests/helper.py
+++ b/lenspyx/tests/helper.py
@@ -21,14 +21,19 @@ def _extend_cl(cl:np.ndarray, lmax):
         ret[:lmax_cl+1] = cl
         return ret
 
-def syn_alms(spin, lmax_unl=5120, ctyp=np.complex128):
+def syn_alms(spin, lmax_unl=5120, ctyp=np.complex128, white=False):
     ncomp = 1 + (abs(spin) > 0)
     mmax_unl = lmax_unl
     rtyp = lenspyx.remapping.deflection_028.rtype[ctyp]
     eblm = np.empty( (ncomp, Alm.getsize(lmax_unl, mmax_unl)), dtype=ctyp)
-    eblm[0] = synalm(_extend_cl(cls_unl['ee' if abs(spin) > 0 else 'tt'][:lmax_unl + 1], lmax_unl), lmax_unl, mmax_unl, rlm_dtype=rtyp)
+    clgg = _extend_cl(cls_unl['ee' if abs(spin) > 0 else 'tt'][:lmax_unl + 1], lmax_unl)
+    clcc = _extend_cl(cls_unl['bb'][:lmax_unl + 1], lmax_unl)
+    if white:
+        clgg = np.ones_like(clgg)
+        clcc = np.ones_like(clcc)
+    eblm[0] = synalm(clgg, lmax_unl, mmax_unl, rlm_dtype=rtyp)
     if ncomp > 1:
-       eblm[1] = synalm(_extend_cl(cls_unl['bb'][:lmax_unl + 1], lmax_unl), lmax_unl, mmax_unl, rlm_dtype=rtyp)
+       eblm[1] = synalm(clcc, lmax_unl, mmax_unl, rlm_dtype=rtyp)
     return eblm
 
 def syn_dlm(lmax_unl=5120, ctyp=np.complex128):
diff --git a/lenspyx/tests/test_adjsyng_cap.py b/lenspyx/tests/test_adjsyng_cap.py
new file mode 100644
index 0000000..c4eb278
--- /dev/null
+++ b/lenspyx/tests/test_adjsyng_cap.py
@@ -0,0 +1,64 @@
+import numpy as np
+from lenspyx.tests.helper import syn_ffi_ducc_29
+import pylab as pl
+from duccjc.sht import adjoint_synthesis_general_ringweight as adj_syng_w
+from ducc0.sht import adjoint_synthesis_general as adj_syng
+from lenspyx.utils_hp import alm_copy, alm2cl
+from lenspyx.remapping.utils_geom import st2mmax
+from lenspyx.utils_cap import fskycap, eps_opti, args_default
+
+if __name__ == '__main__':
+    """This tests adjoint_synthesis_general_cap
+    
+        This generates delfected positions according to LCDM,and compare output alm of full-sky vs cap routines
+    
+    """
+    args = args_default()
+    args.whiten = True
+
+    thta = 0 / 180 * np.pi
+    thtb = 50 / 180 * np.pi # cap size
+    #dtheta = 7. / 180 * np.pi
+
+
+    ffi, geom = syn_ffi_ducc_29(lmax_len=args.lmax_len, dlmax=args.dlmax, dlmax_gl=args.dlmax_gl, nthreads=args.nt,
+                             verbosity=1, epsilon=10 ** (-args.epsilon))
+    lmax_unl, mmax_unl = args.lmax_len + args.dlmax, args.lmax_len + args.dlmax
+
+    sht_mode = 'STANDARD' if (args.spin == 0 or not args.gonly) else 'GRAD ONLY'
+    ncomp = 1 + (args.spin > 0) * (not args.gonly)
+    ptg = ffi._get_ptg() # locations
+    assert ptg.shape[-1] == 2, ptg.shape
+
+
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=False, northsouth_sym=False)
+    slic = slice(np.min(geom_trc.ofs), np.max(geom_trc.ofs + geom_trc.nph))
+    ptg_sliced = ptg[slic, :]
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=True, northsouth_sym=False)
+    m = np.random.standard_normal((ncomp, ptg_sliced.shape[0]))
+    # Full-sky adjoint synthesis_general, at higher accuracy
+    eblm = adj_syng(map=m, lmax=lmax_unl, mmax=mmax_unl, loc=ptg_sliced, spin=args.spin, epsilon=ffi.epsilon * 0.1,
+                  nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity * 1)
+
+    if thta == 0:    # capped synthesis_general
+        thtcap = np.max(ptg_sliced[:, 0]) * 1.0001
+        eps_apo = eps_opti(lmax_unl, thtcap, dl=args.dl)
+        new_mmax = min(int(st2mmax(args.spin, thtcap * (1. + eps_apo), lmax_unl)) + 1, mmax_unl) if args.adapt_mmax else mmax_unl
+        print("Testing capped syng %2.f deg, fsky %.2f"%( (thtcap / np.pi * 180), fskycap(thtcap)))
+        print('tht cap in deg %.0f, eps apo %.2f'%(thtcap/np.pi * 180, eps_apo))
+        print('induced dlmax %s'%(int(lmax_unl *eps_apo)))
+        print('reduction in mmax by a factor %.1f'%(mmax_unl/new_mmax))
+
+        eblm_2= adj_syng_w(map=m, lmax=lmax_unl, mmax=new_mmax, loc=ptg_sliced,
+                spin=args.spin, epsilon=ffi.epsilon,
+                nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity, thtcap=thtcap, eps_apo=eps_apo, apofct=args.apofct)
+        # Resize alms to full-sky adjoint synthesis general
+        eblm_2 = np.array([alm_copy(alm, new_mmax, lmax_unl, mmax_unl) for alm in eblm_2])
+        for i, (ref, diff) in enumerate(zip(eblm, eblm_2 - eblm)):
+            cldiff = alm2cl(diff, diff, lmax_unl, mmax_unl, lmax_unl)
+            clref = alm2cl(ref, ref, lmax_unl, mmax_unl, lmax_unl)
+            pl.semilogy(np.sqrt(cldiff[2:]/clref[2:]), label='comp %s capped, weighted (dl=%s)'%(i, args.dl))
+            pl.xlabel(r'$\ell$')
+            pl.ylabel(r'$C_\ell$')
+    pl.legend()
+    pl.show()
\ No newline at end of file
diff --git a/lenspyx/tests/test_syng_cap.py b/lenspyx/tests/test_syng_cap.py
new file mode 100644
index 0000000..231c1a5
--- /dev/null
+++ b/lenspyx/tests/test_syng_cap.py
@@ -0,0 +1,63 @@
+import numpy as np
+from lenspyx.tests.helper import syn_ffi_ducc_29, syn_alms
+
+import pylab as pl
+from duccjc.sht import synthesis_general_ringweight as syng_w, adjoint_synthesis_general_ringweight as adj_syng_w
+from ducc0.sht import synthesis_general as syng, adjoint_synthesis_general as adj_syng
+from lenspyx.utils_hp import alm_copy
+from lenspyx.remapping import deflection_029, deflection_028 as deflection_28
+from lenspyx.remapping.utils_geom import st2mmax
+from lenspyx.utils_cap import fskycap, eps_opti, examine, args_default
+
+if __name__ == '__main__':
+    """This tests synthesis_general_cap
+
+        This generates delfected positions according to LCDM,and compare output maps of full-sky vs cap routines
+
+    """
+    args = args_default()
+    thta = 0 / 180 * np.pi
+    thtb = 50 / 180 * np.pi # cap size
+#    args.whiten = False
+
+
+    ffi, geom = syn_ffi_ducc_29(lmax_len=args.lmax_len, dlmax=args.dlmax, dlmax_gl=args.dlmax_gl, nthreads=args.nt,
+                             verbosity=1, epsilon=10 ** (-args.epsilon))
+    lmax_unl, mmax_unl = args.lmax_len + args.dlmax, args.lmax_len + args.dlmax
+
+    eblm = syn_alms(args.spin, lmax_unl=lmax_unl, ctyp=np.complex64 if ffi.single_prec else np.complex128, white=args.whiten)
+    eblm = np.atleast_2d(eblm)
+    if args.gonly:
+        eblm = eblm[:1]
+    sht_mode = deflection_28.ducc_sht_mode(eblm, args.spin)
+    ptg = ffi._get_ptg() # locations
+    assert ptg.shape[-1] == 2, ptg.shape
+
+
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=False, northsouth_sym=False)
+    slic = slice(np.min(geom_trc.ofs), np.max(geom_trc.ofs + geom_trc.nph))
+    ptg_sliced = ptg[slic, :]
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=True, northsouth_sym=False)
+
+    # Full-sky synthesis_general, at higher accuracy
+    values = syng(lmax=lmax_unl, mmax=mmax_unl, alm=eblm, loc=ptg, spin=args.spin, epsilon=ffi.epsilon * 0.1,
+                  nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity * 1)
+
+    if thta == 0:    # capped synthesis_general
+        thtcap = np.max(ptg_sliced[:, 0]) * 1.0001
+        eps_apo = eps_opti(lmax_unl, thtcap, dl=args.dl)
+        new_mmax = min(int(st2mmax(args.spin, thtcap * (1. + eps_apo), lmax_unl)) + 1, mmax_unl) if args.adapt_mmax else mmax_unl
+        print("Testing capped syng %2.f deg, fsky %.2f"%( (thtcap / np.pi * 180), fskycap(thtcap)))
+        print('tht cap in deg %.0f, eps apo %.2f'%(thtcap/np.pi * 180, eps_apo))
+        print('induced dlmax %s'%(int(lmax_unl *eps_apo)))
+        print('reduction in mmax by a factor %.1f'%(mmax_unl/new_mmax))
+
+        gclm = np.array([alm_copy(alm, mmax_unl, lmax_unl, new_mmax) for alm in eblm])
+        values_w = syng_w(lmax=lmax_unl, mmax=new_mmax, alm=gclm, loc=ptg_sliced,
+                ringweights=np.array([1.]), spin=args.spin, epsilon=ffi.epsilon,
+                nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity, thtcap=thtcap, eps_apo=eps_apo, apofct=args.apofct)
+        print(np.max(np.abs(values_w - values[:, slic])))
+        for i, (val, diff) in enumerate(zip(values[:, slic], values_w-values[:, slic])):
+            examine(val, diff, geom_trc, ffi.epsilon, label='comp %s capped, weighted (dl=%s)'%(i, args.dl))
+    pl.legend()
+    pl.show()
\ No newline at end of file
diff --git a/lenspyx/utils_cap.py b/lenspyx/utils_cap.py
new file mode 100644
index 0000000..25b4a74
--- /dev/null
+++ b/lenspyx/utils_cap.py
@@ -0,0 +1,126 @@
+import numpy as np
+from lenspyx.remapping.utils_geom import Geom
+import pylab as pl
+
+class args_default:
+    def __init__(self):
+        self.lmax_len = 4000
+        self.dlmax_gl = 500
+        self.spin = 2
+        self.nt = 4
+        self.HL = 0
+        self.alloc = 0
+        self.tracemalloc = False
+        self.epsilon = 7
+        self.gonly = False
+        self.dlmax = 500
+        self.whiten = True # make the spectra white before interpolation
+        self.apofct = 0
+        self.dl=7
+        self.adapt_mmax = True
+
+
+def fskycap(thtcap):
+    return (1. - np.cos(thtcap)) * 0.5
+
+
+class transition01:
+    """Helpers for transition functions equal to 1 at 0 and 0 at 1"""
+    @staticmethod
+    def _eval01(x, version):
+        assert np.all( (x < 1) & (x > 0))
+        if version in ['ES', 'es', 3]:
+            # Exponential of semi-circle, with b giving 1/2 at 1/2
+            beta = np.log(0.5) / (np.sqrt(3/4.) - 1.)
+            return np.exp(beta * (np.sqrt(1. - x ** 2) - 1))
+        if version in ['KB', 'kb', 4]: # Kaiser-bessel
+            b = 5.74
+            return np.i0(b * np.sqrt(1. - x ** 2)) / np.i0(b)
+        if version in [0]: # Smooth typical bump function. this is 1/2 at 1/2
+            fx  = np.exp(-1. / x)
+            f1x = np.exp(-1. / (1 - x))
+            return f1x / (fx + f1x)
+        if version in ['Hann', 1]:
+            return np.cos(x * np.pi * 0.5) ** 2
+        if version in [2]: # another smooth bump
+            return np.exp(-1. /( 1 - x ** 2) + 1.)
+
+    @staticmethod
+    def eval(x, version):
+        ret = np.zeros(x.size)
+        ret[np.where(x <= 0.)] = 1.
+        ret[np.where(x >= 1.)] = 0.
+        i = np.where( (x > 0) & (x < 1))
+        ret[i] = transition01._eval01(x[i], version)
+        return ret
+
+    @staticmethod
+    def eval_fft(npts, version):
+        x = np.arange(npts) * ( (2 * np.pi) / npts ) - np.pi
+        return np.fft.fft(transition01.eval(x, version))
+
+
+    @staticmethod
+    def bump(x, x1, dx, version):
+        """bump, equal to 1 on [-1, 1], and 0 outside of [-1-dx, 1+dx]
+
+        """
+        assert dx >= 0, dx
+        ret = np.zeros(x.size)
+        ax = np.abs(x)
+        ret[np.where(ax <= x1)] = 1.
+        ret[np.where(ax >= (x1 + dx))] = 0.
+        i = np.where( (ax > x1) & (ax < (x1 + dx)) )
+        ret[i] = transition01.eval( (ax[i] - x1) / dx, version)
+        return ret
+
+    @staticmethod
+    def bump_fft(npts, x1, dx, version):
+        x = np.arange(npts) * ( (2 * np.pi) / npts ) - np.pi
+        bp = transition01.bump(x, x1, dx, version)
+        return np.fft.fft(bp)
+
+def Nc(lmax, dl, thetacap, dtheta):
+    """
+
+    Args:
+        lmax: band-limit of the alm array
+        dl: Effective band-limit of the apodization window (for white spectra, dl=7 gives single precision results?)
+        thetacap: co-latitude angle defining the cap
+        dtheta: angular distance dedicated to apodization
+
+        This should have a mimimum at sqrt(dl / lmax * (thetcap /pi))
+
+    Returns:
+
+    """
+    return (lmax + np.pi * dl / dtheta) * (thetacap + dtheta) / np.pi
+
+# parametrize interval
+def eps_opti(lmax, thetacap, dl=7):
+    """Guess of optimum choice of apodization length, given the alm band-limit lmax and coordinate of the cap, and the [0, pi) window bandlimit"""
+    guess = np.sqrt(dl / lmax * np.pi / thetacap * np.pi)
+    return min(guess, (np.pi / thetacap - 1.)*(1.-1e-13)) # Cant overshoot
+
+def examine(ref, diff, geom:Geom, epsilon_ref, label=''):
+    """Look at difference to the ref map ring by ring
+
+
+    """
+    rad2deg = 180 / np.pi
+    # extracts rings and plot rms dev
+    rms = np.zeros(geom.theta.size)
+    thta = np.min(geom.theta)
+    thtb = np.max(geom.theta)
+    dtheta = 0.1 * (thtb - thta)
+    for i, ir in enumerate(np.argsort(geom.theta)):
+        pix = geom.rings2pix(geom, [ir])
+        rms[i] = np.sqrt(np.mean(diff[pix] ** 2)) /np.sqrt(np.mean(ref[pix] ** 2))
+    pl.semilogy(np.sort(geom.theta) * rad2deg, rms, label=label)
+    pl.axhline(epsilon_ref,c='k')
+    pl.axvline(thta * rad2deg, c='grey')
+    pl.axvline(thtb * rad2deg, c='grey')
+    pl.xlim( max((thta-1.2 * dtheta) * rad2deg, 0.), min((thtb + 1.2 * dtheta) * rad2deg, np.pi * rad2deg))
+    pl.xlabel(r'$\theta$ [deg]')
+    pl.ylabel(r'rel dev. (rms)')
+    return rms
\ No newline at end of file
diff --git a/lenspyx/utils_hp.py b/lenspyx/utils_hp.py
index d751844..2a77fc0 100644
--- a/lenspyx/utils_hp.py
+++ b/lenspyx/utils_hp.py
@@ -9,7 +9,7 @@ def almxfl(alm:np.ndarray, fl:np.ndarray, mmax:int or None, inplace:bool):
     Parameters
     ----------
     alm : array
-      The alm to multiply
+      The alm, or alms to multiply
     fl : array
       The function (at l=0..fl.size-1) by which alm must be multiplied.
     mmax : None or int
@@ -24,20 +24,20 @@ def almxfl(alm:np.ndarray, fl:np.ndarray, mmax:int or None, inplace:bool):
       if inplace is True.
 
     """
-    lmax = Alm.getlmax(alm.size, mmax)
+    lmax = Alm.getlmax(alm.shape[-1], mmax)
     if mmax is None or mmax < 0:
         mmax = lmax
     assert fl.size > lmax, (fl.size, lmax)
     if inplace:
         for m in range(mmax + 1):
             b = m * (2 * lmax + 1 - m) // 2 + m
-            alm[b:b + lmax - m + 1] *= fl[m:lmax+1]
+            alm[...,b:b + lmax - m + 1] *= fl[m:lmax+1]
         return
     else:
         ret = np.empty_like(alm)
         for m in range(mmax + 1):
             b = m * (2 * lmax + 1 - m) // 2 + m
-            ret[b:b + lmax - m + 1] = alm[b:b + lmax - m + 1] * fl[m:lmax+1]
+            ret[...,b:b + lmax - m + 1] = alm[...,b:b + lmax - m + 1] * fl[m:lmax+1]
         return ret
 
 
@@ -63,7 +63,7 @@ def gauss_beam(fwhm:float, lmax:int):
     return bl
 
 
-def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64):
+def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64, seed=None, rngen=None, alms_r:np.ndarray=None):
     """Creates a Gaussian field alm from input cl array
 
     Parameters
@@ -83,21 +83,31 @@ def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64):
         harmonic coefficients of Gaussian field with lmax, mmax parameters
 
     """
+    if rngen is not None:
+        assert hasattr(rngen, 'standard_normal'), 'rngen must have standard_normal method'
     assert lmax + 1 <= cl.size
     if mmax is None or mmax < 0:
         mmax = lmax
     alm_size = Alm.getsize(lmax, mmax)
-    alm = rng.standard_normal(alm_size, dtype=rlm_dtype) + 1j * rng.standard_normal(alm_size, dtype=rlm_dtype)
-    almxfl(alm, np.sqrt(cl[:lmax+1] * 0.5), mmax, True)
+    rng = default_rng(seed) if rngen is None else rngen
     real_idcs = Alm.getidx(lmax, np.arange(lmax + 1, dtype=int), 0)
-    alm[real_idcs] = alm[real_idcs].real * np.sqrt(2.)
+    if alms_r is None:
+        alm = rng.standard_normal(alm_size, dtype=rlm_dtype) + 1j * rng.standard_normal(alm_size, dtype=rlm_dtype)
+        almxfl(alm, np.sqrt(cl[:lmax+1] * 0.5), mmax, True)
+    else:
+        assert alms_r.shape == (2, alm_size), alms_r.shape
+        rng.standard_normal((2, alm_size), dtype=rlm_dtype, out=alms_r)
+        almxfl(alms_r[0], np.sqrt(cl[:lmax+1] * 0.5), mmax, True)
+        almxfl(alms_r[1], np.sqrt(cl[:lmax+1] * 0.5), mmax, True)
+        alm = alms_r
+    alm[...,real_idcs] = alm[...,real_idcs].real * np.sqrt(2.)
     return alm
 
 
-def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = np.float64):
+def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np.float64, rngen=None):
     """Creates Gaussian field alms from input cl dictionary
 
-    Parametersseed
+    Parameters
     ----------
     cls : dict
         The power spectra of the maps (e.g. as coming from CAMB)
@@ -105,8 +115,11 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n
         Maximum multipole simulated
     mmax: int
         Maximum m defining the alm layout, defaults to lmax if None or < 0
-    rlm_dtype(optional, defaults to np.float64):
+    seed: (optional, defaults to None)
+        Random generator seed
+    rlm_dtype:(optional, defaults to np.float64)
         Precision of real components of the array (e.g. np.float32 for single precision output array)
+    rngen: (optional, defaults to None) Allows to use a custom random generator, seed is ignored in this case.
 
     Returns
     -------
@@ -114,6 +127,8 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n
         harmonic coefficients of Gaussian field with lmax, mmax parameters
 
     """
+    if rngen is not None:
+        assert hasattr(rngen, 'standard_normal'), 'rngen must have standard_normal method'
     lmax_cls = np.max([len(cl) - 1 for cl in cls.values()])
     if lmax is None:
         lmax = lmax_cls
@@ -156,7 +171,7 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n
         m[:] = np.dot(v, np.dot(np.diag(np.sqrt(t)), v.T))
     # Build phases:
     alm_size = Alm.getsize(lmax, mmax)
-    rng = default_rng(seed)
+    rng = default_rng(seed) if rngen is None else rngen
     phases = 1j * rng.standard_normal((ncomp, alm_size), dtype=rlm_dtype)
     phases += rng.standard_normal((ncomp, alm_size), dtype=rlm_dtype)
     phases *= np.sqrt(0.5)
@@ -174,7 +189,7 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n
     elif rlm_dtype == np.float64:
         dtype_complex = np.complex128
     else:
-        assert 0, "please either choose np.float32 (single precission), or np.float64 (double precission) as rlm_dtype"
+        assert 0, "please either choose np.float32 (single precision), or np.float64 (double precision) as rlm_dtype"
     alms = np.zeros((len(labels_wgrad), phases[0].size), dtype=dtype_complex)
     #    for L in Ls: #L @ L.T is full matrx
     for i, f in enumerate(labels):
diff --git a/setup.py b/setup.py
index ef75415..62c17c8 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@
 setup(
     name='lenspyx',
     version=__version__,
-    packages=['lenspyx', 'lenspyx.remapping', 'lenspyx.tests', 'lenspyx.wigners'],
+    packages=['lenspyx', 'lenspyx.remapping', 'lenspyx.tests', 'lenspyx.wigners', 'lenspyx.qest'],
     url='https://github.com/carronj/lenspyx',
     author='Julien Carron',
     data_files=[('lenspyx/data/cls', ['lenspyx/data/cls/FFP10_wdipole_lensedCls.dat',