Reading lightcone halo catalogues
=================================

Approximate lightcone halo catalogues can be generated by running a
halo finder on many simulation snapshots and using particles (such as
black hole particles) from the particle lightcone outputs as tracers to
determine when each halo crosses the lightcone. This is implemented in
the ``lightcone_io.match_black_holes`` module. See
:ref:`halo_lightcones` for details.

Each halo lightcone file contains halos taken from a single simulation
snapshot, which corresponds to a range of redshifts centred on the
snapshot redshift in the lightcone. Within the file the halos are
sorted in order of their pixel index in a low resolution HEALPix
map. This makes it possible to extract halos in particular parts of
the sky.

The halos can be stored in HDF5 files on the local file system or in
remote files accessed via a web service.

Opening local HDF5 files
------------------------

A halo lightcone file stored on the local file system can be opened
using the :py:class:`lightcone_io.HaloLightconeFile` class::

  import lightcone_io as lc
  halos = lc.HaloLightconeFile(filename)

Opening via the hdfstream service
---------------------------------

This module can also access files stored on a remote server using the
`hdfstream <https://hdfstream-python.readthedocs.io/en/latest>`_
python module. To do this, first open the directory containing the
files::

  import hdfstream
  root = hdfstream.open("cosma", "/")

This returns a :obj:`hdfstream.RemoteDirectory` object. Halo lightcone
outputs on the server can be opened by passing the remote directory
object to the :py:class:`lightcone_io.HaloLightconeFile` class. The
filename is interpreted as a path relative to the remote directory on
the server::

  import lightcone_io as lc
  halos = lc.HaloLightconeFile(filename, remote_dir=root)

Accessing halo lightcone data
-----------------------------

We can get a list of the available halo properties with::

  print(halos.properties)

To read positions and bound masses for all halos in the catalogue, for
example, we can do this::

  # List of halo properties to read
  properties = ("Lightcone/HaloCentre", "BoundSubhalo/TotalMass")

  # Read the data
  halo_props = halos.read_halos(properties)

This returns a dict of unyt arrays, which are the halo properties with
unit information attached. Since the halos are sorted by HEALPix
pixel, it's possible to extract halos close to a specified position on
the sky::

  # Line of sight vector specifying a point on the sky
  vector = (1.0, 0.0, 0.0)

  # Angular radius around this point (in radians)
  import numpy as np
  radius = np.radians(10.0)

  # List of halo properties to read
  properties = ("Lightcone/HaloCentre", "BoundSubhalo/TotalMass")

  # Read the data
  halo_props = halos.read_halos_in_radius(vector, radius, properties)

Again, this returns a dict of unyt arrays with the halo properties,
but only halos in the specified patch of sky are included.

Reading SOAP properties for lightcone halos
-------------------------------------------

Not all SOAP halo properties are copied into the lightcone halo
catalogue because it would make the files unreasonably large. However,
for each halo in the halo lightcone we store the index in the
corresponding SOAP output as ``InputHalos/SOAPIndex``. It is therefore
possible to look up a wide range of SOAP properties for halos in the
lightcone catalogue.

The :py:class:`lightcone_io.HaloLightconeFile` class has an optional
parameter ``soap_filename`` which can be used to do this. For
example::

      import lightcone_io as lc
      halos = lc.HaloLightconeFile(filename="hbt_lightcone_halos/lightcone0/lightcone_halos_0070.hdf5",
                                   soap_filename="SOAP-HBT/halo_properties_0070.hdf5")

This opens one of the halo lightcone files AND the corresponding SOAP
output which the halos were taken from. We can then select a patch of
sky to read in, as above::

      # Line of sight vector specifying a point on the sky
      vector = (1.0, 0.0, 0.0)

      # Angular radius around this point (in radians)
      import numpy as np
      radius = np.radians(10.0)

and read in some halo properties::

    # List of halo properties to read
    properties = ("Lightcone/HaloCentre", "SO/200_crit/TotalMass")

    # Read the data
    halo_props = halos.read_halos(properties)

As before, this returns a dictionary of unyt arrays with the halo
properties but any properties which do not exist in the halo lightcone
file will be read in from SOAP instead. In this case the quantity
``SO/200_crit/TotalMass`` is read from the SOAP file and reordered to
match the lightcone halo catalogue.