diff --git a/README.md b/README.md
index 814d783..e473049 100644
--- a/README.md
+++ b/README.md
@@ -38,4 +38,4 @@ The notebook 'demo_lensit.ipynb' shows an example of iterative lensing map recon
 Other example and tests scripts might follow, or you may just write to me.
 
 ![alt text](https://erc.europa.eu/sites/default/files/content/erc_banner-vertical.jpg)
-![SNSF logo](./docs/SNF_logo_standard_web_color_neg_e.svg)
\ No newline at end of file
+![SNSF logo](./docs/SNF_logo_standard_web_color_neg_e.svg)
diff --git a/lensit/__init__.py b/lensit/__init__.py
index b1943c7..d173335 100644
--- a/lensit/__init__.py
+++ b/lensit/__init__.py
@@ -358,6 +358,21 @@ def get_config(exp):
         Beam_FWHM_amin = 1.
         ellmin=10
         ellmax=3000
+    elif exp == 'S4_sayan':
+        sN_uKamin = 1.
+        Beam_FWHM_amin = 1.
+        ellmin=100
+        ellmax=4000
+    elif exp == 'S4_no_noise':
+        sN_uKamin = 0.
+        Beam_FWHM_amin = 1.
+        ellmin=100
+        ellmax=4000
+    elif exp == 'CMB_HD':
+        sN_uKamin = 0.5
+        Beam_FWHM_amin = 0.5
+        ellmin=100
+        ellmax=4000
     elif exp == 'Planck45':
         sN_uKamin = 45.
         Beam_FWHM_amin = 5
diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py
index b3283ce..e3d2d70 100644
--- a/lensit/clusterlens/lensingmap.py
+++ b/lensit/clusterlens/lensingmap.py
@@ -1,42 +1,61 @@
+from __future__ import print_function
+
 from lensit.clusterlens.profile import  profile
 import numpy as np
 from lensit import ell_mat, pbs
 from lensit.sims import ffs_phas, ffs_cmbs, ffs_maps
-from lensit.misc.misc_utils import enumerate_progress, gauss_beam
+from lensit.misc.misc_utils import enumerate_progress, gauss_beam, pp_to_kk, p_to_k
 from lensit.ffs_deflect import ffs_deflect
 import lensit as li
 from os.path import join as opj
 import os
 from camb import CAMBdata
 
+import os
+import pickle as pk
 
-def get_cluster_libdir(cambinifile, profilename, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp):
-    return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename, 
+import numpy as np
+
+from lensit.sims.ffs_cmbs import sim_cmb_unl, sims_cmb_len
+from lensit.pbs import pbs
+from lensit.misc.misc_utils import npy_hash
+from lensit.sims import ffs_phas, sims_generic
+from lensit.ffs_deflect import ffs_deflect
+
+def get_cluster_libdir(cambinifile, profilename, key, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp):
+    return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename,"lensed_by_%s" % key, 
     'npix%s_lpix_%samin_lmaxsky%s' % (npix, lpix_amin, ellmax_sky), 'M200_%.6E_z%s' % (M200, z), '%s_sims' % nsims)
 
 
 
 class cluster_maps(object):
-    def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, cmb_exp='5muKamin_1amin', cache_maps=False):
+    def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, ellmin_data=10, cmb_exp='5muKamin_1amin', cache_maps=False):
         """Library for flat-sky CMB simulations lensed by a galaxy cluster.
 
         Args:
             lib_dir: various things will be cached there
+            key: a keyword among "lss", "cluster" and "lss_plus_cluster", which is gonna decide how the patch has been lensed
             npix: number of pixels on one side of the square box
             lpix_amin: physical size (in arcmin) of the pixels
             nsims: number of CMB maps to simulate
             cosmo: Instantiated camb.results.CAMBdata object 
-            profparams: dict containing the parameters defining the profile 
+            profparams: dict containing the parameters defining the profile M200, z, xmaxn
             profilename: string defining the density profile of the cluster (e.g. nfw)
             ellmax_sky: maximum multipole of CMB spectra used to generate the CMB maps 
     """
+        assert key in ['lss', 'cluster', 'lss_plus_cluster'], key
         self.libdir = libdir
+        self.key = key
         self.cosmo = cosmo
         if profilename == 'nfw':
             self.M200 = profparams['M200c']
             self.z = profparams['z']
+            self.xmaxn = profparams['xmaxn']
         self.npix = npix
         self.lpix_amin = lpix_amin
+        lbox_amin = npix*lpix_amin
+        lbox_rad = (lbox_amin/60)*(np.pi/180)
+        self.lbox_rad = lbox_rad
         shape = (self.npix, self.npix)
         lpix_rad = self.lpix_amin * np.pi / 180 / 60
         lsides = (lpix_rad*self.npix, lpix_rad*self.npix)
@@ -50,13 +69,31 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd
         self.num_threads=int(os.environ.get('OMP_NUM_THREADS', 1))
         self.lib_skyalm =  ell_mat.ffs_alm_pyFFTW(self.ellmat, num_threads=self.num_threads, filt_func=lambda ell: ell <= self.ellmax_sky)
 
-        nfields = 3 
+        #nfields = 4 
+        if self.key == 'cluster':
+            nfields = 4
+        else:
+            nfields = 4
 
         camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T
-        self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]}
-        for cl in self.cls_unl.values():
-            cl[6000:] = 0
+        #camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T/100
+
+        #camb_cls = cosmo.get_lensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T
+       
+        cpp_fid = cosmo.get_lens_potential_cls(lmax=self.ellmax_sky, raw_cl=True).T[0]
         
