Commits

Britton Smith committed 6bfebf3

Updated get_dndz machinery to do variance fields.

Comments (0)

Files changed (2)

-from lightraytools import *
    files = glob("%s/*.h5" % data_dir)
    # Resolution here is the spectral resolution.  This will determine 
    # the number of ray elements blended together to form one absorber.
    # Range is the range in log column densities to be included. 
    # Format for range is, for example, range=[14, 20].
    # If range is set to None, as it is now, all absorbers are included.
    dn, dn_bins, total_rays, total_path, field_stats = dndl_rays(files, field, bins=bins, 
                                                                 path='dz', resolution=5000,
                                                                 stats_fields=stats_fields,
                                                                 filename=filename)
    output_file = "%s/%s.h5" % (data_dir, field)
    # This is the field in the ray data that has the ion fraction.
    fraction_field = "%sIon_Fraction" % (field[0:field.find('NumberDensity')])
    # Call the above wrapper for getting dn/dz.
    # stats_fields is a list of fields for which we want the mean properties for each column density bin.
    get_dndz(data_dir, field, bins=bins, stats_fields=['Metallicity', 'Temperature', 'Density', fraction_field],
             filename=output_file)
+from lightraytools import *
    files = glob("%s/*.h5" % data_dir)
    # Resolution here is the spectral resolution.  This will determine 
    # the number of ray elements blended together to form one absorber.
    # Range is the range in log column densities to be included. 
    # Format for range is, for example, range=[14, 20].
    # If range is set to None, as it is now, all absorbers are included.
    dn, dn_bins, total_rays, total_path, field_stats = dndl_rays(files, field, bins=bins, 
                                                                 path='dz', resolution=5000,
                                                                 mean_fields=mean_fields,
                                                                 variance_fields=variance_fields,
                                                                 filename=filename)
    output_file = "%s.h5" % field
    # This is the field in the ray data that has the ion fraction.
    fraction_field = "%sIon_Fraction" % (field[0:field.find('NumberDensity')])
    # Call the above wrapper for getting dn/dz.
    # stats_fields is a list of fields for which we want the mean properties for each column density bin.
    get_dndz(data_dir, field, bins=bins, 
             mean_fields=['Metallicity', 'Temperature', 'Density', fraction_field],
             variance_fields=['los_velocity'],
             filename=output_file)
 import numpy as na
 import h5py
 
-def dndl_rays(lightray_files, field, path='dz', bins=None, stats_fields=None,
+def dndl_rays(lightray_files, field, path='dz', bins=None, 
+              mean_fields=None, variance_fields=None,
               filename=None, **kwargs):
     "Calculate dn/dl for multiple light rays."
 
     dn_total = na.zeros((bins.size-1), dtype=float)
     dn_bins_total = None
     field_stats = None
-    if stats_fields is None: stats_fields = []
+    if mean_fields is None: mean_fields = []
+    if variance_fields is None: variance_fields = []
     total_path = 0.0
 
     for file in lightray_files:
-        dn, dn_bins, dpath, f_stats = dndl_ray(file, field, path=path, bins=bins, stats_fields=stats_fields, **kwargs)
+        dn, dn_bins, dpath, f_stats = dndl_ray(file, field, path=path, bins=bins, 
+                                               mean_fields=mean_fields, variance_fields=variance_fields,
+                                               **kwargs)
         bins = dn_bins
 
         if field_stats is None:
             field_stats = {}
-            for s_field in stats_fields + [field]:
+            for s_field in mean_fields + variance_fields + [field]:
                 field_stats[s_field] = [[] for bin in range(dn_bins.size - 1)]
 
-        for s_field in stats_fields + [field]:
+        for s_field in mean_fields + variance_fields + [field]:
             for bin in range(dn_bins.size - 1):
                 field_stats[s_field][bin].extend(f_stats[s_field][bin])
 
                                                       path, total_path)
 
     final_field_stats = {}
