diff --git a/dmpipe/dm_prepare.py b/dmpipe/dm_prepare.py index 785e740..b9448c1 100644 --- a/dmpipe/dm_prepare.py +++ b/dmpipe/dm_prepare.py @@ -9,6 +9,7 @@ from __future__ import absolute_import, division, print_function import os +import sys import copy from shutil import copyfile @@ -167,6 +168,13 @@ def _write_profile_yaml(cls, target, profile_path, targ_ver, spatial): source_model.update(dict(SpatialModel='DiffuseSource', SpatialType='SpatialMap', Spatial_Filename=target.j_map_file)) + elif spatial in ['hpxmap']: + if target.j_map_file is None: + j_map_file = profile_path.replace('.yaml', '.fits') + target.write_jmap_hpx(j_map_file) + source_model.update(dict(SpatialModel='DiffuseSource', + SpatialType='SpatialMap', + Spatial_Filename=target.j_map_file)) elif spatial in ['radial']: target.j_rad_file = profile_path.replace('.yaml', '.dat') target.write_j_rad_file() @@ -180,6 +188,13 @@ def _write_profile_yaml(cls, target, profile_path, targ_ver, spatial): source_model.update(dict(SpatialModel='DiffuseSource', SpatialType='SpatialMap', Spatial_Filename=target.d_map_file)) + elif spatial in ['hpxdmap']: + if target.d_map_file is None: + d_map_file = profile_path.replace('.yaml', '.fits') + target.write_dmap_hpx(d_map_file) + source_model.update(dict(SpatialModel='DiffuseSource', + SpatialType='SpatialMap', + Spatial_Filename=target.d_map_file)) elif spatial in ['dradial']: target.d_rad_file = profile_path.replace('.yaml', '.dat') target.write_d_rad_file() diff --git a/dmpipe/dm_spectral.py b/dmpipe/dm_spectral.py index b4b526c..2a94f27 100644 --- a/dmpipe/dm_spectral.py +++ b/dmpipe/dm_spectral.py @@ -60,12 +60,12 @@ class ConvertCastro(Link): @staticmethod def is_decay_sed(sedfile): tokens = os.path.splitext(os.path.basename(sedfile))[0].split('_') - return tokens[2] in ['point', 'dmap', 'dradial'] + return tokens[2] in ['point', 'dmap', 'dradial', 'hpxdmap'] @staticmethod def is_ann_sed(sedfile): tokens = os.path.splitext(os.path.basename(sedfile))[0].split('_') - return tokens[2] in ['point', 'map', 'radial'] + return tokens[2] in ['point', 'map', 'radial', 'hpxmap'] @staticmethod def select_channels(channels, sedfile): diff --git a/dmpipe/dm_spectral_utils.py b/dmpipe/dm_spectral_utils.py index d750e7c..fb8765c 100644 --- a/dmpipe/dm_spectral_utils.py +++ b/dmpipe/dm_spectral_utils.py @@ -696,6 +696,114 @@ def x_edges(self): return np.power(10, log_masses) + def write_glory_duck_v1(self, basepath, **kwargs): + """Write this object into the format that the 'Glory Duck' collaboration wants + + Paramters + --------- + + basepath : str + Path to output files + """ + kwcopy = kwargs.copy() + sigmav_min = kwcopy.get('logsigmav_min', -30) + sigmav_max = kwcopy.get('logsigmav_max', -24) + #sigmav_nstep = kwcopy.get('sigmav_nstep', 20000) + sigmav_nstep = kwcopy.get('sigmav_nstep', 200) + + sigmav_steps = np.logspace(sigmav_min, sigmav_max, sigmav_nstep)[::-1] + + for imass, mass in enumerate(self.masses): + filepath = "%s_%0.1f_GeV.txt" % (basepath, mass) + + ll_interp = self[imass].interp + ll_vals = ll_interp(sigmav_steps) + + outfile = open(filepath, 'w!') + + sys.stdout.write("Writing %i values to %s: " % (sigmav_nstep, filepath)) + sys.stdout.flush() + for i, (sigmav, ll) in enumerate(zip(sigmav_steps, ll_vals)): + if i % 100 == 0: + sys.stdout.write('.') + sys.stdout.flush() + outfile.write("%6e %6e\n" % (sigmav, ll)) + outfile.close() + sys.stdout.write('!\n') + + + def write_glory_duck_v2(self, filepath, **kwargs): + """Write this object into the format that the 'Glory Duck' collaboration wants + + Paramters + --------- + + basepath : str + Path to output files + """ + kwcopy = kwargs.copy() + sigmav_min = kwcopy.get('logsigmav_min', -30) + sigmav_max = kwcopy.get('logsigmav_max', -24) + #sigmav_nstep = kwcopy.get('sigmav_nstep', 20000) + sigmav_nstep = kwcopy.get('sigmav_nstep', 200) + mass_vals = kwcopy.get('mass_vals', None) + + sigmav_steps = np.logspace(sigmav_min, sigmav_max, sigmav_nstep)[::-1] + + ll_vals_list = [] + if mass_vals is not None: + mass_list = mass_vals + else: + mass_list = self.masses + + # Do linear interpolation in log-space + log_masses = np.log10(self.masses) + log_mass_diffs = log_masses[1:] - log_masses[0:-1] + log_mass_vals = np.log10(mass_list) + + for log_mass in log_mass_vals: + mass_index = np.searchsorted(log_masses, log_mass) + if mass_index < 0: + ll_vals_list.append(np.zeros((sigmav_nstep))) + continue + elif mass_index >= len(log_masses): + ll_vals_list.append(np.zeros((sigmav_nstep))) + continue + elif mass_index == len(log_masses) - 1: + ll_interp = self[mass_index].interp + ll_vals = ll_interp(sigmav_steps) + ll_vals_list.append(ll_vals) + continue + frac_lo = (log_masses[mass_index+1] - log_mass) / log_mass_diffs[mass_index] + frac_hi = 1 - frac_lo + ll_interp_lo = self[mass_index].interp + ll_vals_lo = ll_interp_lo(sigmav_steps) + ll_interp_hi = self[mass_index+1].interp + ll_vals_hi = ll_interp_hi(sigmav_steps) + ll_vals = frac_lo*ll_vals_lo + frac_hi*ll_vals_hi + ll_vals_list.append(ll_vals) + + ll_vals_stack = np.vstack(ll_vals_list).T + + outfile = open(filepath, 'w!') + sys.stdout.write("Writing %i values to %s: " % (sigmav_nstep, filepath)) + + outfile.write("%.2f" % np.log10(self._astro_value)) + for mass in mass_list: + outfile.write(" %.2f" % (mass)) + outfile.write("\n") + sys.stdout.flush() + for i, (sigmav, ll_vals) in enumerate(zip(sigmav_steps, ll_vals_stack)): + if i % 1000 == 0: + sys.stdout.write('.') + sys.stdout.flush() + outfile.write("%6e " % (sigmav)) + for ll_val in ll_vals: + outfile.write("%6e " % (ll_val)) + outfile.write("\n") + outfile.close() + sys.stdout.write('!\n') + class DMSpecTable(object): """ Version of the DM spectral tables in tabular form @@ -857,6 +965,8 @@ def ebin_refs(self): """ return self._e_table["E_REF"].data + + def write_fits(self, filepath, clobber=False): """ Write this object to a FITS file