Skip to content

Commit d748e2f

Browse files
echaussidonechaussidon
authored
Merge branch 'main' into blinding_edmond
2 parents ca06919 + a0ad0c2 commit d748e2f

File tree

3 files changed

+141
-2
lines changed

3 files changed

+141
-2
lines changed

‎py/LSS/cosmodesi_io_tools.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@
2727

2828

2929

30-
logger = logging.getLogger('xirunpc')
30+
logger = logging.getLogger('cosmodesi_io')
3131

3232

3333
def get_scratch_dir():

‎scripts/mock_tools/mkfast_Y1.py

+138
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
from astropy.io import fits # Access to FITS (Flexible Image Transport System) files.
2+
from astropy.table import Table, hstack, vstack, Column # A class to represent tables of heterogeneous data.
3+
import fitsio
4+
import numpy as np
5+
import glob
6+
import os
7+
import healpy as hp
8+
import argparse
9+
import sys
10+
11+
12+
def radec2thphi(ra,dec):
13+
return (-dec+90.)*np.pi/180.,ra*np.pi/180.
14+
15+
16+
if os.environ['NERSC_HOST'] == 'cori':
17+
scratch = 'CSCRATCH'
18+
elif os.environ['NERSC_HOST'] == 'perlmutter':
19+
scratch = 'PSCRATCH'
20+
else:
21+
print('NERSC_HOST is not cori or permutter but is '+os.environ['NERSC_HOST'])
22+
sys.exit('NERSC_HOST not known (code only works on NERSC), not proceeding')
23+
24+
25+
parser = argparse.ArgumentParser()
26+
parser.add_argument("--tracer", help="tracer type to use",choices=['LRG','ELG_LOPnotqso','QSO','BGS_BRIGHT'],default='LRG')
27+
parser.add_argument("--mocktype", help="type of mock to use",choices=['Ab','EZ','EZ2gpc'],default='EZ')
28+
parser.add_argument("--mockversion", help="version of mock type ",choices=['1stgen'],default='1stgen')
29+
#parser.add_argument("--mockpath", help="Location of mock file(s)",default='/global/cfs/cdirs/desi/cosmosim/FirstGenMocks/AbacusSummit/CutSky/')
30+
#parser.add_argument("--mockfile", help="formattable name of mock file(s). e.g. cutsky_{TYPE}_{Z}_AbacusSummit_base_c000_ph{PH}.fits. TYPE will be replaced with tracer type. PH will be replaced with realization number for simulation of mock.",default='cutsky_{TYPE}_{Z}_AbacusSummit_base_c000_ph{PH}.fits')
31+
parser.add_argument("--real", help="number for the realization",default=1,type=int)
32+
#parser.add_argument("--survey", help="points to set of tiles",default='Y1')
33+
parser.add_argument("--reg", help="region",choices=['NGC','SGC'],default='NGC')
34+
parser.add_argument("--dataver", help="points to set of tiles",default='v0.1')
35+
parser.add_argument("--fastver", help="version for output",default='test')
36+
parser.add_argument("--base_output", help="base directory for output",default='/global/cfs/cdirs/desi/survey/catalogs/Y1/mocks/fast/')
37+
38+
args = parser.parse_args()
39+
40+
if args.tracer[:3] == 'BGS':
41+
prog = 'bright'
42+
else:
43+
prog = 'dark'
44+
#select mock data
45+
46+
outdir = args.base_output+args.mockversion+'/'+args.mocktype+'/'+args.fastver+'/'+args.tracer+'/'
47+
if not os.path.exists(outdir):
48+
os.makedirs(outdir)
49+
foutname = outdir + args.tracer+'_'+str(args.real)+'_'+args.reg+'_clustering.dat.fits'
50+
51+
print('output will be written to '+foutname)
52+
53+
if args.mockversion == '1stgen':
54+
zs = {'ELG':'z1.100','LRG':'z0.800','QSO':'z1.400'}
55+
if args.mocktype == 'Ab':
56+
mockpath = '/global/cfs/cdirs/desi/cosmosim/FirstGenMocks/AbacusSummit/CutSky/'
57+
file_name = 'cutsky_'+args.tracer[:3]+'_'+zs[args.tracer[:3]]+'_AbacusSummit_base_c000_ph'+str(args.real)+'.fits'
58+
if args.mocktype == 'EZ':
59+
mockpath = '/global/cfs/cdirs/desi/cosmosim/FirstGenMocks/EZmock/CutSky_6Gpc/'+args.tracer[:3]+'/'+zs[args.tracer[:3]]+'/'
60+
if args.tracer == 'LRG':
61+
file_name = 'cutsky_LRG_z0.800_EZmock_B6000G1536Z0.8N216424548_b0.385d4r169c0.3_seed'+str(args.real)+'_'+args.reg+'.fits'
62+
if args.tracer[:3] == 'ELG':
63+
file_name = 'cutsky_ELG_z1.100_EZmock_B6000G1536Z1.1N648012690_b0.345d1.45r40c0.05_seed'+str(args.real)+'_'+args.reg+'.fits'
64+
if args.tracer == 'QSO':
65+
file_name = 'cutsky_QSO_z1.400_EZmock_B6000G1536Z1.4N27395172_b0.053d1.13r0c0.6_seed'+str(args.real)+'_'+args.reg+'.fits'
66+
def mask(main=0, nz=0, Y5=0, sv3=0):
67+
return main * (2**3) + sv3 * (2**2) + Y5 * (2**1) + nz * (2**0)
68+
thepath = os.path.join(mockpath,file_name)
69+
data = fitsio.read(thepath,columns=['RA','DEC','Z','Z_COSMO','STATUS'])#f[1].data
70+
status = data['STATUS'][()]
71+
idx = np.arange(len(status))
72+
mask_main = mask(main=0, nz=1, Y5=1, sv3=0)
73+
if args.tracer == 'LRG':
74+
mask_main = mask(main=1, nz=1, Y5=1, sv3=0)
75+
idx_main = idx[(status & (mask_main))==mask_main]
76+
data = data[idx_main]
77+
78+
#load geometric mask
79+
80+
healpix_mask = hp.read_map('/global/cfs/cdirs/desi/survey/catalogs/Y1/LSS/iron/LSScats/'+args.dataver+'/healpix_map_ran_comp_'+args.tracer+'.fits')
81+
82+
### SUBSAMPLE ###
83+
84+
th,phi = radec2thphi(data['RA'],data['DEC'])
85+
dpix = hp.ang2pix(1024,th,phi,nest=True)
86+
mask_comp = np.zeros(len(data))
87+
mask_comp = healpix_mask[dpix]
88+
rans = np.random.random(len(data))
89+
mask_keep = (rans < mask_comp)
90+
91+
print('mask will keep '+str(np.sum(mask_keep))+' for Y1 out of '+str(len(data)))
92+
93+
### ADD LINES TO GET RANDOMS AND SUBSAMPLE ###
94+
95+
### COMPLETENESS AS FUNCTION OF NTILE FROM https://desi.lbl.gov/trac/wiki/keyprojects/y1kp3/Y1details#Completenessfootprintstatistics
96+
97+
ntile_dic = {'BGS_BRIGHT':{0:0,1:0.512,2:755,3:0.889,4:0.954},'ELG_LOPnotqso':{0:0,1:0.308,2:0.400,3:0.528,4:0.658,5:0.767,6:0.853,7:0.917},\
98+
'LRG':{0:0,1:0.585,2:0.741,3:0.869,4:0.936,5:0.968,6:0.985,7:0.995},'QSO':{0:0,1:0.794,2:0.957,3:0.989,4:0.994,5:0.996,6:0.998,7:0.999}}
99+
100+
#load ntile map
101+
102+
ntile_map = hp.read_map('/global/cfs/cdirs/desi/survey/catalogs/Y1/LSS/iron/LSScats/'+args.dataver+'/healpix_map_ntile_'+prog+'.fits')
103+
104+
ntile = np.zeros(len(data))
105+
ntile = ntile_map[dpix]
106+
ntile_comp = np.zeros(len(data))
107+
for i in range(0,len(ntile)):
108+
nt = ntile[i]
109+
nti = int(nt)
110+
if nti == nt:
111+
cp = ntile_dic[args.tracer][nti]
112+
else:
113+
cp_low = ntile_dic[args.tracer][nti]
114+
cp_high = ntile_dic[args.tracer][nti+1]
115+
cp = cp_low + (nt-nti)*(cp_high-cp_low) #just linearly interpolate
116+
ntile_comp[i] = cp
117+
118+
wts = 1/ntile_comp
119+
120+
rans = np.random.random(len(data))
121+
comp_keep = (rans < ntile_comp)
122+
123+
print('we will keep '+str(np.sum(mask_keep&comp_keep))+' for Y1 out of '+str(len(data)))
124+
125+
data = Table(data)
126+
127+
data['WEIGHT'] = wts
128+
data = data[mask_keep&comp_keep]
129+
130+
### NEED TO ADD SOME KIND OF n(z) sub-sampling ###
131+
132+
133+
data.write(foutname,overwrite=True,format='fits')
134+
135+
136+
137+
138+
### ADD LINES TO SUBSAMPLE AS FUNCTION OF NTILE AND ADD WEIGHT THAT COMPENSATES

‎scripts/pkrun.py

+2-1
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,8 @@
2323
from pypower import CatalogFFTPower, PowerSpectrumStatistics, CatalogSmoothWindow, PowerSpectrumSmoothWindow, PowerSpectrumOddWideAngleMatrix, PowerSpectrumSmoothWindowMatrix, utils, setup_logging
2424
from LSS.tabulated_cosmo import TabulatedDESI
2525

26-
from xirunpc import read_clustering_positions_weights, concatenate_data_randoms, compute_angular_weights, catalog_dir, get_regions, get_zlims, get_scratch_dir
26+
from xirunpc import compute_angular_weights
27+
from LSS.cosmodesi_io_tools import read_clustering_positions_weights, concatenate_data_randoms, catalog_dir, get_regions, get_zlims, get_scratch_dir
2728

2829

2930
os.environ['OMP_NUM_THREADS'] = os.environ['NUMEXPR_MAX_THREADS'] = '1'

0 commit comments

Comments
 (0)