first push

lya_2pt/correlation.py

@@ -6,6 +6,21 @@
 from lya_2pt.tracer_utils import get_angle
 
 
+import pdb
+class ForkedPdb(pdb.Pdb):
+    """A Pdb subclass that may be used
+    from a forked multiprocessing child
+
+    """
+    def interaction(self, *args, **kwargs):
+        _stdin = sys.stdin
+        try:
+            sys.stdin = open('/dev/stdin')
+            pdb.Pdb.interaction(self, *args, **kwargs)
+        finally:
+            sys.stdin = _stdin
+
+
 def compute_xi(healpix_id):
     """Compute correlation function
 
@@ -38,6 +53,7 @@ def compute_xi(healpix_id):
     rt_grid = np.zeros(total_size)
     z_grid = np.zeros(total_size)
     num_pairs_grid = np.zeros(total_size, dtype=np.int32)
+    gamma_grid = np.zeros(total_size)
 
     for tracer1 in globals.tracers1[healpix_id]:
         with globals.lock:
@@ -55,6 +71,9 @@ def compute_xi(healpix_id):
             globals.rp_max, globals.rt_max
             )
 
+        ForkedPdb().set_trace()
+
+
         for tracer2 in neighbours:
             angle = get_angle(
                 tracer1.x_cart, tracer1.y_cart, tracer1.z_cart, tracer1.ra, tracer1.dec,
@@ -73,8 +92,9 @@ def compute_xi(healpix_id):
     rp_grid[w] /= weights_grid[w]
     rt_grid[w] /= weights_grid[w]
     z_grid[w] /= weights_grid[w]
+    gamma_grid[w] /= weights_grid[w]
 
-    return healpix_id, (xi_grid, weights_grid, rp_grid, rt_grid, z_grid, num_pairs_grid)
+    return healpix_id, (xi_grid, weights_grid, rp_grid, rt_grid, z_grid, num_pairs_grid, gamma_grid)
 
 
 @njit
@@ -117,3 +137,8 @@ def compute_xi_pair(
             rt_grid[bins] += rt * weight12
             z_grid[bins] += (z1[i] + z2[j]) / 2 * weight12
             num_pairs_grid[bins] += 1
+
+            #gamma grid
+            #globals.gamma_fun
+
+            #gamma_grid[bins] += globals.gamma_fun[i] * globals.gamma_fun[j] * weight12     #GAMMA FUNCTION

lya_2pt/global_data.py

@@ -22,3 +22,5 @@
 
 get_old_distortion = None
 rejection_fraction = None
+
+gamma_fun = None

lya_2pt/interface.py

@@ -95,6 +95,7 @@ def __init__(self, config):
         globals.num_bins_rt_model = self.settings.getint('num_bins_rt_model')
         globals.rejection_fraction = self.settings.getfloat('rejection_fraction')
         globals.get_old_distortion = self.settings.getboolean('get-old-distortion')
+        globals.gamma_fun = self.settings.getfloat("gamma_z_error")
 
         self.nside = self.settings.getint("nside")
         self.num_cpu = self.settings.getint("num-cpu")

lya_2pt/output.py

@@ -97,9 +97,9 @@ def write_cf_healpix(self, output, healpix_id, global_config, settings):
             correlation_name += "_BLIND"
 
         results.write(
-            [output[0], output[1]],
-            names=[correlation_name, "WEIGHT_SUM"],
-            comment=['unnormalized correlation', 'Sum of weight'],
+            [output[0], output[1], output[6]],
+            names=[correlation_name, "WEIGHT_SUM", 'GAMMA_GAMMA'],
+            comment=['unnormalized correlation', 'Sum of weight','gamma autocorrelation'],
             header=header2,
             extname='CORRELATION'
         )

lya_2pt/utils.py

@@ -1,6 +1,7 @@
 import os
 from os import mkdir
 from pathlib import Path
+from astropy.table import Table
 
 import numpy as np
 
@@ -142,3 +143,44 @@ def check_dir(dir: Path):
     """
     if not dir.is_dir():
         mkdir(dir)
+
+
+def line_prof(A,mu,sig,wave):
+    return A*(1/(2*np.sqrt(np.pi)*sig))*np.exp(-0.5*(wave-mu)**2/sig**2)
+
+def gen_cont(lrest,dv=1):
+    #tuning the amplitudes of peaks by fitting mocks with 250km/s error
+    amps=[30,1.5,1.5,0.5,1.5,1,1.5,5,25]
+    #emission line means
+    bs=[1025.7,1063,1073,1082,1084,1118,1128,1175,1215.6]
+    #tuning emission line default widths
+    cs=[10,5.5,3.5,5,5,4,4,7,15]
+
+    fdv = np.exp(dv/3e5)    
+    cs = [np.sqrt(c**2+(b*(fdv-1))**2) for b,c in zip(bs,cs)]
+    line_props = Table({'amp':amps,'lambda_line':bs,'width':cs})
+
+    #flux of smooth component
+    smooth_level = 1
+    scale_factor = 1 
+
+    #gaussian peaks of emission lines onto smooth components
+    continuum = smooth_level
+    #lyb
+    continuum += line_prof(*list(line_props)[0],lrest)
+    continuum += line_prof(*list(line_props)[1],lrest)
+    continuum += line_prof(*list(line_props)[2],lrest)
+    continuum += line_prof(*list(line_props)[3],lrest)
+    continuum += line_prof(*list(line_props)[4],lrest)
+    continuum += line_prof(*list(line_props)[5],lrest)
+    continuum += line_prof(*list(line_props)[6],lrest)
+    #CIII]
+    continuum += line_prof(*list(line_props)[7],lrest)
+    #lya
+    continuum += line_prof(*list(line_props)[8],lrest)
+
+    return continuum/scale_factor
+
+def gen_gamma(lrest,sigma_v):
+    gamma_fun = get_cont(lrest,sigma_v)/get_cont(lrest,0) - 1
+    return gamma_fun