clplot.py 2.59 KB
Newer Older
Maude Le Jeune's avatar
Maude Le Jeune committed
1 2 3 4 5 6 7 8 9
""" clplot.py

Make a plot. 
"""

import healpy as hp
import pylab as pl
import spherelib as sp

Maude Le Jeune's avatar
Maude Le Jeune committed
10 11
### Gather all sim_ids
#multiplex cross_prod group_by '0'
Maude Le Jeune's avatar
Maude Le Jeune committed
12

Maude Le Jeune's avatar
Maude Le Jeune committed
13 14 15
### Define some global parameters
lst_par = ['deltal', 'lmin_ratio']
lst_tag = lst_par
Maude Le Jeune's avatar
Maude Le Jeune committed
16

Maude Le Jeune's avatar
Maude Le Jeune committed
17 18 19
### Retrieve some global parameters
load_param('cmb',   globals(), ["input_cl", "lmax", "cmb_unit"])
load_param("noise", globals(), ["noise_unit", "noise_power"])
Maude Le Jeune's avatar
Maude Le Jeune committed
20

Maude Le Jeune's avatar
Maude Le Jeune committed
21
## load cls
Maude Le Jeune's avatar
Maude Le Jeune committed
22
inp_cls = pl.loadtxt(input_cl)[0:lmax+1,0]
Maude Le Jeune's avatar
Maude Le Jeune committed
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
try:
    cl_n = pl.loadtxt(glob_seg('clnoise', 'cl*.txt')[0])
except:
    cl_n = pl.ones((lmax+1))*noise_power

## Convert noise in cmb unit
cv_factor= sp.astro.convfact(fr=noise_unit, to=cmb_unit)
cl_n = cl_n*cv_factor**2

## Get mean cl values and error bars
cl_unbs  = glob_seg('clcmb', 'cls.txt')
cl_mlls  = glob_seg('clcmb', 'cl*mll*.txt')

## Binning
deltal = 10
bins   = sp.uniformbin(0, lmax, deltal)
lmean  = pl.squeeze(sp.get_lmean(bins))
nbins  = len(lmean)
nsims    = len(cl_unbs)
cl_unb   = pl.zeros((nbins, nsims))
cl_mll   = pl.zeros((nbins, nsims))
for sim_id in range(nsims):
    cl_unb[:,sim_id] = sp.bin_cl(pl.loadtxt(cl_unbs[sim_id]), bins)
    cl_mll[:,sim_id] = sp.bin_cl(pl.loadtxt(cl_mlls[sim_id]), bins)
er_unb = pl.std(cl_unb, axis=1)
er_mll = pl.std(cl_mll, axis=1)
cl_unb = pl.mean(cl_unb, axis=1)
cl_mll = pl.mean(cl_mll, axis=1)

Maude Le Jeune's avatar
Maude Le Jeune committed
52 53

## make a plot
Maude Le Jeune's avatar
Maude Le Jeune committed
54
smooth_win = 0.05
Maude Le Jeune's avatar
Maude Le Jeune committed
55 56
ell = pl.arange(lmax+1)
pl.figure()
Maude Le Jeune's avatar
Maude Le Jeune committed
57 58 59 60 61 62 63 64 65 66 67 68 69 70
pl.errorbar(lmean, lmean*(lmean+1)*cl_unb/(2*pl.pi), yerr=lmean*(lmean+1)*er_unb /(2*pl.pi),color="b")
pl.errorbar(lmean, lmean*(lmean+1)*cl_mll/(2*pl.pi), yerr=lmean*(lmean+1)*er_mll /(2*pl.pi),color="c" )
pl.plot(ell, ell*(ell+1)*cl_n /(2*pl.pi),color="r")
pl.plot(ell, ell*(ell+1)*inp_cls /(2*pl.pi),color="k")
pl.ylabel("$\ell (\ell+1) c_\ell / 2 \pi$", fontsize=18)
pl.xlabel("$\ell$", fontsize=18)
pl.title("Power spectra in %s$^2$"%cmb_unit)

line1 = pl.Line2D(lmean, lmean, linestyle='-', color='b',linewidth=1)
line2 = pl.Line2D(lmean, lmean, linestyle='-',color='c', linewidth=1)
line3 = pl.Line2D(lmean, lmean, linestyle='-',color='r', linewidth=1)
line4 = pl.Line2D(lmean, lmean, linestyle='-',color='k', linewidth=1)
pl.legend((line1, line2,line3,line4),["unbiased", "mll", "noise", "input"], loc="upper center")

Maude Le Jeune's avatar
Maude Le Jeune committed
71 72
pl.savefig(get_data_fn("cls.png"))
pl.savefig(get_data_fn("cls.pdf"))
Maude Le Jeune's avatar
Maude Le Jeune committed
73 74 75 76 77 78 79 80 81 82

# lmin_ratio = 200
# pl.figure()
# pl.plot(pl.arange(lmin_ratio,lmax+1), sp.smooth_cl(cl_unb[lmin_ratio:]/cl_mll[lmin_ratio:],dloverl=smooth_win))
# pl.legend([ "unbiased/mll"], loc="upper center")
# pl.ylabel("$c_\ell^{unb} /c_\ell^{mll} $", fontsize=18)
# pl.xlabel("$\ell$", fontsize=18)
# pl.title("Power spectra ratio")
# pl.savefig(get_data_fn("ratio.png"))