+        #if self.key == 'cluster':
+        #    self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]}
+        #else:
+        #    self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid}
+        self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid}
+        
+        # Set the high ell to zero except for the lensing potential cpp_fid
+        for key, cl in zip(self.cls_unl.keys(), self.cls_unl.values()):
+            if key != 'pp':
+                cl[6000:] = 0
+                
+                
         # Generate the CMB random phases
         skypha_libdir = opj(self.libdir,  'len_alms', 'skypha')
         skypha = ffs_phas.ffs_lib_phas(skypha_libdir, nfields, self.lib_skyalm, nsims_max=nsims)
@@ -66,19 +103,34 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd
         pbs.barrier()
 
         self.haloprofile = profile(self.cosmo, profilename)
+        self.xmax = self.xmaxn * self.haloprofile.get_concentration(self.M200, self.z)
+        self.kappa0 = self.haloprofile.get_kappa0(self.M200, self.z, self.xmax)
 
         lencmbs_libdir = opj(self.libdir, 'len_alms')
 
         # Get the noise and beam of the experiment
         sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(self.cmb_exp)
-        self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha)
-        
+        #self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha)
+        #self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha)
+        #self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha)
+        #self.len_cmbs = sim_cmb_unl(self.cls_unl, lib_pha=skypha)
+        if self.key == 'lss':
+            self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha)
+        elif self.key == 'cluster':
+            self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.xmaxn, self.haloprofile, lib_pha=skypha)
+        elif self.key == 'lss_plus_cluster':
+            self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.xmaxn, self.haloprofile, lib_pha=skypha)
+        else:
+            assert 0
+            
         # Set the lmax of the data maps
         self.ellmax_data = ellmax_data
+        self.ellmin_data = ellmin_data
         self.cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=self.ellmax_data)
 
         # The resolution of the data map could be lower than the resolution of the sky map (=the true map)
-        self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads)
+        #self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads)
+        self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: np.logical_and(ell <= self.ellmax_data , ell>=self.ellmin_data), num_threads=self.num_threads)
         
         vcell_amin2 = np.prod(self.lib_datalm.ell_mat.lsides) / np.prod(self.lib_datalm.ell_mat.shape) * (180 * 60. / np.pi) ** 2
         nTpix = sN_uKamin / np.sqrt(vcell_amin2)
