Sandydev: fixing velocity bias along x and y direction

abacusnbody/hod/GRAND_HOD.py

@@ -227,11 +227,11 @@ def gen_cent(pos, vel, mass, ids, multis, randoms, vdev, deltac, fenv, shear,
             if keep[i] == 1:
                 # loop thru three directions to assign galaxy velocities and positions
                 lrg_x[j1] = pos[i,0]
-                lrg_vx[j1] = vel[i,0] + alpha_c_L * vdev[i] # velocity bias
+                lrg_vx[j1] = vel[i,0] + alpha_c_L * vdev[i,0] # velocity bias
                 lrg_y[j1] = pos[i,1]
-                lrg_vy[j1] = vel[i,1] + alpha_c_L * vdev[i] # velocity bias
+                lrg_vy[j1] = vel[i,1] + alpha_c_L * vdev[i,1] # velocity bias
                 lrg_z[j1] = pos[i,2]
-                lrg_vz[j1] = vel[i,2] + alpha_c_L * vdev[i] # velocity bias
+                lrg_vz[j1] = vel[i,2] + alpha_c_L * vdev[i,2] # velocity bias
                 # rsd only applies to the z direction
                 if rsd and origin is not None:
                     nx = lrg_x[j1] - origin[0]
@@ -253,11 +253,11 @@ def gen_cent(pos, vel, mass, ids, multis, randoms, vdev, deltac, fenv, shear,
             elif keep[i] == 2:
                 # loop thru three directions to assign galaxy velocities and positions
                 elg_x[j2] = pos[i,0]
-                elg_vx[j2] = vel[i,0] + alpha_c_E * vdev[i] # velocity bias
+                elg_vx[j2] = vel[i,0] + alpha_c_E * vdev[i,0] # velocity bias
                 elg_y[j2] = pos[i,1]
-                elg_vy[j2] = vel[i,1] + alpha_c_E * vdev[i] # velocity bias
+                elg_vy[j2] = vel[i,1] + alpha_c_E * vdev[i,1] # velocity bias
                 elg_z[j2] = pos[i,2]
-                elg_vz[j2] = vel[i,2] + alpha_c_E * vdev[i] # velocity bias
+                elg_vz[j2] = vel[i,2] + alpha_c_E * vdev[i,2] # velocity bias
                 # rsd only applies to the z direction
                 if rsd and origin is not None:
                     nx = elg_x[j2] - origin[0]
@@ -279,11 +279,11 @@ def gen_cent(pos, vel, mass, ids, multis, randoms, vdev, deltac, fenv, shear,
             elif keep[i] == 3:
                 # loop thru three directions to assign galaxy velocities and positions
                 qso_x[j3] = pos[i,0]
-                qso_vx[j3] = vel[i,0] + alpha_c_Q * vdev[i] # velocity bias
+                qso_vx[j3] = vel[i,0] + alpha_c_Q * vdev[i,0] # velocity bias
                 qso_y[j3] = pos[i,1]
-                qso_vy[j3] = vel[i,1] + alpha_c_Q * vdev[i] # velocity bias
+                qso_vy[j3] = vel[i,1] + alpha_c_Q * vdev[i,1] # velocity bias
                 qso_z[j3] = pos[i,2]
-                qso_vz[j3] = vel[i,2] + alpha_c_Q * vdev[i] # velocity bias
+                qso_vz[j3] = vel[i,2] + alpha_c_Q * vdev[i,2] # velocity bias
                 # rsd only applies to the z direction
                 if rsd and origin is not None:
                     nx = qso_x[j3] - origin[0]

abacusnbody/hod/abacus_hod.py

@@ -123,6 +123,7 @@ def __init__(self, sim_params, HOD_params, clustering_params = None, chunk=-1, n
         self.want_ranks = HOD_params.get('want_ranks', False)
         self.want_AB = HOD_params.get('want_AB', False)
         self.want_shear = HOD_params.get('want_shear', False)
+        self.want_expvel = HOD_params.get('want_expvel', False)
         self.want_rsd = HOD_params['want_rsd']
 
         if clustering_params is not None:
@@ -250,7 +251,7 @@ def staging(self):
         hid = np.empty([Nhalos_tot], dtype = int)
         hmultis = np.empty([Nhalos_tot])
         hrandoms = np.empty([Nhalos_tot])
-        hveldev = np.empty([Nhalos_tot])
+        hveldev = np.empty((Nhalos_tot, 3))
         hsigma3d = np.empty([Nhalos_tot])
         hc = np.empty([Nhalos_tot])
         hrvir = np.empty([Nhalos_tot])
@@ -308,7 +309,10 @@ def staging(self):
             halo_ids = maskedhalos["id"].astype(int) # halo IDs
             halo_pos = maskedhalos["x_L2com"] # halo positions, Mpc / h
             halo_vels = maskedhalos['v_L2com'] # halo velocities, km/s
-            halo_vel_dev = maskedhalos["randoms_gaus_vrms"] # halo velocity dispersions, km/s
+            if self.want_expvel:
+                halo_vel_dev = maskedhalos["randoms_exp"] # halo velocity dispersions, km/s
+            else:
+                halo_vel_dev = maskedhalos["randoms_gaus_vrms"] # halo velocity dispersions, km/s
             halo_sigma3d = maskedhalos["sigmav3d_L2com"] # 3d velocity dispersion
             halo_c = maskedhalos['r98_L2com']/maskedhalos['r25_L2com'] # concentration
             halo_rvir = maskedhalos['r98_L2com'] # rvir but using r98
@@ -531,7 +535,19 @@ def run_hod(self, tracers = None, want_rsd = True, want_nfw = False, NFW_draw =
             # np.random.seed(reseed)
             mtg = MTGenerator(np.random.PCG64(reseed))
             r1 = mtg.random(size=len(self.halo_data['hrandoms']), nthread=Nthread, dtype=np.float32)
-            r2 = mtg.standard_normal(size=len(self.halo_data['hveldev']), nthread=Nthread, dtype=np.float32)
+            if self.want_expvel:
+                rt0 = mtg.random(size=len(self.halo_data['hrandoms']), nthread=Nthread, dtype=np.float32)
+                rt1 = mtg.random(size=len(self.halo_data['hrandoms']), nthread=Nthread, dtype=np.float32)
+                rt2 = mtg.random(size=len(self.halo_data['hrandoms']), nthread=Nthread, dtype=np.float32)
+                rt = np.vstack((rt0, rt1, rt2)).T
+                r2 = np.zeros((len(rt), 3), dtype=np.float32)
+                r2[rt >= 0.5] = -np.log(2*(1-rt[rt >= 0.5]))
+                r2[rt < 0.5] = np.log(2*rt[rt < 0.5])
+            else:
+                r20 = mtg.standard_normal(size=len(self.halo_data['hveldev']), nthread=Nthread, dtype=np.float32)
+                r21 = mtg.standard_normal(size=len(self.halo_data['hveldev']), nthread=Nthread, dtype=np.float32)
+                r22 = mtg.standard_normal(size=len(self.halo_data['hveldev']), nthread=Nthread, dtype=np.float32)
+                r2 = np.vstack((r20, r21, r22)).T
             r3 = mtg.random(size=len(self.particle_data['prandoms']), nthread=Nthread, dtype=np.float32)
             self.halo_data['hrandoms'] = r1
             self.halo_data['hveldev'] = r2*self.halo_data['hsigma3d']/np.sqrt(3)

abacusnbody/hod/prepare_sim.py

@@ -624,8 +624,10 @@ def prepare_slab(i, savedir, simdir, simname, z_mock, z_type, tracer_flags, MT,
     halos['npstartA'] = halos_pstart_new
     halos['npoutA'] = halos_pnum_new
     halos['randoms'] = np.random.random(len(halos)) # attaching random numbers
+    halos['randoms_exp'] = (np.random.randint(0, 2, size = (len(halos), 3))*2-1)\
+    *np.random.exponential(scale = np.repeat(halos["sigmav3d_L2com"], 3).reshape((-1, 3))/np.sqrt(3), size = (len(halos), 3)) # attaching random numbers
     halos['randoms_gaus_vrms'] = np.random.normal(loc = 0,
-        scale = halos["sigmav3d_L2com"]/np.sqrt(3), size = len(halos)) # attaching random numbers
+        scale = np.repeat(halos["sigmav3d_L2com"], 3).reshape((-1, 3))/np.sqrt(3), size = (len(halos), 3)) # attaching random numbers
 
     # output halo file
     print("outputting new halo file ")

tests/ref_hod/AbacusSummit_base_c000_ph001-abridged/z2.250/galaxies_rsd/ELGs.dat

@@ -20,4 +20,4 @@ x y z vx vy vz mass id
 -989.9779939780151 2595.872235514386 174.58946130616962 279.78515625 405.76171875 682.470703125 594760988767.7638 1364044
 -989.8491016637607 2700.924001600403 128.59961886324763 192.24267578125 -99.3597640991211 -54.88894271850586 4730669850376.221 1374058
 -989.9346596474115 2712.400836645169 79.9785797238952 346.61865234375 -133.11767578125 -166.9921875 491415994265.56366 1375036
--989.411052509211 2702.3862712400905 129.18043285068367 17.578125 143.5 -284.125 2389589362673.3203 1374059
+-989.8232582237715 2546.1908703460344 572.2576660921198 618.0 278.25 -720.5 15543930805739.074 1360581