clcmb.py 1.49 KB
Newer Older
Maude Le Jeune's avatar
Maude Le Jeune committed
1 2 3 4
""" clcmb.py

Coadd cmb and noise map.
Compute the power spectrum 
Maude Le Jeune's avatar
Maude Le Jeune committed
5
Correct it from the noise bias
Maude Le Jeune's avatar
Maude Le Jeune committed
6 7 8
Correct it from the mode coupling matrix
"""

Maude Le Jeune's avatar
Maude Le Jeune committed
9 10 11 12
### Gather inputs by sim_id
#multiplex cross_prod group_by 'cmb'
sim_id   = seg_input 

Maude Le Jeune's avatar
Maude Le Jeune committed
13 14 15 16
import healpy as hp
import pylab as pl
import spherelib as sp

Maude Le Jeune's avatar
Maude Le Jeune committed
17 18 19 20
### Define some global parameters
lst_par = ['sim_id']
lst_tag = lst_par

Maude Le Jeune's avatar
Maude Le Jeune committed
21 22
### Retrieve some global parameters
load_param("cmb", globals(), ["lmax", "cmb_unit"])
Maude Le Jeune's avatar
Maude Le Jeune committed
23
load_param("noise", globals(), ["noise_unit", "noise_power"])
Maude Le Jeune's avatar
Maude Le Jeune committed
24 25

### Get inputs
Maude Le Jeune's avatar
Maude Le Jeune committed
26 27
cmb_fn   = glob_seg("cmb", "map*.fits")[0]
cmb_map  = hp.read_map(cmb_fn)
Maude Le Jeune's avatar
Maude Le Jeune committed
28 29
noise_fn = glob_seg("noise", "map*.fits")[0]
noise    = hp.read_map(noise_fn)
Maude Le Jeune's avatar
Maude Le Jeune committed
30 31 32

### Coadd cmb and noise maps
cv_factor= sp.astro.convfact(fr=noise_unit, to=cmb_unit)
Maude Le Jeune's avatar
Maude Le Jeune committed
33
noise    = noise*cv_factor
Maude Le Jeune's avatar
Maude Le Jeune committed
34 35 36 37 38 39
cmb_map  = cmb_map + noise

### Make a plot
coadd_fn = get_data_fn("noisy-cmb.png")
hp.mollview(cmb_map, title="noisy cmb in %s"%cmb_unit)
pl.savefig(coadd_fn)
Maude Le Jeune's avatar
Maude Le Jeune committed
40 41 42 43


### Compute power spectrum
cl = hp.anafast (cmb_map, lmax)
Maude Le Jeune's avatar
Maude Le Jeune committed
44
cl = cl - (noise_power*cv_factor*cv_factor) ## debias
Maude Le Jeune's avatar
Maude Le Jeune committed
45

Maude Le Jeune's avatar
Maude Le Jeune committed
46
### save unbiased cl to disk
Maude Le Jeune's avatar
Maude Le Jeune committed
47 48 49 50 51 52 53
cl_fn = get_data_fn("cls.txt")
pl.savetxt(cl_fn, cl)

## load mll'
mll_fn = glob_seg('noise', 'mll*')[0]
load_products (mll_fn , globals(), ["mll"])

Maude Le Jeune's avatar
Maude Le Jeune committed
54
## correct cl from coupling matrix
Maude Le Jeune's avatar
Maude Le Jeune committed
55
cl = pl.linalg.solve(mll , cl)
Maude Le Jeune's avatar
Maude Le Jeune committed
56
cl = cl * sum(mll[0,:])
Maude Le Jeune's avatar
Maude Le Jeune committed
57

Maude Le Jeune's avatar
Maude Le Jeune committed
58
### save mode couplit corrected cl to disk
Maude Le Jeune's avatar
Maude Le Jeune committed
59 60 61
cl_fn = get_data_fn("cls-mll.txt")
pl.savetxt(cl_fn, cl)

Maude Le Jeune's avatar
Maude Le Jeune committed
62 63
### Set output : forward as many childs as sim ids
seg_output = [sim_id] 
Maude Le Jeune's avatar
Maude Le Jeune committed
64 65