From 70e88c592ff91d8008211bc62fa9a8c2c94d52e3 Mon Sep 17 00:00:00 2001 From: Josh Borrow <joshua.borrow@durham.ac.uk> Date: Fri, 9 Nov 2018 09:45:43 +0000 Subject: [PATCH] Added code snippet to show how to extract particles from a given halo --- .../source/VELOCIraptorInterface/output.rst | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/doc/RTD/source/VELOCIraptorInterface/output.rst b/doc/RTD/source/VELOCIraptorInterface/output.rst index b903a65012..cf0e386d4a 100644 --- a/doc/RTD/source/VELOCIraptorInterface/output.rst +++ b/doc/RTD/source/VELOCIraptorInterface/output.rst @@ -58,6 +58,64 @@ Besides the ``.catalog_particles`` file, there is also a but only for the unbound particles, a particle can only be present in one of these two lists. +Extracting the particles in a given halo +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The ``.catalog_particles`` file returns particle IDs that need to be matched +with those in your snapshot to find the particles in the file that you +wish to extract. The python snippet below should give you an idea of how to +go about doing this for the bound particles. + +First, we need to extract the offset from the ``.catalog_group`` file, and +work out how many _bound_ partices are in our halo. We can do this by +looking at the next offset. Then, we can ID match those with the snapshot +file, and get the mask for the _positions_ in the file that correspond +to our bound particles. (Note this requires ``numpy > 1.15.0``). + +.. code-block:: python + :linenos: + + import numpy as np + import h5py + + snapshot_file = h5py.File("swift_snapshot.hdf5", "r") + group_file = h5py.File("velociraptor_output.catalog_group", "r") + particles_file = h5py.File("velociraptor_output.catalog_particles", "r") + + halo = 100 + # Grab the start position in the particles file to read from + halo_start_position = group_file["Offset"][halo] + halo_end_position = group_file["Offset"][halo + 1] + # We're done with that file now, best to close earlier rather than later + group_file.close() + + # Get the relevant particle IDs for that halo; this includes particles + # of _all_ types. + particle_ids_in_halo = particles_file["Particle_IDs"][ + halo_start_position:halo_end_position + ] + # Again, we're done with that file. + particles_file.close() + + # Now, the tricky bit. We need to create the correspondance between the + # positions in the snapshot file, and the ids. + + # Let's look for the dark matter particles in that halo. + particle_ids_from_snapshot = snapshot_file["PartType1/ParticleIDs"][...] + + _, indices_v, indices_p = np.intersect1d( + particle_ids_in_halo, + particle_ids_from_snapshot, + assume_unique=True, + return_indices=True, + ) + + # indices_p gives the positions in the particle file where we will find + # the co-ordinates that we're looking for! To get the positions of all of + # those particles, + particle_positions_in_halo = snapshot_file["PartType1/Coordinates"][indices_p] + + Catalog_parttypes file ---------------------- -- GitLab