diff --git a/post_process/py/extract_skewers.py b/post_process/py/extract_skewers.py index 6037c608..d924ff93 100644 --- a/post_process/py/extract_skewers.py +++ b/post_process/py/extract_skewers.py @@ -9,12 +9,14 @@ import camb_cosmo import thermal_model -def get_skewers_filename(num,n_skewers,width_Mpc,scale_T0=None, +def get_skewers_filename(num,n_skewers,width_Mpc,axis=None,scale_T0=None, scale_gamma=None): """Filename storing skewers for a particular temperature model""" filename='skewers_'+str(num)+'_Ns'+str(n_skewers) filename+='_wM'+str(int(1000*width_Mpc)/1000) + if axis is not None: + filename+='_axis'+str(axis) if scale_T0: filename+='_sT'+str(int(1000*scale_T0)/1000) if scale_gamma: @@ -23,11 +25,13 @@ def get_skewers_filename(num,n_skewers,width_Mpc,scale_T0=None, return filename -def get_snapshot_json_filename(num,n_skewers,width_Mpc): +def get_snapshot_json_filename(num,n_skewers,width_Mpc,axis=None): """Filename describing the set of skewers for a given snapshot""" filename='snap_skewers_'+str(num)+'_Ns'+str(n_skewers) filename+='_wM'+str(int(1000*width_Mpc)/1000) + if axis is not None: + filename+='_axis'+str(axis) filename+='.json' return filename @@ -55,8 +59,9 @@ def thermal_broadening_Mpc(T_0,dkms_dMpc): return sigma_T_Mpc -def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50, - width_Mpc=0.1,scales_T0=None,scales_gamma=None): +def rescale_write_skewers_z(simdir,num,skewers_dir=None, + n_skewers=50,width_Mpc=0.1,axis=None, + scales_T0=None,scales_gamma=None): """Extract skewers for a given snapshot, for different temperatures.""" # don't rescale unless asked to @@ -82,6 +87,8 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50, 'width_Mpc':width_Mpc, 'width_kms':width_kms, 'T0_ini':T0_ini, 'gamma_ini':gamma_ini, 'scales_T0':scales_T0, 'scales_gamma':scales_gamma} + if axis is not None: + sim_info['axis']=axis # will also store measured values sim_T0=[] @@ -96,19 +103,20 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50, for scale_gamma in scales_gamma: T0=T0_ini*scale_T0 gamma=gamma_ini*scale_gamma - sk_filename=get_skewers_filename(num,n_skewers,width_Mpc, - scale_T0,scale_gamma) + sk_filename=get_skewers_filename(num,n_skewers=n_skewers, + width_Mpc=width_Mpc,axis=axis, + scale_T0=scale_T0,scale_gamma=scale_gamma) # avoid (if possible) to use set_T0, might break fake_spectra if (scale_T0==1.0) and (scale_gamma==1.0): skewers=get_skewers_snapshot(simdir,skewers_dir,num, n_skewers=n_skewers,width_kms=width_kms, - skewers_filename=sk_filename) + axis=axis,skewers_filename=sk_filename) else: skewers=get_skewers_snapshot(simdir,skewers_dir,num, n_skewers=n_skewers,width_kms=width_kms, set_T0=T0,set_gamma=gamma, - skewers_filename=sk_filename) + axis=axis,skewers_filename=sk_filename) # call mean flux, so that the skewers are really computed mf=skewers.get_mean_flux() @@ -131,7 +139,7 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50, sim_info['sim_scale_gamma']=sim_scale_gamma sim_info['sk_files']=sk_files - snapshot_filename=get_snapshot_json_filename(num,n_skewers,width_Mpc) + snapshot_filename=get_snapshot_json_filename(num,n_skewers,width_Mpc,axis) sim_info['snapshot_filename']=snapshot_filename json_file = open(skewers_dir+'/'+snapshot_filename,"w") json.dump(sim_info,json_file) @@ -143,33 +151,33 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50, def get_skewers_snapshot(simdir,skewers_dir,snap_num,n_skewers=50,width_kms=10, - set_T0=None,set_gamma=None,skewers_filename=None): + axis=None,set_T0=None,set_gamma=None, + skewers_filename='test_skewers.hdf5'): """Extract skewers for a particular snapshot""" - if not skewers_filename: - skewers_filename="skewers_"+str(snap_num)+"_"+str(n_skewers) - skewers_filename+="_"+str(width_kms) - if set_T0: - skewers_filename+='_T0_'+str(set_T0) - if set_gamma: - skewers_filename+='_gamma_'+str(set_gamma) - skewers_filename+='.hdf5' - # check that spectra file does not exist yet (should crash) if os.path.exists(skewers_dir+'/'+skewers_filename): print(skewers_filename,'already exists in',skewers_dir) - # avoid (if possible) to use set_T0, not always present in fake_spectra - if (set_T0 is None) and (set_gamma is None): + # new functionality to specify simulation axis for skewers + if axis is not None: + # for now only working without temperature rescalings + assert (set_T0 is None) and (set_gamma is None), 'get_skewers_snapshot' skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/', - nspec=n_skewers,res=width_kms,savefile=skewers_filename, + nspec=n_skewers,res=width_kms,axis=axis, + savefile=skewers_filename, savedir=skewers_dir,reload_file=True) else: - skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/', - nspec=n_skewers,res=width_kms,savefile=skewers_filename, - savedir=skewers_dir,reload_file=True, - set_T0=set_T0,set_gamma=set_gamma) - + # avoid (if possible) to use set_T0, not always present in fake_spectra + if (set_T0 is None) and (set_gamma is None): + skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/', + nspec=n_skewers,res=width_kms,savefile=skewers_filename, + savedir=skewers_dir,reload_file=True) + else: + skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/', + nspec=n_skewers,res=width_kms,savefile=skewers_filename, + savedir=skewers_dir,reload_file=True, + set_T0=set_T0,set_gamma=set_gamma) return skewers diff --git a/post_process/py/snapshot_admin.py b/post_process/py/snapshot_admin.py index cf3809a8..2ed404cd 100644 --- a/post_process/py/snapshot_admin.py +++ b/post_process/py/snapshot_admin.py @@ -34,20 +34,20 @@ def __init__(self,snap_json,scales_tau=None,kF_Mpc=None): self.scales_tau=[1.0] - def get_all_flux_power(self,simdir=None): + def get_all_flux_power(self,simdir=None,skewers_dir=None): """Loop over all skewers, and return flux power for each""" - if simdir is not None: - # get box size from GenIC file, to normalize power - if 'simdir' in self.data: - simdir = self.data['simdir'] - elif 'basedir' in self.data: - simdir = self.data['basedir'] + if simdir is None: + simdir = self.data['simdir'] + if skewers_dir is None: + skewers_dir=self.data['skewers_dir'] + + # get box size from GenIC file, to normalize power genic_file=simdir+'/paramfile.genic' L_Mpc=read_genic.L_Mpc_from_paramfile(genic_file,verbose=True) - skewers_dir=simdir+"output/skewers/" + # get snapshot number from JSON data snap_num=self.data['snap_num'] # will loop over all temperature models in snapshot @@ -96,6 +96,10 @@ def get_p1d_json_filename(self,p1d_label): filename=p1d_label+'_'+str(num)+'_Ns'+str(n_skewers) filename+='_wM'+str(int(1000*width_Mpc)/1000) + + if 'axis' in self.data: + filename+='_axis'+str(self.data['axis']) + filename+='.json' return filename @@ -107,7 +111,8 @@ def write_p1d_json(self,p1d_label='p1d'): if p1d_label is None: p1d_label='p1d' - filename=self.data['simdir']+'/'+self.get_p1d_json_filename(p1d_label) + filename=self.data['skewers_dir']+'/'+self.get_p1d_json_filename(p1d_label) + print('will print P1d to file',filename) # make sure we have already computed P1D diff --git a/post_process/scripts/arxiv_flux_power.py b/post_process/scripts/arxiv_flux_power.py index a1472660..38a88543 100644 --- a/post_process/scripts/arxiv_flux_power.py +++ b/post_process/scripts/arxiv_flux_power.py @@ -14,9 +14,10 @@ parser = argparse.ArgumentParser() parser.add_argument('--simdir', type=str, help='Base simulation directory',required=True) parser.add_argument('--snap_num', type=int, help='Snapshop number',required=True) -parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=False) +parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=True) parser.add_argument('--n_skewers', type=int, default=10, help='Number of skewers per side',required=False) parser.add_argument('--width_Mpc', type=float, default=0.1, help='Cell width (in Mpc)',required=False) +parser.add_argument('--axis', type=int, help='Axis to use to extract skewers (1,2,3)', required=False) parser.add_argument('--scales_tau', type=str, default='1.0', help='Comma-separated list of optical depth scalings to use.',required=False) parser.add_argument('--p1d_label', type=str, default=None, help='String identifying P1D measurement and / or tau scaling.',required=False) parser.add_argument('--verbose', action='store_true', help='Print runtime information',required=False) @@ -25,19 +26,12 @@ verbose=args.verbose print('verbose =',verbose) -# main simulation folder -simdir=args.simdir -if args.skewers_dir: - skewers_dir=args.skewers_dir -else: - skewers_dir=simdir+'/output/skewers/' - scales_tau=[float(scale) for scale in args.scales_tau.split(',')] if verbose: print('will scale tau by',scales_tau) # try to read information about filtering length in simulation -kF_json=simdir+'/filtering_length.json' +kF_json=args.simdir+'/filtering_length.json' if os.path.isfile(kF_json): # read json file with filtering data with open(kF_json) as json_data: @@ -48,9 +42,9 @@ kF_Mpc=None # read file containing information of all temperature rescalings in snapshot -snap_filename=skewers_dir+'/'+extract_skewers.get_snapshot_json_filename( +snap_filename=args.skewers_dir+'/'+extract_skewers.get_snapshot_json_filename( num=args.snap_num,n_skewers=args.n_skewers, - width_Mpc=args.width_Mpc) + width_Mpc=args.width_Mpc,axis=args.axis) if verbose: print('setup snapshot admin from file',snap_filename) diff --git a/post_process/scripts/extract_default_skewers.py b/post_process/scripts/extract_default_skewers.py index 67ccdf01..8678e935 100644 --- a/post_process/scripts/extract_default_skewers.py +++ b/post_process/scripts/extract_default_skewers.py @@ -11,7 +11,8 @@ parser = argparse.ArgumentParser() parser.add_argument('--simdir', type=str, help='Base simulation directory',required=True) parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=True) -parser.add_argument('--snap_num', type=int, default=8, help='Snapshot number',required=False) +parser.add_argument('--snap_num', type=int, help='Snapshot number',required=True) +parser.add_argument('--axis', type=int, help='Axis to use to extract skewers (1,2,3)',required=False) parser.add_argument('--n_skewers', type=int, default=10, help='Number of skewers per side',required=False) parser.add_argument('--width_Mpc', type=float, default=0.1, help='Cell width (in Mpc)',required=False) parser.add_argument('--verbose', action='store_true', help='Print runtime information',required=False) @@ -24,7 +25,7 @@ # extract skewers for one snapshot, and return some info info=extract_skewers.rescale_write_skewers_z(simdir=args.simdir, num=args.snap_num, skewers_dir=args.skewers_dir, - n_skewers=args.n_skewers,width_Mpc=args.width_Mpc) + axis=args.axis,n_skewers=args.n_skewers,width_Mpc=args.width_Mpc) if args.verbose: print('print extra info about skewers')