@@ -176,18 +228,184 @@ def get_kappa_map(self, M200, z, xmax=None):
             convergence map: numpy array of shape self.lib_skyalm.shape
         """
         return self.haloprofile.kappa_map(M200, z, self.lib_skyalm.shape, self.lib_skyalm.lsides, xmax)
+    
+    def get_kappa0_from_sim(self, lmin, lmax, phi_obs_lm, NL):
+        """ Get the kappa_0 estimate from one simulation
+        Args: 
+            lmin, lmax: the lrange for the sum in the numerator and denominator of the estimator
+            phi_obs_lm: phi (lensing potential) estimate of a single simulation/observation
+            NL: reconstruction niose of the phi estimate
+        Returns:
+            kappa_0 estiamte for the single phi simulation/observation
+            
+        P.S. the arguments is taken as the phi (not kappa) of the observation and accordingly the noise
+        """
+        ell = np.where(self.lib_skyalm.get_Nell()[:lmax+1] > 1)[0]
+        ell = ell[ell >= lmin]
+        #kappa_l = kappa_ell_lensit(ell)*lbox_rad
+        kappa_l = self.haloprofile.analitic_kappa_ft(self.M200, self.z, self.xmax, ell)
+        num_l = self.lib_skyalm.get_Nell()[ell]
+        phi_obs_el = self.lib_skyalm.bin_realpart_inell(phi_obs_lm) #the average 
+        denom_int = (num_l)*(kappa_l/self.lbox_rad/self.kappa0)**2 / (NL[ell]*pp_to_kk(ell))
+        num_int = (num_l)*(kappa_l/self.lbox_rad/self.kappa0) * phi_obs_el[ell]*p_to_k(ell) / (NL[ell]*pp_to_kk(ell)) ## since the observed quantity is phi and we need kappa
+        denom1 = np.sum(denom_int)
+        numer1 = np.sum(num_int)
+        return numer1/denom1     
+
+verbose = False
 
+def get_fields(cls):
+    fields = ['p', 't', 'e', 'b', 'o']
+    ret = ['p', 't', 'e', 'b', 'o']
+    for f in fields:
+        if not ((f + f) in cls.keys()): ret.remove(f)
+    for k in cls.keys():
+        for f in k:
+            if f not in ret: ret.append(f)
+    return ret
 
 
 class sim_cmb_len_cluster(ffs_cmbs.sims_cmb_len):
-    def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False):
+    def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, xmaxn:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False):
         super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False)
 
-        self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides)
+        self.conc_par = haloprofile.get_concentration(M200,z)
+        self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= xmaxn * self.conc_par)
+        #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides)
+        #self.kappa_map = 1e2*haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par)
         self.defl_map = haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) 
 
-        assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors."
+        #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors."
 
     def _get_f(self, idx=None):
         """We refine the displacement field using the convergence map of the cluster"""
         return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)
+    
+class sim_cmb_len_lss_plus_cluster:
+    def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, xmaxn:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False):
+        #super(sim_cmb_len_lss_plus_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False)
+
+        self.M200 = M200
+        self.z = z
+        self.haloprofile = haloprofile
+        self.conc_par = self.haloprofile.get_concentration(M200,z)
+        self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= xmaxn * self.conc_par)
+        #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides)
+        self.defl_map = self.haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) 
+        
+        if not os.path.exists(lib_dir) and pbs.rank == 0:
+            os.makedirs(lib_dir)
+        pbs.barrier()
+        self.lib_skyalm = lib_skyalm
+        fields = get_fields(cls_unl)
+        if lib_pha is None and pbs.rank == 0:
+            lib_pha = ffs_phas.ffs_lib_phas(os.path.join(lib_dir, 'phas'), len(fields), lib_skyalm)
+        else:  # Check that the lib_alms are compatible :
+            assert lib_pha.lib_alm == lib_skyalm
+        pbs.barrier()
+
+        self.unlcmbs = sim_cmb_unl(cls_unl, lib_pha)
+        self.Pool = use_Pool
+        self.cache_lens = cache_lens
+        fn_hash = os.path.join(lib_dir, 'sim_hash.pk')
+        if not os.path.exists(fn_hash) and pbs.rank == 0:
+            pk.dump(self.hashdict(), open(fn_hash, 'wb'), protocol=2)
+        pbs.barrier()
+        sims_generic.hash_check(self.hashdict(), pk.load(open(fn_hash, 'rb')))
+        self.lib_dir = lib_dir
+        self.fields = fields
+
+        #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors."
+
+    def _get_f_cluster(self, idx=None):
+        """We refine the displacement field using the convergence map of the cluster"""
+        return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)
+    
+    def hashdict(self):
+        return {'unl_cmbs': self.unlcmbs.hashdict(),
+                'M200': self.M200,
+                'z': self.z}#,
+                #'haloprofile': self.haloprofile}
+
+    def is_full(self):
+        return self.unlcmbs.lib_pha.is_full()
+
+    def get_sim_plm(self, idx):
+        return self.unlcmbs.get_sim_plm(idx)
+
+    def get_sim_olm(self, idx):
+        return self.unlcmbs.get_sim_olm(idx)
+
+    def _get_f(self, idx):
+        if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields:
+            plm = self.get_sim_plm(idx)
+            olm = self.get_sim_olm(idx)
+            return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False)
+        elif 'p' in self.unlcmbs.fields:
+            plm = self.get_sim_plm(idx)
+            return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm)
+        elif 'o' in self.unlcmbs.fields:
+            olm = self.get_sim_olm(idx)
+            return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False)
+        else:
+            assert 0
+            
+    def _get_f_lss_plus_cluster(self, idx):
+        if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields:
+            plm = self.get_sim_plm(idx)
+            olm = self.get_sim_olm(idx)
+            return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)
+        elif 'p' in self.unlcmbs.fields:
+            plm = self.get_sim_plm(idx)
+            return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)
+        elif 'o' in self.unlcmbs.fields:
+            olm = self.get_sim_olm(idx)
+            return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)
+        else:
+            assert 0
+
+    def get_sim_alm(self, idx, field):
+        if field == 't':
+            return self.get_sim_tlm(idx)
+        elif field == 'p':
+            return self.get_sim_plm(idx)
+        elif field == 'o':
+            return self.get_sim_olm(idx)
+        elif field == 'q':
+            return self.get_sim_qulm(idx)[0]
+        elif field == 'u':
+            return self.get_sim_qulm(idx)[1]
+        elif field == 'e':
+            return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[0]
+        elif field == 'b':
+            return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[1]
+        else:
+            assert 0, (field, self.fields)
+
+    def get_sim_tlm(self, idx):
+        fname = os.path.join(self.lib_dir, 'sim_%04d_tlm.npy'%idx)
+        if not os.path.exists(fname):
+            Tlm = self._get_f(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool)
+            Tlm = self._get_f_cluster(idx).lens_alm(self.lib_skyalm, Tlm, use_Pool=self.Pool)
+            #Tlm = self._get_f_lss_plus_cluster(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool)
+            if not self.cache_lens: return Tlm
+            np.save(fname, Tlm)
+        return np.load(fname)
+
+    def get_sim_qulm(self, idx):
+        fname = os.path.join(self.lib_dir, 'sim_%04d_qulm.npy'%idx)
+        if not os.path.exists(fname):
+            Qlm, Ulm = self.lib_skyalm.EBlms2QUalms(
+                np.array([self.unlcmbs.get_sim_elm(idx), self.unlcmbs.get_sim_blm(idx)]))
+            f = self._get_f(idx)
+            #f = self._get_f_lss_plus_cluster(idx)
+            Qlm = f.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool)
+            Ulm = f.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool)
+            f_cluster = self._get_f_cluster(idx)
+            Qlm = f_cluster.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool)
+            Ulm = f_cluster.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool)
+            if not self.cache_lens: return np.array([Qlm, Ulm])
+            np.save(fname, np.array([Qlm, Ulm]))
+        return np.load(fname)
+    
+    
\ No newline at end of file
diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py
index 6053b11..cbeb360 100644
--- a/lensit/clusterlens/profile.py
+++ b/lensit/clusterlens/profile.py
@@ -41,7 +41,7 @@ def get_rho_s(self, M200, z, const_c=None):
     def get_rs(self, M200, z, const_c=None):
         return self.get_r200(M200, z) / self.get_concentration(M200, z, const_c)
     
-    def get_thetas_amin(M200,z, const_c=None):
+    def get_thetas_amin(self, M200,z, const_c=None):
         return self.r_to_theta(z, self.get_rs(M200, z, const_c))
 
     def get_r200(self, M200, z):
@@ -70,11 +70,10 @@ def sigma_crit(self, z):
         Dang_LS = self.cosmo.angular_diameter_distance2(z, self.zcmb)
         return const.c_Mpcs**2 /(4*np.pi*const.G_Mpc3_pMsol_ps2) *  Dang_OS/Dang_OL/Dang_LS
     
-    def get_kappa0(self, M200, z):
+    def get_kappa0(self, M200, z, xmax=None):
         """Return the value of the lensing convergence profile for x = 1, i.e. R = rs"""
-        rs = self.get_rs(M200, z)
-        rhos = self.get_rho_s(M200, z)
-        return 2*rhos*rs/3 / self.sigma_crit(z)
+        thsc = self.get_thetas_amin(M200,z)
+        return self.kappa_theta(M200, z, thsc, xmax=xmax)[0]
 
     def sigma_nfw(self, M200, z, R, xmax=None):
         """Analytic expression for the surface mass desinty of a NFW profile 
