True merging of the branch fit_p1d into the master branch

.gitignore

@@ -1,2 +1,3 @@
 **/*.pyc
 *~
+etc/pkmiss_interp.fits

analysis/bias_beta_vs_rmin.py

@@ -0,0 +1,119 @@
+import glob
+import h5py
+import fitsio
+import numpy as np
+import matplotlib.pyplot as plt
+import picca.wedgize
+from SaclayMocks import util
+
+
+'''
+This code plots the evolution of bias and beta when we change rmin in the fitter
+with all the other parameters fixed to the value found with the standard fit
+'''
+
+SMALL_SIZE = 18
+MEDIUM_SIZE = 20
+BIGGER_SIZE = 22
+plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
+plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
+plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
+plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
+
+plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
+plt.rc('figure', figsize=(9,7))
+
+# Parameters
+# indir = "/global/cscratch1/sd/tetourne/Analysis/redo_dr16/Fits/cf/Bias_beta_vs_rmin/z_0_10"
+# indir = "/global/cscratch1/sd/tetourne/Out/mock_0.2/Fits"
+indir = "/global/cscratch1/sd/tetourne/Out/v4.7.22/from_transmissions/Fit_pred"
+# title = 'dr16'
+title = 'prediction'
+
+files = np.sort(glob.glob(indir+"/*0.h5"))
+print(files)
+rmin_list = []
+bias_eta_list = []
+bias_eta_err_list = []
+beta_list = []
+beta_err_list = []
+bias_list = []
+bias_err_list = []
+beff_list = []
+beff_err_list = []
+
+for f in files:
+    pars = util.extract_h5file(f)
+    i = f.rfind('rmin')
+    rmin_list.append(np.int32(f[i+4:i+6]))
+    bias_eta = pars[2]['bias_eta_LYA']
+    beta = pars[2]['beta_LYA']
+    f = pars[2]['growth_rate']
+    bias_eta_err = pars[3]['bias_eta_LYA']
+    beta_err = pars[3]['beta_LYA']
+    cov = pars[4]['cov[beta_LYA, bias_eta_LYA]']
+    bias = bias_eta * f / beta
+    # bias_err = np.sqrt((f*bias_eta_err/beta)**2 + (bias_eta*f*beta_err/beta**2)**2 - 2*cov/bias_eta/beta)
+    bias_err = np.sqrt((f*bias_eta_err/beta)**2 + (bias_eta*f*beta_err/beta**2)**2 - 2*cov*bias_eta*f**2/beta**3)
+    beff = bias_eta*f * np.sqrt(1+2/3*beta+1/5*beta**2) / beta
+    # db_dbiaseta = f*np.sqrt(1+2/3*beta+1/5*beta**2)/beta
+    # db_dbeta = (bias_eta*f/beta)*np.sqrt(1+2/3*beta+1/5*beta**2)*((1/3+1/5*beta)/(1+2/3*beta+1/5*beta**2) - 1/beta)
+    # beff_err = np.sqrt((db_dbiaseta*bias_eta_err)**2 + (db_dbeta*beta_err)**2 + 2*db_dbiaseta*db_dbeta*cov)
+    beff_err = util.beff_err(bias_eta, bias_eta_err, beta, beta_err, cov, f)
+    bias_eta_list.append(bias_eta)
+    bias_eta_err_list.append(bias_eta_err)
+    beta_list.append(beta)
+    beta_err_list.append(beta_err)
+    bias_list.append(bias)
+    bias_err_list.append(bias_err)
+    beff_list.append(beff)
+    beff_err_list.append(beff_err)
+
+print("\nbias_eta:")
+print(bias_eta_list)
+print(bias_eta_err_list)
+print("\beta:")
+print(beta_list)
+print(beta_err_list)
+print("\nbias:")
+print(bias_list)
+print(bias_err_list)
+print("\beff:")
+print(beff_list)
+print(beff_err_list)
+
+f1, ax1 = plt.subplots()
+ax1.errorbar(rmin_list, bias_eta_list, yerr=bias_eta_err_list, fmt='o')
+ax1.set_xlabel(r'$r_{min}$ [Mpc/h]')
+ax1.set_ylabel(r'$b_{\eta}$')
+ax1.grid()
+ax1.set_title(title)
+plt.tight_layout()
+
+f2, ax2 = plt.subplots()
+ax2.errorbar(rmin_list, beta_list, yerr=beta_err_list, fmt='o')
+ax2.set_xlabel(r'$r_{min}$ [Mpc/h]')
+ax2.set_ylabel(r'$\beta$')
+ax2.grid()
+ax2.set_title(title)
+plt.tight_layout()
+
+f3, ax3 = plt.subplots()
+ax3.errorbar(rmin_list, bias_list, yerr=bias_err_list, fmt='o')
+ax3.set_xlabel(r'$r_{min}$ [Mpc/h]')
+ax3.set_ylabel(r'$b$')
+ax3.grid()
+ax3.set_title(title)
+plt.tight_layout()
+
+f4, ax4 = plt.subplots()
+ax4.errorbar(rmin_list, beff_list, yerr=beff_err_list, fmt='o')
+ax4.set_xlabel(r'$r_{min}$ [Mpc/h]')
+ax4.set_ylabel(r'$b_{eff}$')
+ax4.grid()
+ax4.set_title(title)
+plt.tight_layout()
+
+plt.show()

analysis/compare_cf.py

@@ -1,73 +1,134 @@
+import fitsio
+import picca.wedgize
+import argparse
+import matplotlib.pyplot as plt
+import numpy as np
+import h5py
 
 
-r_pow = 2
-rtmax = 200
-rpmax = 200
+# Plot options
+SMALL_SIZE = 18
+MEDIUM_SIZE = 20
+BIGGER_SIZE = 22
+plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
+plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
+plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
+plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
 
-f, ax = plt.subplots(figsize=(12,8))
+plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
+plt.rc('figure', figsize=(9,7))
 
-data = f1
-da = data[1]['DA'][:]
-co = data[1]['CO'][:]
-mu0, mu1, mu2, mu3 = 0, 0.5, 0.8, 1
-w1 = picca.wedgize.wedge(mumin=mu0,mumax=mu1, rtmax=rtmax, rpmax=rpmax)
-w2 = picca.wedgize.wedge(mumin=mu1,mumax=mu2, rtmax=rtmax, rpmax=rpmax)
-w3 = picca.wedgize.wedge(mumin=mu2,mumax=mu3, rtmax=rtmax, rpmax=rpmax)
-data_wedge1 = w1.wedge(da,co)
-coef1 = data_wedge1[0]**r_pow
-data_wedge2 = w2.wedge(da,co)
-coef2 = data_wedge2[0]**r_pow
-data_wedge3 = w3.wedge(da,co)
-coef3 = data_wedge3[0]**r_pow
+parser = argparse.ArgumentParser()
+parser.add_argument("-i", type=str, nargs="*", help="list of files")
+parser.add_argument("--labels", type=str, nargs='*', help='list of labels', default=None)
+parser.add_argument("--rt-min", type=float, default=0.)
+parser.add_argument("--rt-max", type=float, default=200.)
+parser.add_argument("--nrt", type=int, default=50)
+parser.add_argument("--rp-min", type=float, default=0.)
+parser.add_argument("--rp-max", type=float, default=200.)
+parser.add_argument("--nrp", type=int, default=50)
+parser.add_argument("--mu-min", type=float, default=-1.)
+parser.add_argument("--mu-max", type=float, default=1.)
+parser.add_argument("--r-pow", type=int, default=2)
+parser.add_argument("--error-bar", action='store_true')
+parser.add_argument("--title", default="")
 
-ax.errorbar(data_wedge1[0],coef1*data_wedge1[1],yerr=coef1*sp.sqrt(sp.diag(data_wedge1[2])),marker='.', label="1", color='r')
-ax.errorbar(data_wedge2[0],coef2*data_wedge2[1],yerr=coef2*sp.sqrt(sp.diag(data_wedge2[2])),marker='.', color='r')
-ax.errorbar(data_wedge3[0],coef3*data_wedge3[1],yerr=coef3*sp.sqrt(sp.diag(data_wedge3[2])),marker='.', color='r')
+args = parser.parse_args()
 
-data = f2
-da = data[1]['DA'][:]
-co = data[1]['CO'][:]
-mu0, mu1, mu2, mu3 = 0, 0.5, 0.8, 1
-w1 = picca.wedgize.wedge(mumin=mu0,mumax=mu1, rtmax=rtmax, rpmax=rpmax)
-w2 = picca.wedgize.wedge(mumin=mu1,mumax=mu2, rtmax=rtmax, rpmax=rpmax)
-w3 = picca.wedgize.wedge(mumin=mu2,mumax=mu3, rtmax=rtmax, rpmax=rpmax)
-data_wedge1 = w1.wedge(da,co)
-coef1 = data_wedge1[0]**r_pow
-data_wedge2 = w2.wedge(da,co)
-coef2 = data_wedge2[0]**r_pow
-data_wedge3 = w3.wedge(da,co)
-coef3 = data_wedge3[0]**r_pow
+files = args.i
+labels = args.labels
+rpmin = args.rp_min
+rpmax = args.rp_max
+nrp = args.nrp
+rtmin = args.rt_min
+rtmax = args.rt_max
+nrt = args.nrt
+mumin = args.mu_min
+mumax = args.mu_max
+r_pow = args.r_pow
 
-ax.errorbar(data_wedge1[0],coef1*data_wedge1[1],yerr=coef1*sp.sqrt(sp.diag(data_wedge1[2])),marker='.', label="2", color='b')
-ax.errorbar(data_wedge2[0],coef2*data_wedge2[1],yerr=coef2*sp.sqrt(sp.diag(data_wedge2[2])),marker='.', color='b')
-ax.errorbar(data_wedge3[0],coef3*data_wedge3[1],yerr=coef3*sp.sqrt(sp.diag(data_wedge3[2])),marker='.', color='b')
+if labels is None:
+    labels = np.arange(len(files))
 
-data = f3
-da = data[1]['DA'][:]
-co = data[1]['CO'][:]
-mu0, mu1, mu2, mu3 = 0, 0.5, 0.8, 1
-w1 = picca.wedgize.wedge(mumin=mu0,mumax=mu1, rtmax=rtmax, rpmax=rpmax)
-w2 = picca.wedgize.wedge(mumin=mu1,mumax=mu2, rtmax=rtmax, rpmax=rpmax)
-w3 = picca.wedgize.wedge(mumin=mu2,mumax=mu3, rtmax=rtmax, rpmax=rpmax)
-data_wedge1 = w1.wedge(da,co)
-coef1 = data_wedge1[0]**r_pow
-data_wedge2 = w2.wedge(da,co)
-coef2 = data_wedge2[0]**r_pow
-data_wedge3 = w3.wedge(da,co)
-coef3 = data_wedge3[0]**r_pow
+# fmt = ['.', '.', '.', 'x', '+', 'o', '.', 'x', 'o', '-.', ':']
+fmt = ['.', '.', '.', '.', '.', '.', '.', '.', '.', '.']
+# colors = ['black', 'darkblue', 'darkgreen', 'red', 'darkorange', 'darkviolet', 'saddlebrown', 'dodgerblue', 'deeppink']
+# colors = ['b', 'darkorange', 'r']
+# colors = ['black', 'royalblue', 'r']
+colors = ['darkblue', 'darkblue', 'r', 'r']
+linestyles = ['-', '--', '-', '--']
 
-ax.errorbar(data_wedge1[0],coef1*data_wedge1[1],yerr=coef1*sp.sqrt(sp.diag(data_wedge1[2])),marker='.', label="3", color='g')
-ax.errorbar(data_wedge2[0],coef2*data_wedge2[1],yerr=coef2*sp.sqrt(sp.diag(data_wedge2[2])),marker='.', color='g')
-ax.errorbar(data_wedge3[0],coef3*data_wedge3[1],yerr=coef3*sp.sqrt(sp.diag(data_wedge3[2])),marker='.', color='g')
+mu0, mu1, mu2, mu3, mu4 = 0, 0.5, 0.8, 0.95, 1
+w = picca.wedgize.wedge(mumin=mumin,mumax=mumax, rtmax=rtmax, rpmax=rpmax, rtmin=rtmin, rpmin=rpmin, nrt=nrt, nrp=nrp,absoluteMu=True)
+w1 = picca.wedgize.wedge(mumin=mu0,mumax=mu1, rtmax=rtmax, rpmax=rpmax, rtmin=rtmin, rpmin=rpmin, nrt=nrt, nrp=nrp,absoluteMu=True)
+w2 = picca.wedgize.wedge(mumin=mu1,mumax=mu2, rtmax=rtmax, rpmax=rpmax, rtmin=rtmin, rpmin=rpmin, nrt=nrt, nrp=nrp,absoluteMu=True)
+w3 = picca.wedgize.wedge(mumin=mu2,mumax=mu3, rtmax=rtmax, rpmax=rpmax, rtmin=rtmin, rpmin=rpmin, nrt=nrt, nrp=nrp,absoluteMu=True)
+w4 = picca.wedgize.wedge(mumin=mu3,mumax=mu4, rtmax=rtmax, rpmax=rpmax, rtmin=rtmin, rpmin=rpmin, nrt=nrt, nrp=nrp,absoluteMu=True)
 
-ax.set_xlabel(r"$r \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
-if r_pow == 2:
-    ax.set_ylabel(r"$r^{2}\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
-if r_pow == 1:
-    ax.set_ylabel(r"$r\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
-if r_pow == 0:
-    ax.set_ylabel(r"$\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
-ax.grid()
-ax.legend()
-plt.show()
+f1, ax1 = plt.subplots(figsize=(12,8))
+f2, ax2 = plt.subplots(figsize=(12,8))
+for i, f in enumerate(files):
+    if '.fits' in f:
+        h5_flag = False
+        data = fitsio.FITS(f)
+        da = data[1]['DA'][:]
+        co = data[1]['CO'][:]
+    if '.h5' in f:
+        print("Reading {} file ; using cov mat from {}".format(f, files[i-1]))
+        h5_flag = True
+        ff = h5py.File(f, 'r')
+        if "LYA(LYA)xLYA(LYA)" in ff.keys():
+            da = ff["LYA(LYA)xLYA(LYA)/fit"][...]
+        elif "cf_z_0_10" in ff.keys():
+            da = ff["cf_z_0_10/fit"][...]
+        else:
+            idx1 = int(f.rfind("/"))+1
+            idx2 = int(f.find(".h5"))
+            da = ff[f[idx2:idx2]+"/fit"][...]
+        ff.close()
+
+    data_wedge = w.wedge(da,co)
+    coef = data_wedge[0]**r_pow
+    data_wedge1 = w1.wedge(da,co)
+    coef1 = data_wedge1[0]**r_pow
+    data_wedge2 = w2.wedge(da,co)
+    coef2 = data_wedge2[0]**r_pow
+    data_wedge3 = w3.wedge(da,co)
+    coef3 = data_wedge3[0]**r_pow
+    data_wedge4 = w4.wedge(da,co)
+    coef4 = data_wedge4[0]**r_pow
+    if h5_flag:
+        ax1.plot(data_wedge[0],coef*data_wedge[1], label=labels[i], color=colors[i], linestyle=linestyles[i])
+        ax2.plot(data_wedge1[0],coef1*data_wedge1[1], label=labels[i], color=colors[i], linestyle=linestyles[i])
+        ax2.plot(data_wedge2[0],coef2*data_wedge2[1], color=colors[i], linestyle=linestyles[i])
+        ax2.plot(data_wedge3[0],coef3*data_wedge3[1], color=colors[i], linestyle=linestyles[i])
+        ax2.plot(data_wedge4[0],coef4*data_wedge4[1], color=colors[i], linestyle=linestyles[i])
+    elif args.error_bar:
+        ax1.errorbar(data_wedge[0],coef*data_wedge[1],yerr=coef*np.sqrt(np.diag(data_wedge[2])), label=labels[i], color=colors[i], fmt=fmt[i])
+        ax2.errorbar(data_wedge1[0],coef1*data_wedge1[1],yerr=coef1*np.sqrt(np.diag(data_wedge1[2])), label=labels[i], color=colors[i], fmt=fmt[i])
+        ax2.errorbar(data_wedge2[0],coef2*data_wedge2[1],yerr=coef2*np.sqrt(np.diag(data_wedge2[2])), color=colors[i], fmt=fmt[i])
+        ax2.errorbar(data_wedge3[0],coef3*data_wedge3[1],yerr=coef3*np.sqrt(np.diag(data_wedge3[2])), color=colors[i], fmt=fmt[i])
+        ax2.errorbar(data_wedge4[0],coef4*data_wedge4[1],yerr=coef4*np.sqrt(np.diag(data_wedge4[2])), color=colors[i], fmt=fmt[i])
+    else:
+        ax1.plot(data_wedge[0],coef*data_wedge[1], label=labels[i], color=colors[i], marker=fmt[i])
+        ax2.plot(data_wedge1[0],coef1*data_wedge1[1], label=labels[i], color=colors[i], marker=fmt[i])
+        ax2.plot(data_wedge2[0],coef2*data_wedge2[1], color=colors[i], marker=fmt[i])
+        ax2.plot(data_wedge3[0],coef3*data_wedge3[1], color=colors[i], marker=fmt[i])
+        ax2.plot(data_wedge4[0],coef4*data_wedge4[1], color=colors[i], marker=fmt[i])
 
+for ax in [ax1, ax2]:
+    ax.grid()
+    ax.legend()
+    ax.set_title("CF - {}".format(args.title), fontsize=20)
+    ax.set_xlabel(r"$r \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
+    if r_pow == 2:
+        ax.set_ylabel(r"$r^{2}\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
+    if r_pow == 1:
+        ax.set_ylabel(r"$r\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
+    if r_pow == 0:
+        ax.set_ylabel(r"$\xi(r) \, [\mathrm{Mpc \, h^{-1}}]$",fontsize=20)
+    plt.tight_layout()
+
+plt.show()

analysis/dla_to_drq.py

@@ -2,39 +2,41 @@
 import argparse
 import numpy as np
 
+"""
+This codes convert the DLA catalog produced by the SaclayMocks to a 
+DLA catalog readable by picca (to compute DLA auto-correlation, or
+to remove DLA from delta files for instance
+"""
 
 parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
-parser.add_argument("--dla-cat", type=str, default=None, required=True, help='path to dla catalog file')
-parser.add_argument("--master", type=str, default=None, required=True, help='path to master file')
-parser.add_argument("--out", type=str, default=None, required=True, help='output drq')
+parser.add_argument("-i", type=str, default=None, required=True, help='path to input dla catalog file')
+parser.add_argument("-o", type=str, default=None, required=True, help='path to output dla catalog file')
+# parser.add_argument("--master", type=str, default=None, required=True, help='path to master file')
+parser.add_argument("--nhi-min", type=float, default=None, required=False, help='remove DLA with log(n_HI) < NHI_MIN')
 args = parser.parse_args()
 
 ### Data
-print("Reading {}".format(args.dla_cat))
-h = fitsio.FITS(args.dla_cat)
-master = fitsio.read(args.master, ext=1)
-data = {}
-data['Z'] = h[1]['Z_DLA'][:]
-data['THING_ID'] = h[1]['MOCKID'][:]
-data['PLATE'] = h[1]['MOCKID'][:]
-data['MJD'] = h[1]['MOCKID'][:]
-data['FIBERID'] = h[1]['MOCKID'][:]
-qso_id = master['THING_ID']
-data['RA'] = []
-data['DEC'] = []
-for i in data['THING_ID']:
-    msk = (i == qso_id)
-    data['RA'].append(master['RA'][msk])
-    data['DEC'].append(master['DEC'][msk])
-data['RA'] = np.array(data['RA'])
-data['DEC'] = np.array(data['DEC'])
-h.close()
+print("Reading {}".format(args.i))
+data = fitsio.read(args.i, ext='DLACAT')
+if args.nhi_min is not None:
+    msk = data['N_HI_DLA'] > args.nhi_min
+    print("Removing {} DLAs with log(n_HI) < {} out of {} total DLAs".format((~msk).sum(), args.nhi_min, data.size))
+else:
+    msk = np.bool_(np.ones(data.size))
 
-### Save data
-print(data['RA'].size)
-print("Writting out {}".format(args.out))
-out = fitsio.FITS(args.out,'rw',clobber=True)
-cols = [ v for k,v in data.items() ]
-names = [ k for k in data.keys() ]
-out.write(cols,names=names)
-out.close()
+outfits = fitsio.FITS(args.o, 'rw', clobber=True)
+names = ['RA', 'DEC', 'THING_ID', 'Z', 'PLATE', 'MJD', 'FIBERID', 'NHI', 'ZQSO']
+table = []
+table.append(data['RA'][msk])
+table.append(data['DEC'][msk])
+table.append(data['MOCKID'][msk])
+table.append(data['Z_DLA_RSD'][msk])
+table.append(data['MOCKID'][msk])
+table.append(data['MOCKID'][msk])
+table.append(data['MOCKID'][msk])
+table.append(data['N_HI_DLA'][msk])
+table.append(data['Z_QSO_RSD'][msk])
+
+print("Writting out {}".format(args.o))
+outfits.write(table, names=names, extname='DLACAT')
+outfits.close()