-    for s_field in stats_fields + [field]:
+    for s_field in mean_fields + variance_fields + [field]:
         final_field_stats[s_field] = na.zeros(((dn_bins.size - 1),2))
         for bin in range(dn_bins.size - 1):
             if s_field == field:
         output['dndl'].attrs['total_path'] = total_path
         output['dndl'].attrs['total_rays'] = len(lightray_files)
         output.create_dataset('dn_bins', data=dn_bins_total)
-        if stats_fields is not None:
-            for s_field in stats_fields + [field]:
-                output.create_dataset(s_field, data=final_field_stats[s_field])
+        for s_field in mean_fields + variance_fields + [field]:
+            output.create_dataset(s_field, data=final_field_stats[s_field])
         output.close()
 
     return (dn, dn_bins_total, len(lightray_files), total_path, final_field_stats)
 
 def dndl_ray(lightray_file, field, path='dz', resolution=None, cosmology=None, 
-             stats_fields=None, **kwargs):
+             mean_fields=None, variance_fields=None, **kwargs):
     "Calculte dn/dl for a single ray."
 
     field_stats = None
+    if mean_fields is None: mean_fields = []
+    if variance_fields is None: variance_fields = []
+
+    print lightray_file
 
     # Get data.
     input = h5py.File(lightray_file,'r')
     ray_data['z'] = input['redshift'].value
     ray_data['dz'] = input['dredshift'].value
     ray_data['dl'] = input['dl'].value
-    if stats_fields is not None:
-        for s_field in stats_fields:
-            ray_data[s_field] = input[s_field].value
+    for s_field in mean_fields + variance_fields:
+        ray_data[s_field] = input[s_field].value
     input.close()
 
     ray_data[field].size
     if resolution is None:
         ray_data[field] *= ray_data['dl']
     else:        
-        new_data = ray_resample(resolution, ray_data, fields=[field], weight_field=field, mean_fields=stats_fields)
+        new_data = ray_resample(resolution, ray_data, fields=[field], weight_field=field, 
+                                mean_fields=mean_fields, variance_fields=variance_fields)
         ray_data.update(new_data)
         del new_data
     del ray_data['dl']
 
     field_data_nz = ray_data[field] > 0
     ray_data[field] = ray_data[field][field_data_nz]
-    if stats_fields is not None:
-        for s_field in stats_fields:
-            ray_data[s_field] = ray_data[s_field][field_data_nz]
+    for s_field in mean_fields + variance_fields:
+        ray_data[s_field] = ray_data[s_field][field_data_nz]
     ray_data[field] = na.log10(ray_data[field])
     del field_data_nz
 
     dn, dn_bins = na.histogram(ray_data[field], **kwargs)
 
-    if stats_fields is not None:
-        field_stats = {field: []}
-        for s_field in stats_fields:
-            field_stats[s_field] = []
-        digits = na.digitize(ray_data[field], dn_bins)
-        digits[digits == dn_bins.size] -= 1
-        digits -= 1
-        for bin in range(dn_bins.size - 1):
-            indices = digits == bin
-            field_stats[field].append(ray_data[field][indices])
-            for s_field in stats_fields:
-                field_stats[s_field].append(ray_data[s_field][indices])
+    field_stats = {field: []}
+    for s_field in mean_fields + variance_fields:
+        field_stats[s_field] = []
+    digits = na.digitize(ray_data[field], dn_bins)
+    digits[digits == dn_bins.size] -= 1
+    digits -= 1
+    for bin in range(dn_bins.size - 1):
+        indices = digits == bin
+        field_stats[field].append(ray_data[field][indices])
+        for s_field in mean_fields + variance_fields:
+            field_stats[s_field].append(ray_data[s_field][indices])
 
     print "%s - %s: path (%s): %.2f, min = %.2f, max = %.2f." % \
         (lightray_file, field, path, dpath, ray_data[field].min(), ray_data[field].max())