@@ -101,26 +100,36 @@ def sigma_nfw(self, M200, z, R, xmax=None):
     def fx(self, x, tol=6):
         """This integral of the NFW profile along the line of sight is integrated up to infinity
         See Bartelmann 1996 or Wright et al 1999"""
+        ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7
         f = np.zeros_like(x)
         xp = np.where(x>1)
         xo = np.where(np.abs(x-1.) < 10**(-tol)) 
-        xm = np.where(x<1)
+        xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        #xm = np.where(x<1)
+        xd = np.where(x<self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        idx = self.theta_amin_to_x(2*1e14, 0.7, 0.12)
         f[xp] = (1 - 2/np.sqrt(x[xp]**2 - 1) * np.arctan(np.sqrt((x[xp] - 1)/(x[xp] + 1))))/(x[xp]**2-1)
         f[xm] = (1 - 2/np.sqrt(1 - x[xm]**2) * np.arctanh(np.sqrt((1 - x[xm])/(1 + x[xm]))))/(x[xm]**2-1)
         f[xo] = 1/3
+        f[xd] = (1 - 2/np.sqrt(1 - idx**2) * np.arctanh(np.sqrt((1 - idx)/(1 + idx))))/(idx**2-1)
         return f
 
 
     def gx(self, x, xmax, tol=5):
         """We apply a cutoff in the halo profile for x>xmax, i.e. R>rs*xmax
         See equation 27 of Takada and Jain 2003"""
+        ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7
         g = np.zeros_like(x)
-        xp = np.where(np.logical_and(x>1, x<xmax))
+        xp = np.logical_and(x>1, x<=xmax)
         xo = np.where(np.abs(x-1.) < 10**(-tol))       
-        xm = np.where(x<1)
+        xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        #xm = np.where(x<1)
+        xd = np.where(x<self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        idx = self.theta_amin_to_x(2*1e14, 0.7, 0.12)
         g[xp] = (np.sqrt(xmax**2 - x[xp]**2)/(1+xmax) - 1/np.sqrt(x[xp]**2 - 1) * np.arccos((x[xp]**2 + xmax) / (x[xp] * (1+xmax))) )/(x[xp]**2-1)
         g[xm] = (np.sqrt(xmax**2 - x[xm]**2)/(1+xmax) - 1/np.sqrt(1-x[xm]**2) * np.arccosh((x[xm]**2 + xmax) / (x[xm] * (1+xmax))) )/(x[xm]**2-1)
         g[xo] = np.sqrt(xmax**2 - 1)/(3*(1+xmax)) * (1+ 1/ (1+xmax))
+        g[xd] = (np.sqrt(xmax**2 - idx**2)/(1+xmax) - 1/np.sqrt(1-idx**2) * np.arccosh((idx**2 + xmax) / (idx * (1+xmax))) )/(idx**2-1)
         return g
 
 
@@ -166,20 +175,20 @@ def sigma_int(self, M200, z, R, xmax = 100, npoints=1000):
                 sigma[i] = np.trapz( 2 * r_arr * rho_arr / np.sqrt(r_arr**2 - iR**2), r_arr)
         return sigma
 
-    def analitic_kappa_ft(self, M200, z, ell, const_c=None):
+    def analitic_kappa_ft(self, M200, z, xmax, ell, const_c=None):
         """Analytic Fourier transform of the convergence fiels for a NFW profile
             from Oguri&Takada 2010, Eq.28"""
-        c = self.get_concentration(M200, z, const_c)
-        mu_nfw = np.log(1. + c) - c / (1. + c)
+        c1 = self.get_concentration(M200, z, const_c)
+        mu_nfw = np.log(1. + c1) - c1 / (1. + c1)
         rs = self.get_rs(M200, z)
         chi = self.cosmo.comoving_radial_distance(z)
         k = ell / chi
         x = ((1. + z) * k * rs)
         Six, Cix = sici(x)
-        Sixpc, Cixpc = sici(x * (1. + c))
+        Sixpc, Cixpc = sici(x * (1. + xmax))
         Sidiff = Sixpc - Six
         Cidiff = Cixpc - Cix
-        u0 = np.sin(x) * Sidiff + np.cos(x) * Cidiff - np.sin(x * c) / (x * (1. + c))
+        u0 = np.sin(x) * Sidiff + np.cos(x) * Cidiff - np.sin(x * xmax) / (x * (1. + xmax))
         ufft = 1. / mu_nfw * u0
 
         kappaft = M200 * ufft * (1+z)**2 /chi**2 / self.sigma_crit(z)
@@ -224,6 +233,18 @@ def theta_amin_to_r(self, z, theta_amin):
         R = Dang * theta_amin * np.pi/180/60
         return R
 
+    def theta_amin_to_x(self, M200, z, th_amin, const_c=None):
+        """Angle substended at chararcteric scale x = r/rs
+
+            Args:
+                M200: cluster M200 mass in solar masses (?)
+                z: cluster redshift
+                x:  dimensionless R / Rs
+
+        """
+        return th_amin/self.get_thetas_amin(M200, z, const_c=const_c)
+
+
 
     def x_to_theta_amin(self,M200, z, x, const_c=None):
         """Angle substended at chararcteric scale x = r/rs
@@ -262,9 +283,18 @@ def kappa_map(self, M200, z, shape, lsides, xmax=None):
         dtheta_x = lsides[0]/shape[0] * 180/np.pi*60
         dtheta_y = lsides[1]/shape[1] * 180/np.pi*60
         # Center of the cluster
-        x0 = (shape[0]+1.)/2. 
-        y0 = (shape[1]+1.)/2
-        X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
+        #x0 = shape[0]/2
+        #y0 = shape[1]/2       
+        x0 = 0.
+        y0 = 0.
+
+        #X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
+        if shape[0] % 2 == 0 and shape[1] % 2 == 0:
+            X, Y =  np.meshgrid(np.concatenate((np.arange(0,shape[0]//2), np.arange(-shape[0]//2,0))), np.concatenate((np.arange(0,shape[1]//2), np.arange(-shape[1]//2,0))))
+        else:
+            assert 0, "pixel size is not right"
+        
+        #X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
         theta_amin = self.pix_to_theta(X, Y, (dtheta_x,dtheta_y),  (x0, y0))
 
         return self.kappa_theta(M200, z, theta_amin, xmax=xmax)
diff --git a/lensit/ffs_covs/ffs_cov.py b/lensit/ffs_covs/ffs_cov.py
index 33299d6..1658bd5 100644
--- a/lensit/ffs_covs/ffs_cov.py
+++ b/lensit/ffs_covs/ffs_cov.py
@@ -672,7 +672,7 @@ def apply_conddiagcl(self, typ, alms, use_cls_len=True, use_Pool=0):
     def apply_condpseudiagcl(self, typ, alms, use_Pool=0):
         return self.apply_conddiagcl(typ, alms, use_Pool=use_Pool)
 
-    def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, iblms2=None, **kwargs):
+    def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, iblms2=None,  ellmax_gradleg=None, **kwargs):
         r"""Unormalized quadratic estimates (potential and curl).
 
         Note:
@@ -688,6 +688,7 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, ib
             use_cls_len: use lensed or unlensed cls in QE weights (numerator), defaults to lensed cls
             use_cls_grad: use grad cls in QE weights (numerator), defaults to use_cls_len
             iblms2: second leg inverse variance filtered CMB maps (useful for RDN0 estimates)
+            ellmax_gradleg: cut scales above this value on the gradient leg (used for unbiased cluster lensing reconstruction)
 
         """
         assert iblms.shape == self._skyalms_shape(typ), (iblms.shape, self._skyalms_shape(typ))
@@ -714,6 +715,13 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, ib
                 for _j in range(len(typ)):
                     clms2[_i] += pmat.get_unlPmat_ij(typ, self.lib_skyalm, weights_cls, _i, _j) * iblms2[_j]
 
+        if ellmax_gradleg is not None:
+            fl = np.ones(self.lib_skyalm.ellmax+1)
+            fl[ellmax_gradleg+1:] = 0
+            print("Using grad_cut for T")
+            if 'T' in typ:
+                lib_qlm.almxfl(clms[0], fl, inplace=True)
+
         t.checkpoint("  get_qlms::mult with %s Pmat" % ({True: 'len', False: 'unl'}[use_cls_len]))
 
         _map = lambda alm: self.lib_skyalm.alm2map(alm)
diff --git a/lensit/ffs_iterators/ffs_iterator.py b/lensit/ffs_iterators/ffs_iterator.py
index 939e45e..b346daa 100644
--- a/lensit/ffs_iterators/ffs_iterator.py
+++ b/lensit/ffs_iterators/ffs_iterator.py
@@ -88,10 +88,10 @@ def __init__(self, lib_dir, typ, filt, dat_maps, lib_qlm, Plm0, H0, cpp_prior,
 
         def newton_step_length(it, norm_incr):  # FIXME
             # Just trying if half the step is better for S4 QU
-            if filt.Nlev_uKamin('t') > 2.1: return 1.0
+            if filt.Nlev_uKamin('t') > 2.1: return 0.5
             if filt.Nlev_uKamin('t') <= 2.1 and norm_incr >= 0.5:
-                return 0.5
-            return 0.5
+                return 0.1
+            return 0.1
 
         self.newton_step_length = newton_step_length
         self.soltn0 = soltn0
@@ -662,9 +662,15 @@ def __init__(self, lib_dir, typ, filt, dat_maps, lib_qlm, Plm0, H0, cpp_prior,
         #             'u': (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_sky_ivf + 1)}
         lmax_ivf = filt.lib_datalm.ellmax
         iso_libdat = filt.lib_datalm
-        cls_noise = {'t': (filt.Nlev_uKamin('t') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1),
-                     'q': (filt.Nlev_uKamin('q') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1),
-                     'u': (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1)}
+        cl_t = (filt.Nlev_uKamin('t') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1)
+        cl_q = (filt.Nlev_uKamin('q') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1)
+        cl_u = (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1)        
+        #cl_t[3000:] = 1e9
+        #cl_q[3000:] = 1e9
+        #cl_u[3000:] = 1e9
+        cls_noise = {'t': cl_t,
+                     'q': cl_q,
+                     'u': cl_u}
         self.isocov = ffs_cov.ffs_diagcov_alm(os.path.join(lib_dir, 'isocov'),
                                               iso_libdat, filt.cls, filt.cls, filt.cl_transf, cls_noise,
                                               lib_skyalm=filt.lib_skyalm, init_rank=init_rank,
diff --git a/lensit/ffs_iterators/ffs_iterator_nufft.py b/lensit/ffs_iterators/ffs_iterator_nufft.py
index 0dc4e25..cd17298 100644
--- a/lensit/ffs_iterators/ffs_iterator_nufft.py
+++ b/lensit/ffs_iterators/ffs_iterator_nufft.py
@@ -91,10 +91,10 @@ def __init__(self, lib_dir, typ, filt: filtr, dat_maps, lib_qlm, Plm0, H0, cpp_p
 
         def newton_step_length(it, norm_incr):  # FIXME
             # Just trying if half the step is better for S4 QU
-            if filt.Nlev_uKamin('t') > 2.1: return 1.0
+            if filt.Nlev_uKamin('t') > 2.1: return 0.5
             if filt.Nlev_uKamin('t') <= 2.1 and norm_incr >= 0.5:
-                return 0.5
-            return 0.5
+                return 0.1
+            return 0.1
 
         self.newton_step_length = newton_step_length
         self.soltn0 = soltn0
diff --git a/lensit/misc/misc_utils.py b/lensit/misc/misc_utils.py
index e9f990c..295122f 100644
--- a/lensit/misc/misc_utils.py
+++ b/lensit/misc/misc_utils.py
@@ -119,11 +119,23 @@ def pp_to_kk(ls):
     """
     return ls ** 2 * (ls+1) ** 2 * 0.25 
 
+def p_to_k(ls):
+    """ Converts lensing potential el into lensing convergence el.
+    :math:`C_\ell^{\kappa \kappa} = (\ell (\ell +1))^2 /4 C_\ell^{\phi \phi}`
+    
+    """
+    return ls * (ls+1.) * 0.5 
+
 def kk_to_pp(ls):
     """Converts lensing convergence power spectrum into lensing potential power spectrum.
     """
     return cl_inverse(pp_to_kk(ls))
 
+def k_to_p(ls):
+    """Converts lensing convergence el into lensing potential el.
+    """
+    return cl_inverse(p_to_k(ls))
+
 class timer:
     def __init__(self, verbose, prefix='', suffix=''):
         self.t0 = time.time()
diff --git a/lensit/sims/ffs_maps.py b/lensit/sims/ffs_maps.py
index e0bb310..2d4052e 100644
--- a/lensit/sims/ffs_maps.py
+++ b/lensit/sims/ffs_maps.py
@@ -233,6 +233,19 @@ def _loadQnoise(self):
     def _loadUnoise(self):
         return np.load(self.nUpix)
 
+    def _build_sim_tmap_unl(self, idx):
+        tmap = self.lib_skyalm.almxfl(self.lencmbs.unlcmbs.get_sim_tlm(idx), self.cl_transf)
+        tmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, tmap))
+        return tmap + self.get_noise_sim_tmap(idx)
+
+    def _build_sim_qumap_unl(self, idx):
+        qmap, umap = self.lencmbs.unlcmbs.get_sim_qulm(idx)
+        self.lib_skyalm.almxfl(qmap, self.cl_transf, inplace=True)
+        self.lib_skyalm.almxfl(umap, self.cl_transf, inplace=True)
+        qmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, qmap))
+        umap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, umap))
+        return np.array([qmap + self.get_noise_sim_qmap(idx), umap + self.get_noise_sim_umap(idx)])
+
     def _build_sim_tmap(self, idx):
         tmap = self.lib_skyalm.almxfl(self.lencmbs.get_sim_tlm(idx), self.cl_transf)
         tmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, tmap))
