Commit 82a9b633 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

demo cmb

parent bc31724f
""" This is the docstring of the segment script.
This segment script illustrates how to use the various utilities
available from the default environment.
"""
### How to organize the segment script ?
###
### A segment script contains some python code which will be applied on
### some data set.
### The pipelet software has no demand on the code, but it
### offers a lot of utilities which may set the user straight
### about some good practice.
###
### From the pipeline point of view the segment script corresponds to
### some processing applied to the data (input/output) with some tuned parameters.
### The pipelet software offers some default utilities with respect of that
### 3 entities :
### - processing,
### - data,
### - parameters.
### The data and files utilities :
###
### Data are usually stored on files.
### The pipelet software offers a specific location for all data files
### that can be parsed from the web interface.
### To get one on those files :
output_fn = get_data_fn ('mydata.dat')
### To retrieve some data file from an upstream segment named 'first':
input_fn = glob_seg ('*.dat', 'first')[0]
### One may need some temporary files also :
temporary_fn = get_tmp_fn()
### To read and write data to files, one may use some default routines
### based on the pickle python module :
load_products (input_fn , globals(), ["D"])
save_products (output_fn, globals(), ["D"])
### The parameters utilities :
###
### The pipelet software saves automatically all parameters listed in
### the variable named 'lst_par'.
### Some of them can be made visible from the interface using the
### variable 'lst_tag'
lst_par = ['my_dim', 'my_option', 'my_nb_of_iter']
lst_tag = ['my_dim', 'my_nb_of_iter']
### To retrieve some parameters of an upstream segment :
load_param ('first', globals(), ['my_previous_option'])
### The processing utilities :
###
### A standard logging.Logger object is available from the segment
### context.
log.info("Until now, nothing's wrong\n")
### It may occur that the processing applied to the data contains some
### generic code you want to recycle from one segment to another.
### This portion of code can be factorized in a hook script names
### 'seg_segname_hookname.py' which can be called with:
hook ('preproc', globals())
### It is also possible to call an arbitrary subprocess without
### loosing the log facility:
logged_subprocess(['my_sub_process', 'my_dim', 'my_nb_of_iter'])
### The pipe scheme utilities :
###
### The pipe scheme can be controlled from the segment environment
### also. The default behaviour is set by the segment input and
### output variables :
seg_output = [1,2,3] ## will generate 3 instances of the downstream
## segment. Each of them will receive one element of
## the list as input.
seg_output = seg_input ## will generate 1 instance of the downstream
## segment, per current instance.
seg_output = None ## will generate 1 instance of the downstream
## segment for all current instance.
### This default behavious can be altered by specifying an @multiplex
### directive (see documentation).
""" clplot.py
Correct cls from mask effect.
Make a plot.
"""
import healpy as hp
import pylab as pl
import spherelib as sp
## Gather all cls
#multiplex cross_prod group_by 0
load_param('mask', globals(), ["mask_types", "lmax", "fsky"])
load_param('cmb', globals(), ["input_cl"])
## load cls
apo_cls= pl.loadtxt(glob_seg('cls', 'apodized*.txt')[0])
bin_cls= pl.loadtxt(glob_seg('cls', 'binary*.txt')[0])
## load mll'
mll_fn = glob_seg('mask', 'mll*')[0]
load_products (mll_fn , globals(), ["mll"])
## corrected cls
binc_cls = pl.linalg.solve(mll , bin_cls)
bin_cls = bin_cls/fsky
apo_cls = apo_cls/fsky
## smoothed cls
dloverl = 0.05
binc_cls = sp.smooth_cl(binc_cls, dloverl=dloverl)
bin_cls = sp.smooth_cl(bin_cls, dloverl=dloverl)
apo_cls = sp.smooth_cl(apo_cls, dloverl=dloverl)
## input cls
inp_cls = pl.loadtxt(input_cl)[0:lmax+1,0]
## make a plot
ell = pl.arange(lmax+1)
pl.figure()
pl.plot(ell, ell*(ell+1)*bin_cls/(2*pl.pi))
pl.plot(ell, ell*(ell+1)*binc_cls/(2*pl.pi))
pl.plot(ell, ell*(ell+1)*apo_cls/(2*pl.pi))
pl.plot(ell, ell*(ell+1)*inp_cls/(2*pl.pi), color="black")
pl.legend([ "fsky", "mll", "apodized mask", "input"], loc="upper right")
pl.savefig(get_data_fn("cls.png"))
pl.savefig(get_data_fn("cls.pdf"))
""" cls.py
Apply the mask to the cmb map
Compute the power spectrum
Correct it from the mode coupling matrix
Same for apodized mask
"""
import healpy as hp
import pylab as pl
### Retrieve some global parameters
load_param("cmb", globals(), ["lmax"])
seed = seg_input["cmb"]
mask_type = seg_input["mask"]
### Get inputs
cmb_fn = glob_seg("cmb", "map*.fits")[0]
cmb_map = hp.read_map(cmb_fn)
mask_fn = glob_seg("mask", "%s*.fits"%mask_type)[0]
mask = hp.read_map(mask_fn)
### Mask the input map
cmb_mask = cmb_map*mask
### Compute power spectrum
cl = hp.anafast (cmb_mask, lmax)
### save it to disk
cl_fn = get_data_fn("%s_cls.txt"%mask_type)
pl.savetxt(cl_fn, cl)
### Set output
seg_output = [mask_type]
""" cmb.py
Generate a cmb map from lambda-CDM power spectrum.
"""
import healpy as hp
import pylab as pl
### Define some global parameters
lst_par = ['lmax', 'nside', 'cmb_unit', 'seed', 'input_cl']
lst_tag = lst_par[:-2]
nside = seg_input.values()[0]
lmax = 2*nside
cmb_unit = "uK_CMB"
seed = 0 ## actually not use by synfast
### Generate a cmb map
input_cl = "lambda_best_fit.txt"
cmb_cl = pl.loadtxt(input_cl)[0:lmax+1,0]
cmb_map = hp.synfast(cmb_cl, nside, lmax=lmax)
### Save it to disk
cmb_cl_fn = get_data_fn ('cls_cmb.txt')
cmb_map_fn = get_data_fn ('map_cmb.fits')
hp.write_map(cmb_map_fn, cmb_map)
pl.savetxt (cmb_cl_fn , cmb_cl)
### Make a plot
cmb_map_fig = cmb_map_fn.replace('.fits', '.png')
hp.mollview(cmb_map)
pl.savefig (cmb_map_fig)
### Set output
seg_output = [seed]
""" mask.py
Compute the mode coupling matrix of an arbitraty binary mask
Compute the apodized version of this mask
"""
import healpy as hp
import pylab as pl
import spherelib as sp
### Define some global parameters
lst_par = ['lmax', 'nside', 'radius', 'fsky', 'mask_types']
lst_tag = lst_par[:-2]
nside = seg_input.values()[0]
lmax = 2*nside
radius = 60 ## arcmin
mask_types = ["binary", "apodized"]
### Read the input mask
mask = hp.read_map("WG2_large_galactic_mask_256.fits")
mask = hp.ud_grade(mask, nside) ## adapt to nside
mask[mask!=1] = 0 ## make a binary mask
### Save it to disk
mask_fn = get_data_fn("%s_mask_%d.fits"%(mask_types[0],nside))
hp.write_map(mask_fn, mask)
### Compute mll'
mll = sp.mask2mll (mask, lmax)
mll_fn = get_data_fn("mll.dat")
save_products (mll_fn, globals(), ["mll"])
### Apodizing
sp.apodize_mask (mask, radius)
fsky = pl.mean(mask**2)
mask_fn = get_data_fn("%s_mask_%d.fits"%(mask_types[1],nside))
hp.write_map(mask_fn, mask)
### Make a plot
mask_fig = mask_fn.replace('.fits', '.png')
hp.mollview(mask)
pl.savefig (mask_fig)
### Set output
seg_output = mask_types
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment