Code for running DES+SPT Year 3 6x2pt analysis

boltzmann/camb/camb_interface.py

@@ -93,7 +93,10 @@ def setup(options):
     config['want_zdrag'] = mode != MODE_BG
     config['want_zstar'] = config['want_zdrag']
 
-
+    more_config['want_chistar'] = options.get_bool(opt, 'want_chistar', default=True)
+    more_config['n_logz'] = options.get_int(opt, 'n_logz', default=0)
+    more_config['zmax_logz'] = options.get_double(opt, 'zmax_logz', default = 1100.)
+
     more_config["lmax_params"] = get_optional_params(options, opt, ["max_eta_k", "lens_potential_accuracy",
                                                                     "lens_margin", "k_eta_fac", "lens_k_eta_reference",
                                                                     #"min_l", "max_l_tensor", "Log_lvalues", , "max_eta_k_tensor"
@@ -440,6 +443,10 @@ def save_distances(r, p, block, more_config):
     z_background = np.linspace(
         more_config["zmin_background"], more_config["zmax_background"], more_config["nz_background"])
 
+    #If desired, append logarithmically distributed redshifts
+    log_z = np.geomspace(more_config["zmax_background"], more_config['zmax_logz'], num = more_config['n_logz'])
+    z_background = np.append(z_background, log_z[1:])
+
     # Write basic distances and related quantities to datablock
     block[names.distances, "nz"] = len(z_background)
     block[names.distances, "z"] = z_background
@@ -481,6 +488,10 @@ def save_distances(r, p, block, more_config):
         block[names.distances, "rs_DV"] = rs_DV
         block[names.distances, "F_AP"] = F_AP
 
+    if more_config['want_chistar']:
+        chistar = (r.conformal_time(0)- r.tau_maxvis)
+        block[names.distances, "CHISTAR"] = chistar
+
 def compute_growth_factor(r, block, P_tot, k, z, more_config):
     if P_tot is None:
         # If we don't have it already, get the default matter power interpolator,
@@ -624,8 +635,7 @@ def sigma8_to_As(block, config, more_config, cosmology_params, dark_energy, reio
     r_temp = camb.get_results(p_temp)
     temp_sig8 = r_temp.get_sigma8()[-1] #what sigma8 comes out from using temp_As?
     As = temp_As*(sigma_8_input/temp_sig8)**2
-    block[cosmo,'A_s'] = As
-
+    block[cosmo,'A_s'] = As    
 
 def execute(block, config):
     config, more_config = config

cmb_lensing/kappa_beam/kappa_beam.py

@@ -10,47 +10,83 @@ def get_nbins(block, section):
         n_b=n_a
     return n_a, n_b
 
-def apply_beam(block, section, beam_sigma):
+def apply_beam(block, section, beam_sigma, is_auto=False, output_section = None, save_name = None):
     n_other, n_cmb = get_nbins(block, section)
     if (n_cmb != 1):
         raise ValueError("Found more than 1 CMB bin - there is definitely only 1 CMB.")
 
     ell = block[section, 'ell']
     B_ell = np.exp(-0.5*ell*(ell+1.)*beam_sigma**2.)
-    for bini in xrange(0,n_other):
+    for bini in range(0,n_other):
         C_ell_orig = block[section, 'bin_{}_1'.format(bini+1)]
-        beamed_C_ell = C_ell_orig*B_ell
-        block[section, 'bin_{}_1'.format(bini+1)] = beamed_C_ell
+        if (is_auto):
+            beamed_C_ell = C_ell_orig*(B_ell**2.)
+        else:
+            beamed_C_ell = C_ell_orig*B_ell
+        if (output_section is None):
+            block[section, 'bin_{}_1'.format(bini+1)] = beamed_C_ell
+        else:
+            block[output_section, 'bin_{}_1'.format(bini+1)] = beamed_C_ell
+
+    if  (section != output_section and (not output_section is None)):
+        if (block.has_value(section,"nbin")):
+            block[output_section, "nbin"] = block[section, "nbin"]
+        if (block.has_value(section,"nbin_a")):
+            block[output_section, "nbin_a"] = block[section, "nbin_a"]
+        if (block.has_value(section,"nbin_b")):
+            block[output_section, "nbin_b"] = block[section, "nbin_b"]
+        block[output_section, "save_name"] =  block[section, "save_name"]
+        block[output_section, "ell"] = block[section, "ell"]
+        block[output_section, "sep_name"] = block[section, "sep_name"]
+        block[output_section, "sample_a"] = block[section, "sample_a"]
+        block[output_section, "sample_b"] = block[section, "sample_b"]
+        block[output_section, "is_auto"] = block[section, "is_auto"]
+        if (save_name is not None):
+            block[output_section, "save_name"] = save_name
 
 def setup(options):
     shearkappa_section = options.get_string(option_section, "shearkappa_section", default="shear_cmbkappa_cl")
     galkappa_section = options.get_string(option_section, "galkappa_section", default="galaxy_cmbkappa_cl")
+    kappakappa_section = options.get_string(option_section, "kappakappa_section", default="cmbkappa_cl")
+    shearkappa_output_section = options.get_string(option_section, "shearkappa_output_section", default=None)
+    galkappa_output_section = options.get_string(option_section, "galkappa_output_section", default=None)
+    kappakappa_output_section = options.get_string(option_section, "kappakappa_output_section", default=None)
+    save_name = options.get_string(option_section, "save_name", default=None)
 
     if options.has_value(option_section, 'beam_sigma_arcmin'):
         beam_sigma_arcmin = options[option_section, "beam_sigma_arcmin"]
         # radians
-        beam_sigma = np.radians(beam_sigma_arcmin/60.)
-
+        beam_sigma = np.radians(beam_sigma_arcmin/60.)   
     elif options.has_value(option_section, 'beam_fwhm_arcmin'):
         beam_fwhm_arcmin = options[option_section, "beam_fwhm_arcmin"]
+        # radians
         beam_sigma = (1./np.sqrt(8.*np.log(2.)))*np.radians(beam_fwhm_arcmin/60.)
     else:
         raise ValueError("No beam_sigma_arcmin or beam_fwhm_armcin supplied!")
 
 
-    return {"shearkappa_section":shearkappa_section, "galkappa_section":galkappa_section, "beam_sigma":beam_sigma}
+    return {"shearkappa_section":shearkappa_section, "galkappa_section":galkappa_section, "kappakappa_section":kappakappa_section,\
+            "shearkappa_output_section":shearkappa_output_section, "galkappa_output_section":galkappa_output_section, "kappakappa_output_section":kappakappa_output_section,\
+            "beam_sigma":beam_sigma, "save_name":save_name}
 
 def execute(block, config):
     shearkappa_section = config['shearkappa_section']
     galkappa_section = config['galkappa_section']
+    kappakappa_section = config['kappakappa_section']
+    shearkappa_output_section = config['shearkappa_output_section']
+    galkappa_output_section = config['galkappa_output_section']
+    kappakappa_output_section = config['kappakappa_output_section']
     beam_sigma = config['beam_sigma']
+    save_name = config['save_name']
 
     if block.has_section(galkappa_section):
-        apply_beam(block, galkappa_section, beam_sigma)
-
+        apply_beam(block, galkappa_section, beam_sigma, is_auto = False, output_section  = galkappa_output_section, save_name = save_name)
 
     if block.has_section(shearkappa_section):
-        apply_beam(block, shearkappa_section, beam_sigma)
+        apply_beam(block, shearkappa_section, beam_sigma, is_auto = False, output_section = shearkappa_output_section, save_name = save_name)
+
+    if block.has_section(kappakappa_section):
+        apply_beam(block, kappakappa_section, beam_sigma, is_auto = True, output_section = kappakappa_output_section, save_name = save_name)
 
     return 0
 

cmb_lensing/kappa_ell_cut/kappa_ell_cut.py

@@ -12,7 +12,7 @@ def get_nbins(block, section):
     return n_a, n_b
 
 
-def apply_lcut(block, section, Lmin, Lmax):
+def apply_lcut(block, section, Lmin, Lmax, output_section = None, save_name = None):
     n_other, n_cmb = get_nbins(block, section)
     if (n_cmb != 1):
         raise ValueError("Found more than 1 CMB bin - there is definitely only 1 CMB.")
@@ -21,14 +21,42 @@ def apply_lcut(block, section, Lmin, Lmax):
     w   = np.ones(ell.shape[0])
     w[ell<Lmin]=0
     w[ell>Lmax+1]=0
+
     for bini in range(0,n_other):
         C_ell_orig = block[section, 'bin_{}_1'.format(bini+1)]
-        beamed_C_ell = C_ell_orig*w
-        block[section, 'bin_{}_1'.format(bini+1)] = beamed_C_ell
+        filtered_C_ell = C_ell_orig*w
+        if (not output_section is None):
+            block[output_section, 'bin_{}_1'.format(bini+1)] = filtered_C_ell
+        if (output_section is None):
+            block[section, 'bin_{}_1'.format(bini+1)] = filtered_C_ell
+
+    if  (section != output_section and (not output_section is None)):
+        if (block.has_value(section,"nbin")):
+            block[output_section, "nbin"] = block[section, "nbin"]
+        if (block.has_value(section,"nbin_a")):
+            block[output_section, "nbin_a"] = block[section, "nbin_a"]
+        if (block.has_value(section,"nbin_b")):
+            block[output_section, "nbin_b"] = block[section, "nbin_b"]
+        block[output_section, "save_name"] =  block[section, "save_name"]
+        block[output_section, "ell"] = block[section, "ell"]
+        block[output_section, "sep_name"] = block[section, "sep_name"]
+        block[output_section, "sample_a"] = block[section, "sample_a"]
+        block[output_section, "sample_b"] = block[section, "sample_b"]
+        block[output_section, "is_auto"] = block[section, "is_auto"]
+        if (save_name is not None):
+            block[output_section, "save_name"] = save_name
 
 def setup(options):
     shearkappa_section = options.get_string(option_section, "shearkappa_section", default="shear_cmbkappa_cl")
     galkappa_section = options.get_string(option_section, "galkappa_section", default="galaxy_cmbkappa_cl")
+    kappakappa_section = options.get_string(option_section, "kappakappa_section", default="cmbkappa_cl")
+    shearkappa_output_section = options.get_string(option_section, "shearkappa_output_section", default=None)
+    galkappa_output_section = options.get_string(option_section, "galkappa_output_section", default=None)
+    kappakappa_output_section = options.get_string(option_section, "kappakappa_output_section", default=None)
+    try:
+        save_name = options.get_string(option_section, "save_name", default=None)
+    except:
+        save_name = None
 
     if options.has_value(option_section,'Lmin'):
         Lmin = options[option_section, 'Lmin']
@@ -41,20 +69,30 @@ def setup(options):
         warnings.warn("No Lmax supplied! Setting Lmax=999999.")
         Lmax = 999999
 
-    return {"shearkappa_section":shearkappa_section, "galkappa_section":galkappa_section, "Lmin":Lmin, "Lmax":Lmax}
+    return {"shearkappa_section":shearkappa_section, "galkappa_section":galkappa_section, "kappakappa_section":kappakappa_section, \
+            "shearkappa_output_section":shearkappa_output_section, "galkappa_output_section":galkappa_output_section, "kappakappa_output_section":kappakappa_output_section,
+            "Lmin":Lmin, "Lmax":Lmax, "save_name":save_name}
 
 def execute(block, config):
     shearkappa_section = config['shearkappa_section']
     galkappa_section = config['galkappa_section']
+    kappakappa_section = config['kappakappa_section']
+    shearkappa_output_section = config['shearkappa_output_section']
+    galkappa_output_section = config['galkappa_output_section']
+    kappakappa_output_section = config['kappakappa_output_section']
+    save_name = config['save_name']
+
     Lmin = config['Lmin']
     Lmax = config['Lmax']
 
     if block.has_section(galkappa_section):
-        apply_lcut(block, galkappa_section, Lmin, Lmax)
-
+        apply_lcut(block, galkappa_section, Lmin, Lmax, output_section = galkappa_output_section, save_name = save_name)
 
     if block.has_section(shearkappa_section):
-        apply_lcut(block, shearkappa_section, Lmin, Lmax)
+        apply_lcut(block, shearkappa_section, Lmin, Lmax,  output_section = shearkappa_output_section, save_name = save_name)
+
+    if block.has_section(kappakappa_section):
+        apply_lcut(block, kappakappa_section, Lmin, Lmax, output_section = kappakappa_output_section, save_name = save_name)
 
     return 0
 

examples/des-y3-6x2-values.ini

@@ -0,0 +1,66 @@
+[cosmological_parameters]
+omega_m      =  0.1     0.3     0.9
+h0           =  0.55    0.69    0.91
+omega_b      =  0.03    0.048   0.07
+n_s          =  0.87    0.97    1.07
+A_s          =  0.5e-09 2.19e-9 5.0e-09
+w            =  -1.0
+
+; Old parameter names from original cosmosis version
+; omnuh2       =  0.0006  0.00083 0.00644
+; massive_nu   =  3
+; massless_nu  =  0.046
+
+; New parametrization and names:
+mnu = 0.06  0.0773   0.6
+num_massive_neutrinos = 3
+nnu = 3.046
+
+omega_k      =  0.0
+tau          =  0.0697186
+
+[shear_calibration_parameters]
+m1 = -0.1  0.0  0.1
+m2 = -0.1  0.0  0.1 
+m3 = -0.1  0.0  0.1
+m4 = -0.1  0.0  0.1
+
+[wl_photoz_errors]
+bias_1 = -0.1   0.0   0.1
+bias_2 = -0.1   0.0   0.1
+bias_3 = -0.1   0.0   0.1
+bias_4 = -0.1   0.0   0.1
+
+[lens_photoz_errors]
+bias_1 = -0.05  0.0  0.05
+bias_2 = -0.05  0.0  0.05
+bias_3 = -0.05  0.0  0.05
+bias_4 = -0.05  0.0  0.05
+bias_5 = -0.05  0.0  0.05
+width_1 = 1.0
+width_2 = 1.0
+width_3 = 1.0
+width_4 = 1.0
+width_5 = 0.1 1.0 1.9
+
+[bias_lens]
+b1 = 0.8  1.7   3.0
+b2 = 0.8  1.7   3.0
+b3 = 0.8  1.7   3.0
+b4 = 0.8  2.0   3.0
+b5 = 0.8  2.0   3.0
+
+[mag_alpha_lens]
+alpha_1 =  1.31
+alpha_2 = -0.52
+alpha_3 =  0.34
+alpha_4 =  2.25
+alpha_5 =  1.97
+
+[intrinsic_alignment_parameters]
+z_piv   =  0.62
+A1      = -5.0   0.7  5.0
+A2      = -5.0  -1.36 5.0 
+alpha1  = -5.0  -1.7  5.0
+alpha2  = -5.0  -2.5  5.0
+bias_ta =  0.0   1.0  2.0