@@ -255,6 +268,28 @@ def get_noise_sim_qmap(self, idx):
     def get_noise_sim_umap(self, idx):
         return self._loadUnoise() * self.pix_pha.get_sim(idx, 2)
 
+    def get_sim_tmap_unl(self, idx):
+        fn = os.path.join(self.lib_dir,  'sim_tmap_unl_%05d.npy'%idx)
+        if self.cache_sims and not os.path.exists(fn):
+            if pbs.rank == 0:
+                np.save(fn, self._build_sim_tmap_unl(idx))
+            pbs.barrier()
+        if self.cache_sims:
+            return np.load(fn)
+        else:
+            return self._build_sim_tmap_unl(idx)
+
+    def get_sim_qumap_unl(self, idx):
+        fn = os.path.join(self.lib_dir,  'sim_qumap_unl_%05d.npy'%idx)
+        if self.cache_sims and not os.path.exists(fn):
+            if pbs.rank == 0:
+                np.save(fn, np.array(self._build_sim_qumap_unl(idx)))
+            pbs.barrier()
+        if self.cache_sims:
+            return np.load(fn)
+        else:
+            return self._build_sim_qumap_unl(idx)
+
     def get_sim_tmap(self, idx):
         fn = os.path.join(self.lib_dir,  'sim_tmap_%05d.npy'%idx)
         if self.cache_sims and not os.path.exists(fn):
