""" clplot.py

Make a plot. 
"""

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

### Gather all sim_ids
#multiplex cross_prod group_by '0'

### Define some global parameters
lst_par = ['deltal', 'lmin_ratio']
lst_tag = lst_par

### Retrieve some global parameters
load_param('cmb',   globals(), ["input_cl", "lmax", "cmb_unit"])
load_param("noise", globals(), ["noise_unit", "noise_power"])

## load cls
inp_cls = pl.loadtxt(input_cl)[0:lmax+1,0]
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)


## make a plot
smooth_win = 0.05
ell = pl.arange(lmax+1)
pl.figure()
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")

pl.savefig(get_data_fn("cls.png"))
pl.savefig(get_data_fn("cls.pdf"))

# 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"))