diff --git a/analysis.py b/analysis.py index 598536c..3d3405f 100644 --- a/analysis.py +++ b/analysis.py @@ -39,11 +39,11 @@ def __init__(self, self.debug = debug - # Load template - thetas,kappa_1h,kappa_2h = np.loadtxt(c.template_file,unpack=True) - kappa = kappa_1h + kappa_2h - self.thetas = thetas - self.template_kappa_1d = kappa + # Load template for cross_correlations + thetas_cross,kappa_1h_cross,kappa_2h_cross = np.loadtxt(c.template_file,unpack=True) + kappa_cross = kappa_1h_cross + kappa_2h_cross + self.thetas_cross = thetas_cross + self.kappa_cross = kappa_cross if atype=='flatsky_sim': # If we are doing a padded high-res flat sim @@ -53,9 +53,24 @@ def __init__(self, # will be slightly different from that constructed # directly - self.csim = olensing.FixedLens(thetas, kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) + # Make template profile for lensing + mass = c.mass + z=c.z + delta = c.delta + + thetas,kappa_1h,kappa_2h,_,_,_,_,_ = olensing.kappa_nfw_profiley(mass=mass,conc=None, + z=z,z_s=1100.,background='critical',delta=delta, R_off_Mpc = None, + apply_filter=False) + + kappa = kappa_1h + kappa_2h + self.thetas = thetas + self.kappa = kappa + + self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact, + skip_lensing=c.skip_lensing) _,_,dummy = self.csim.generate_sim(0) # FIXME: is this necessary self.shape, self.wcs = dummy.shape, dummy.wcs + else: if c.periodic: raise ValueError self.thumbnail_r = c.width_deg * u.deg / 2. @@ -66,7 +81,7 @@ def __init__(self, self.rbinner = stats.bin2D(self.modrmap,self.rbin_edges) self.rcents = self.rbinner.cents self.modlmap = enmap.modlmap(self.shape,self.wcs) - self.template = enmap.enmap(maps.interp(thetas,kappa)(self.modrmap),self.wcs) + self.template = enmap.enmap(maps.interp(self.thetas_cross,self.kappa_cross)(self.modrmap),self.wcs) self.reconstructor = Recon(self.shape,self.wcs, @@ -97,9 +112,11 @@ def __init__(self, self.Lmin = c.Lmin self.Lmax = c.Lmax + # Auto-correlation for template self.ktemplate = enmap.fft(self.template * self.ctaper,normalize='phys') cents,p1d,p2d = self.power(self.ktemplate,self.ktemplate) self.template_auto_p1d = p1d + if debug and self.rank==0: io.plot_img(self.ctaper,self._out('ctaper.png'),arc_width=c.width_deg*60., xlabel='$\\Theta$ (arcmin)',ylabel='$\\Theta$ (arcmin)') @@ -110,7 +127,24 @@ def __init__(self, pl.hline(y=0) pl.done(self._out('template_auto_p1d.png')) - + # Cross-correlation of template and lensing template + self.lensing_template = enmap.enmap(maps.interp(self.thetas,self.kappa)(self.modrmap),self.wcs) + self.klensing = enmap.fft(self.lensing_template * self.ctaper,normalize='phys') + cents, p1d_cross, p2d_cross = self.power(self.ktemplate, self.klensing) + self.template_cross_p1d = p1d_cross + + if debug and self.rank==0: + cents, p1d_auto, p2d_auto = self.power(self.klensing, self.klensing) + + io.plot_img((np.fft.fftshift(self.modlmap**2 * p2d_cross)),self._out('template_cross_p2d.png')) + pl = io.Plotter(xyscale='linlin',xlabel='L',ylabel='L^2 C_L') + pl.add(cents,cents**2 * p1d_cross,marker='o',label='cross') + pl.add(cents,cents**2 * p1d_auto,marker='o',label='lensing template auto') + pl.add(cents,cents**2 * p1d,marker='o',label='correlation template auto') + pl.vline(x=500) + pl.hline(y=0) + pl.done(self._out('template_cross_p1d.png')) + self.fit_Kmask = maps.mask_kspace(self.shape, self.wcs, lmin=c.Lmin, lmax=c.Lmax) if self.rank==0: @@ -238,7 +272,9 @@ def finish(self, save_transfer=False): lcents = self.lcents Lmin = self.Lmin Lmax = self.Lmax - p1d = self.template_auto_p1d + if self.atype=='flatsky_sim': + p1d = self.template_cross_p1d + else: p1d = self.template_auto_p1d cmean = self.s.stats['cp1d']['mean'] cerr = self.s.stats['cp1d']['errmean'] ccov = self.s.stats['cp1d']['covmean'] @@ -306,9 +342,10 @@ def finish(self, save_transfer=False): pl = io.Plotter() pl.add_err(rcents[rsel]/u.arcmin,rmean,yerr=rerr,marker='o') - pl.add(self.thetas/u.arcmin, self.template_kappa_1d) + pl.add(self.thetas/u.arcmin, self.kappa) #FIXME: only for flatsky sim? pl._ax.set_xlim(0,10) - pl._ax.set_ylim(-0.01,rmean.max()*1.2) + pl._ax.set_ylim(rmean.min()-0.1,rmean.max()+0.1) + # pl._ax.set_ylim(-0.01,rmean.max()*1.2) pl.hline(y=0) pl.done(self._out('recon_profile.png')) diff --git a/choices.yaml b/choices.yaml index b1efdb0..38717f4 100644 --- a/choices.yaml +++ b/choices.yaml @@ -42,7 +42,13 @@ ybeam_noise_uk_arcmin: 10.0 # ================================================ # How the template and reconstruction are treated # ================================================ +# Set whether lensing will be used in the flatsky sim case +skip_lensing: False +# This is the profile used for lensing reconstruction in the flatsky sim case +mass: 2.0e+13 +delta: 500 #'200' or '500' for critical mass (profiley will calculate m200c or m500c) +z: 1.0 theta_max_arc: 20.0 # All templates and reconstructions are masked beyond this radius so that we are not sensitive to the modeling there roll_arc: 20.0 # with a cosine taper that rolls of to zero in this distance # The above two decide what ctaper looks like (saved to ctaper.png) diff --git a/kapstack.py b/kapstack.py index 49b386b..28b7578 100644 --- a/kapstack.py +++ b/kapstack.py @@ -93,4 +93,4 @@ pl.done('profile_rat.png') -enmap.write_map("stack.fits",out) +enmap.write_map("stack.fits",out) \ No newline at end of file