diff --git a/t/hej_test.hh b/t/hej_test.hh
index bd53558..d8b1009 100644
--- a/t/hej_test.hh
+++ b/t/hej_test.hh
@@ -1,79 +1,95 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <array>
 #include <string>
 #include <vector>
 
 #include "HEJ/Event.hh"
 #include "HEJ/exceptions.hh"
 
 namespace HEJ {
   struct Particle;
 }
 
 //! throw error if condition not fulfilled
 #define ASSERT(x) if(!(x)) { \
     throw std::logic_error("Assertion '" #x "' failed."); \
   }
 
 //! throw error if prop is different between ev1 and ev2
 #define ASSERT_PROPERTY(ev1,ev2,prop) ASSERT(ev1.prop == ev2.prop)
 
 //! throw error if condition not fulfilled
 #define ASSERT_THROW(x, exception) try { \
       x; \
       std::cerr << "'" #x "' did not throw an exception.\n"; \
       throw; \
     } catch(exception const &){} \
     catch (...) { \
       std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \
       throw; \
     }
 
 /** @brief get specific Phase Space Points for njets with boson at pos_boson
  *
  * if pos_boson = -1 (or not implemented) -> no boson
  *
  * njet==7 is special: has less jets, i.e. multiple parton in one jet,
  *                     all partons are massive (4 GeV) -> can be boson/decay
  *                     pos_boson < 0 to select process (see list for details)
  */
 HEJ::Event::EventData get_process(int njet, int pos_boson);
 
 //! select process from string input (see also get_process)
 //!
 //! overwrite_boson to force a specific boson position, indepentent from input
 //! (useful for njet == 7)
 HEJ::Event::EventData parse_configuration(
     std::array<std::string,2> const & in, std::vector<std::string> const & out,
     int overwrite_boson = 0
 );
 
-/** @brief generate a rapidity ordered event
-*/
+/**
+ * @brief Generate a rapidity ordered event (with decays)
+ * @param in             Incoming particle flavours
+ * @param out            Outgoing particle flavours
+ * @param reconstruct    Should call reconstruct_intermediate?
+ * @param decays         Any decays associated with particles from out
+ * @returns              EventData for in -> out with rapidity ordered final state
+ *
+ * Generates EventData for the process in -> out. The outgoing particles
+ * are generated in rapidity order. The incoming flavours are ordered:
+ * {backward, forward}.
+ *
+ * Decays are provided by decays, if specificed, and via reconstruct_intermediate
+ * if requested.
+ *
+ * Note: Specified decays are handled before reconstruct_intermediate is called.
+ */
+
 HEJ::Event::EventData rapidity_order_ps(
     std::array<std::string,2> const & in,
     std::vector<std::string> const & out,
     bool reconstruct = false,
     std::unordered_map<size_t, std::vector<std::string> > decays = {}
   );
 
 //! shuffle particles around
 void shuffle_particles(HEJ::Event::EventData & ev);
 
 //! Helper function to couple quarksfor flavour-changing bosons
 bool couple_quark(std::string const & boson, std::string & quark);
 
 //! Decay kinematics for 1->2
 std::vector<HEJ::Particle> decay_kinematics( HEJ::Particle const & parent );
 
 //! Decay W boson to lepton & neutrino
 std::vector<HEJ::Particle> decay_W( HEJ::Particle const & parent );
 
 //! Decay Z to electron-positron
 std::vector<HEJ::Particle> decay_Z( HEJ::Particle const & parent );