diff --git a/notebooks/demo_kappa0.ipynb b/notebooks/demo_kappa0.ipynb
new file mode 100644
index 0000000..9bc1848
--- /dev/null
+++ b/notebooks/demo_kappa0.ipynb
@@ -0,0 +1,1025 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "\n",
+    "import numpy as np \n",
+    "import camb \n",
+    "from numpy.fft import fftshift\n",
+    "from tqdm import tqdm\n",
+    "\n",
+    "import lensit as li\n",
+    "from lensit.clusterlens import lensingmap, profile \n",
+    "from lensit.misc.misc_utils import gauss_beam, pp_to_kk, p_to_k, cl_inverse\n",
+    "from lensit.ffs_covs import ffs_cov, ell_mat\n",
+    "from plancklens.wigners import wigners\n",
+    "from plancklens import n0s, nhl\n",
+    "from plancklens.n1 import n1\n",
+    "\n",
+    "import os\n",
+    "import os.path as op\n",
+    "import matplotlib as mpl\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "from scipy.interpolate import UnivariateSpline as spline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# We start by intiating CAMB which will give us the relevant cosmology \n",
+    "cambinifile = 'planck_2018_acc'\n",
+    "\n",
+    "pars = camb.read_ini(op.join(op.dirname('/Users/sayan/CMB_WORK/CAMB-1.1.3/inifiles'),  'inifiles', cambinifile + '.ini'))\n",
+    "results = camb.get_results(pars)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The size of the data patch is 307.2 X 307.2 arcmin central box\n",
+      "/Users/sayan/Project_Geneva/modules/lensit_sims/temp/clustermaps/camb_planck_2018_acc/cmbexp_S4_sayan/nfw_profile/lensed_by_cluster/npix1024_lpix_0.3amin_lmaxsky6002/M200_2.000000E+14_z0.7/1000_sims\n"
+     ]
+    }
+   ],
+   "source": [
+    "# We define here the parameters for the profile of the cluster\n",
+    "M200, z, xmaxn = 2 * 1e14, 0.7, 3\n",
+    "profname = 'nfw'\n",
+    "key = \"cluster\" # \"lss\"/\"cluster\"/\"lss_plus_cluster\"\n",
+    "profparams={'M200c':M200, 'z':z, 'xmaxn': xmaxn}\n",
+    "hprofile = profile.profile(results, profname)\n",
+    "xmax = xmaxn*hprofile.get_concentration(M200, z)\n",
+    "# Define here the map square patches\n",
+    "npix = 1024  # Number of pixels\n",
+    "lpix_amin = 0.3 # Physical size of a pixel in arcmin (There is bug when <0.2 amin, due to low precision in Cl_TE at )\n",
+    "\n",
+    "print(\"The size of the data patch is %0.1f X %0.1f arcmin central box\"%(npix*lpix_amin, npix*lpix_amin))\n",
+    "\n",
+    "# Maximum multipole used to generate the CMB maps from the CMB power spectra\n",
+    "# ellmaxsky = 6000 # (bug when ellmax>6300 because of low precision in Cl_TE of CAMB )\n",
+    "ellmaxsky = 6002 \n",
+    "\n",
+    "# Set the maximum ell observed in the CMB data maps\n",
+    "ellmaxdat = 4000\n",
+    "ellmindat = 100\n",
+    "\n",
+    "# Number of simulated maps \n",
+    "nsims = 1000\n",
+    "\n",
+    "# Set CMB experiment for noise level and beam\n",
+    "cmb_exp='S4_sayan'\n",
+    "\n",
+    "# We will cache things in this directory \n",
+    "libdir = lensingmap.get_cluster_libdir(cambinifile,  profname, key, npix, lpix_amin, ellmaxsky, M200, z, nsims, cmb_exp)\n",
+    "#libdir = op.join(libdir,\"trunc\")\n",
+    "print(libdir)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lmax = ellmaxsky\n",
+    "cpp_fid = results.get_lens_potential_cls(lmax=lmax, raw_cl=True).T[0]\n",
+    "\n",
+    "camb_cls = results.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=lmax).T\n",
+    "cls_unl_fid = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid}\n",
+    "\n",
+    "camb_cls_len = results.get_lensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=lmax).T\n",
+    "cls_len_fid = {'tt':camb_cls_len[0], 'ee':camb_cls_len[1], 'bb':camb_cls_len[2], 'te':camb_cls_len[3], 'pp':cpp_fid}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/sayan/Project_Geneva/git_repos/LensIt/lensit/clusterlens/profile.py:320: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  dx_lm = 2 * rfft_kappa * 1.j * KX / (KX**2+KY**2)\n",
+      "/Users/sayan/Project_Geneva/git_repos/LensIt/lensit/clusterlens/profile.py:321: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  dy_lm = 2 * rfft_kappa * 1.j * KY / (KX**2+KY**2)\n"
+     ]
+    }
+   ],
+   "source": [
+    "np.random.seed(seed=20)\n",
+    "clustermaps = lensingmap.cluster_maps(libdir, key, npix, lpix_amin, nsims, results, profparams, profilename=profname,  ellmax_sky = ellmaxsky, ellmax_data=ellmaxdat, ellmin_data=ellmindat, cmb_exp=cmb_exp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ellmax_sky = clustermaps.ellmax_sky\n",
+    "sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(clustermaps.cmb_exp)\n",
+    "\n",
+    "cls_noise = {'t': (sN_uKamin * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1),\n",
+    "            'q':(sN_uKaminP * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1),\n",
+    "            'u':(sN_uKaminP * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1)}  # simple flat noise Cls\n",
+    "# cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=ellmax_sky)\n",
+    "# lib_alm = ell_mat.ffs_alm_pyFFTW(get_ellmat(LD_res, HD_res=HD_res),\n",
+    "                    # filt_func=lambda ell: (ell >= ellmin) & (ell <= ellmax), num_threads=pyFFTWthreads)\n",
+    "# lib_skyalm = ell_mat.ffs_alm_pyFFTW(clustermaps.ellmat,\n",
+    "                    # filt_func=lambda ell: (ell <= ellmax_sky), num_threads=clustermaps.num_threads)\n",
+    "\n",
+    "cl_transf = clustermaps.cl_transf\n",
+    "lib_skyalm = clustermaps.lib_skyalm\n",
+    "\n",
+    "typ = 'T'\n",
+    "\n",
+    "lib_dir = op.join(clustermaps.dat_libdir, typ)\n",
+    "# isocov = ffs_cov.ffs_diagcov_alm(lib_dir, clustermaps.lib_datalm, clustermaps.cls_unl, cls_len, cl_transf, cls_noise, lib_skyalm=lib_skyalm)\n",
+    "isocov = ffs_cov.ffs_diagcov_alm(lib_dir, clustermaps.lib_datalm, cls_unl_fid, cls_unl_fid, cl_transf, cls_noise, lib_skyalm=lib_skyalm)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:00])   inverse len Pmats  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:01]  (total [00:00:01])   Fxx , part 1  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:01])   Fyy , part 1  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:01])   Fxy , part 1  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:01])   Fxx , part 2  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:01])   Fyy , part 2  curvpOlm \n",
+      "  lensit.ffs_covs.ffs_cov   [00:00:00]  (total [00:00:01])   Fxy , part 2  curvpOlm \n"
+     ]
+    }
+   ],
+   "source": [
+    "H0len =  isocov.get_response(typ, lib_skyalm, use_cls_len=True)[0]\n",
+    "def get_starting_point(idx, typ, clustermaps, keyword=None): \n",
+    "    \"\"\"\n",
+    "    This returns initial data for simulation index 'idx' from a CMB-S4 simulation library.\n",
+    "    On first call the simulation library will generate all simulations phases, hence might take a little while.\n",
+    "    \"\"\" \n",
+    "\n",
+    "    print(\" I will be using data from ell=%s to ell=%s only\"%(isocov.lib_datalm.ellmin, isocov.lib_datalm.ellmax))\n",
+    "    print(\" The sky band-limit is ell=%s\"%(isocov.lib_skyalm.ellmax))\n",
+    " \n",
+    "    lib_qlm = lib_skyalm #: This means we will reconstruct the lensing potential for all unlensed sky modes.\n",
+    "    ellmax_sky = lib_skyalm.ellmax\n",
+    "    ell = np.arange(ellmax_sky+1)\n",
+    "\n",
+    "    if typ=='QU':\n",
+    "        datalms = np.array([isocov.lib_datalm.map2alm(m) for m in clustermaps.maps_lib.get_sim_qumap(idx)]) \n",
+    "    elif typ =='T':\n",
+    "        datalms = np.array([isocov.lib_datalm.map2alm(clustermaps.maps_lib.get_sim_tmap(idx))]) \n",
+    "    elif typ =='TQU':\n",
+    "        datalms = np.array([isocov.lib_datalm.map2alm(m) for m in np.array([clustermaps.maps_lib.get_sim_tmap(idx), clustermaps.maps_lib.get_sim_qumap(idx)[0], clustermaps.maps_lib.get_sim_qumap(idx)[1]])]) \n",
+    "   \n",
+    "    use_cls_len = True\n",
+    " \n",
+    "    plm1 = 0.5 * isocov.get_qlms(typ,  isocov.get_iblms(typ, datalms, use_cls_len=use_cls_len)[0], lib_qlm, \n",
+    "                                 use_cls_len=use_cls_len)[0]\n",
+    "    \n",
+    "    plmqe = lib_skyalm.almxfl(plm1, cl_inverse(H0len), inplace=False)\n",
+    "\n",
+    "    return  plmqe"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "1\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "2\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "3\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "4\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "5\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "6\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "7\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "8\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "9\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "10\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "11\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "12\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "13\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "14\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "15\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "16\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "17\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "18\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "19\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "20\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "21\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "22\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "23\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "24\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "25\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "26\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "27\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "28\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "29\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "30\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "31\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "32\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "33\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "34\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "35\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "36\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "37\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "38\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "39\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "40\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "41\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "42\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "43\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "44\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "45\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "46\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "47\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "48\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "49\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "50\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "51\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "52\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "53\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "54\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "55\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "56\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "57\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "58\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "59\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "60\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "61\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "62\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "63\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "64\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "65\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "66\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "67\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "68\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "69\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "70\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "71\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "72\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "73\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "74\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "75\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "76\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "77\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "78\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "79\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "80\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "81\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "82\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "83\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "84\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "85\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "86\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "87\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "88\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "89\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "90\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "91\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "92\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "93\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "94\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "95\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "96\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "97\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "98\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n",
+      "99\n",
+      " I will be using data from ell=140 to ell=3999 only\n",
+      " The sky band-limit is ell=6002\n",
+      "Building nufft fwd plan\n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::mult with len Pmat  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::cartesian gradients  \n",
+      "     [00:00:00]  (total [00:00:00])   get_qlms::rotation to phi Omega  \n"
+     ]
+    }
+   ],
+   "source": [
+    "nmaps = 100\n",
+    "plmqes = [None]*nmaps\n",
+    "\n",
+    "if nsims >1:\n",
+    "    for idx in range(nmaps):\n",
+    "        print(idx)\n",
+    "        plmqes[idx] = get_starting_point(idx, typ, clustermaps)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 100/100 [00:00<00:00, 133.07it/s]\n"
+     ]
+    }
+   ],
+   "source": [
+    "auto_corr_1 = nmaps*[None]\n",
+    "for idx in tqdm(range(nmaps)):\n",
+    "    auto_corr_1[idx] = lib_skyalm.alm2cl(plmqes[idx])\n",
+    "\n",
+    "phi_var_1 = np.mean(auto_corr_1, axis=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 100/100 [00:01<00:00, 88.48it/s]\n"
+     ]
+    }
+   ],
+   "source": [
+    "kappa_0_sims_QE_T_1 = [None]*nmaps\n",
+    "lmax, lmin= 5000, ellmindat\n",
+    "ell = np.arange(ellmaxsky+1)\n",
+    "for idx in tqdm(range(nmaps)):\n",
+    "    kappa_0_sims_QE_T_1[idx] = clustermaps.get_kappa0_from_sim(lmin=lmin, lmax=lmax, phi_obs_lm=plmqes[idx], NL=phi_var_1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "========================QE Results========================\n",
+      "input kappa_0 0.128488\n",
+      "mean kappa 0 from 100 simulations is 0.120446\n",
+      "error from 100 simulations 0.019439\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"========================QE Results========================\")\n",
+    "print(\"input kappa_0 %f\"%(hprofile.get_kappa0(M200, z, xmax)))\n",
+    "print(\"mean kappa 0 from %i simulations is %f\"%(nmaps,np.mean(kappa_0_sims_QE_T_1)))\n",
+    "print(\"error from %i simulations %f\"%(nmaps,np.std(kappa_0_sims_QE_T_1)/np.sqrt(nmaps)))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.9.13 ('test')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "2a0a0d92d61eb81c1bfacf4848f33179e421b61b01830eeaa00bb11be584f4b6"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}