Reap output sequence depends upon inputs-- How do I specify the order of reap output?

Question

I am simulating two groups of interacting oscillators, swings and metronomes. When I perform the NDSolve, I also extract the peaks of the swings and the metronomes. This was going fine until I changed my definition of phi'[0]. With one definition, the swing peaks are assigned as sol[[2,1]] and with another definition, the metronome peaks are assigned as sol[[2,1]] (with the other set of peaks assigned to sol[[2,2]] in each case).

I tried to assign the swing peaks as the first assignment (that is, sol[[2,1]]), but obviously this isn't working. Code is below, sorry that it is longer than I would prefer but I was having trouble reproducing the issue with simpler code.

Why does the Reap seem to assign the WhenEvent outputs in an arbitrary sequence, rather than the order listed in the NDSolve function? Can I hardcode the sequence to avoid this problem (as I tried to, but failed)? I assume I'm missing some fundamental behavior of Reap, if so, what is it? I saw a post referring to the second input of Sow as a "tag" (Defining Tags in Reap) but I'm having trouble seeing the way through still.

Options[NDSolve`IDA] (* Define the IDA method of NDSolve *)
tmax = 500; (* maximum running time for the simulation *)
SW = 2; (* number of swings *)
MET = 5; (* number of metronomes *)

randNxM = RandomReal[{-.33, .33}, MET];
syncNxM =  ConstantArray[.66, MET];
dsyncNxM =  ConstantArray[0, MET];

(* assemble the initial conditions from vectors by swing *)
inits = Transpose[Join[{randNxM, syncNxM}, 2]];
dinits = Transpose[Join[{randNxM, dsyncNxM}, 2]];

(*------ Dynamic equations --------------------------------*)

swingEqns = Table[
Theta[r]''[tt] == 
0.426 Sum[(Theta[j][tt] - Theta[r][tt]), {j, SW}] - 
 0.606 Theta[r][tt] - 0.00016 Theta[r]'[tt] -
 Sum[D[Sin[phi[l][r][tt]], {tt, 2}], {l, MET}], {r, SW}];
swingInits = Table[{Theta[r][0] == .66, Theta[r]'[0] == 0}, {r, SW}];

metrEqns = Table[
phi[q][r]''[tt] == -Sin[phi[q][r][tt]] - 
 0.0105 phi[q][r]'[tt] ((phi[q][r][tt]/.33)^2 - 1) - 
 0.000489 Cos[phi[q][r][tt]] Theta[r]''[tt], {q, MET}, {r, SW}];

(* two initial conditions here *)
metrInitCondZero = 
Table[{phi[q][r][0] == inits[[q, r]], phi[q][r]'[0] == 0}, {q, 
MET}, {r, SW}];
metrInitCondDinits = 
Table[{phi[q][r][0] == inits[[q, r]], 
phi[q][r]'[0] == dinits[[q, r]]}, {q, MET}, {r, SW}];

varsMetr = Flatten[Table[phi[q][r], {q, MET}, {r, SW}]];
varsSwing = Table[Theta[r], {r, SW}];

(* defineEvent defines the peak finding conditional. 
The 1 and 2 at the end of the Sow functions are meant to 
assign SW as the first Reap output and MET as the second Reap output *)
defineEventsSW = Table[With[{j = j},
WhenEvent[Theta[j]'[tt] == 0 && Theta[j][tt] > 0,
 Sow[{j, tt, Theta[j][tt]}, 1]]],
{j, SW}]; 
defineEventsMET = Table[With[{i = i, j = j},
WhenEvent[phi[i][j]'[tt] == 0 && phi[i][j][tt] > 0,
 Sow[{i, j, tt, phi[i][j][tt]}, 2]]],
{i, MET}, {j, SW}]; 

(* solution with phi'[t=0]==0 *)
solZero = Reap@NDSolve[{
 {swingEqns, metrEqns},
 {swingInits, metrInitCondZero},
 {defineEventsSW,
  defineEventsMET}
 },
Union[varsSwing, varsMetr],
{tt, 0, tmax}];

(* solution with phi'[t=0]==dinits *)
solDinits = Reap@NDSolve[{
 {swingEqns, metrEqns},
 {swingInits, metrInitCondDinits},
 {defineEventsSW,
  defineEventsMET}
 },
Union[varsSwing, varsMetr],
{tt, 0, tmax}];

(* Check dimensions of inputs and outputs*)
(* INITIAL CONDITIONS OF THE SAME DIMENSION PRODUCE DIFFERENT \ SOLUTION OUTPUTS!!! *)
(* 2 DIFFERENT INITIAL CONDITION MATRICES *)
Dimensions[metrInitCondZero]
Dimensions[metrInitCondDinits]

Which gives {5, 2, 2} and {5, 2, 2}

(* THE FIRST PART OF THE REAPED RESULT *)
Dimensions@solZero[[2, 1]]
Dimensions@solDinits[[2, 1]]

Which gives {165, 3} and {777, 4} (exact numbers may vary, but note that the swing peaks should have a 3 in the second place and the metronomes should have a 4.)

(* THE SECOND PART OF THE REAPED RESULT *)
Dimensions@solZero[[2, 2]]
Dimensions@solDinits[[2, 2]]

Which gives {777, 4} and {155, 3}. Again note the second value.

Answer 1

As a much more simple illustration of the problem consider:

 Reap[Do[ x = RandomReal[]; 
    If[x < 1/2, Sow[x, 1], Sow[x, 2]], {10}]][[2]]

{{0.482192, 0.0468871, 0.0211528, 0.0476038, 0.424523, 0.152544}, {0.917763, 0.510109, 0.604724, 0.557907}}

If you evaluate this repeatedly you will see the ordering of the Reap results changes at random depending on which Sow is encountered first.

To ensure a consistent ordering, when you use the two argument form of Sow you should also give the second argument to Reap as:

Reap[Do[ x = RandomReal[]; 
   If[x < 1/2, Sow[x, 1], Sow[x, 2]], {10}], {1, 2}][[2]]

or sometimes its more clear if you use some expression for the Reap/Sow pattern:

Reap[Do[ x = RandomReal[]; 
   If[x < 1/2, Sow[x,"low"], Sow[x,"high"]], {10}], {"low", "high"}][[2]]