Hyunmi Song
• 5
• 5 Apr '23

Hello,

To determine if a subhalo has undergone a major merger, I am attempting to measure the mass ratio using the MergerTree.
Then I found out about numMergers function in Sublink, and I have a question about how it works.
When I referred to this discussion forum( Major mergers identifications ), I could see that this function was as follows.
" There, a particular kind of mass ratio (defined as the maximum mass of some type, for both objects) is computed. "

When examining the definitions of fpMass and npMass, it appears that the stellar mass of each object is defined as the value at the time of its maximum. I don't understand why it calculates the mass ratio with the stellar mass at different times. In the Sublink , Rodriguez-Gomez+ (2015) that I referred to, it seems that fpMass uses the stellar mass of the FirstProgenitor at the same time(t_max of nextProgenitor) as when the mass of the nextProgenitor is at its maximum. Is there a reason why we don't use it this way?

The following is the result obtained by running the modified version of the existing code. In addition to the maximum mass returned by maxPastmass, it now also informs us of the time (Snapshot) when the maximum occurs.

< MaxPastMass >

``````def maxPastMass2(tree, index, partType='stars'):
""" Get maximum past mass (of the given partType) along the main branch of a subhalo
specified by index within this tree. """
ptNum =4
branchSize = tree['MainLeafProgenitorID'][index] - tree['SubhaloID'][index] + 1
masses = tree['SubhaloMassType'][index: index + branchSize, ptNum]

return np.max(masses), tree['SnapNum'][index: index + branchSize][np.argmax(masses)]
``````

< NumMergs >

``````numMergers=0
numNp=0
numMid=[]
index=0
minMassRatio=1/10
invMassRatio = 1.0 / minMassRatio
rootID=tree['SubhaloID'][index]
fpID=tree['FirstProgenitorID'][index]
massPartType='stars

while fpID != -1 :
fpIndex= index + (fpID-rootID)
fpMass,fpmax=maxPastMass2(tree,fpIndex,massPartType)
npID= tree['NextProgenitorID'][fpIndex]
numNp=0
while npID!= -1 :
npIndex = index + (npID - rootID)
npMass,npmax=maxPastMass2(tree,npIndex,massPartType)
if fpMass >0.0 and npMass >0.0 :
ratio=npMass/fpMass
if ratio >= minMassRatio and ratio <= invMassRatio:
numMergers += 1
print("#Sn=",tree['SnapNum'][fpIndex])
print("FpMass: ",fpMass,"at SN:",fpmax)
print("NpMass: ",npMass,"at SN:",npmax,"\n")

npID = tree['NextProgenitorID'][npIndex]
fpID = tree['FirstProgenitorID'][fpIndex]
``````
``````#Sn= 82
FpMass: 5.257069 at SN: 82
NpMass 0.5819944 at SN: 74

#Sn= 58
FpMass: 3.2475514 at SN: 58
NpMass: 0.60512143 at SN: 55

#Sn= 35
FpMass: 0.60444576 at SN: 35
NpMass: 0.4369604 at SN: 35

#Sn= 23
FpMass: 0.0460811 at SN: 23
NpMass: 0.0068318914 at SN: 22

#Sn= 19
FpMass:  0.008799356 at SN: 19
NpMass: 0.0015845126 at SN: 18

#Sn= 3
FpMass:  3.4521076e-05 at SN: 3
NpMass: 5.304361e-06 at SN: 2
``````

As seen in the above results, I can observe that the snapshots of fnMass and npMass used are different from each other.

I'll be grateful for your help .

Dylan Nelson
• 5 Apr '23

Hi Hyunmi,

You are correct, the `numMergers()` example code is just meant as an example, and it does not implement exactly the procedure described in Rodriguez-Gomez+15. I also agree that a better approach is to take the peak masses at the same time. I found an old email thread on this topic and paste it below.

Hi all,

I have been identifying interacting pairs within IllustrisTNG. Last week, I realized I could modify the `illustris_python.sublink.numMergers()` code to give me the information I needed. While I was modifying this piece of code, I noticed that there was a small error in the definition of the first progenitor's mass. I'm not sure how impactful this catch will be in the long run, but Lars thought I should share this with you all.

According to the Rodriguez-Gomez+2015 paper, the mass ratio is "based on the stellar masses of the two progenitors taken at t_max, i.e., at the time when the secondary progenitor reaches its maximum stellar mass." The way the code was originally written, the mass ratio is taken as the ratio of the maximum mass of the first progenitor to the maximum mass of the second progenitor. More often than not, these occur at different times. The correct way to do this would be to determine when the second progenitor has its maximum mass, and then determine the first progenitor's mass at that time to find the mass ratio.

I wrote a quick patch to the code (attached) that fixes the problem. To see exactly how different the two codes are, I ran the tree walk-through example from http://www.illustris-project.org/data/docs/scripts/ for the IllustrisTNG z=0 box.

For the 100-105th GroupFirstSub's, the original code gives:

``````for i in range(100, 105):
print GroupFirstSub[i], numMergers
....:
264412 6
265267 6
265926 4
266559 6
267323 4
``````

My code, run on the same set of galaxies, and using the same minimum mass ratio gives:

``````for i in range(100, 105):
print GroupFirstSub[i], numMergers
....:
264412 8
265267 8
265926 13
266559 10
267323 6
``````

To get an idea of how this works for a larger sample of galaxies, I also looked at the number of mergers that each version of the code found as a function of group mass for the first 10,000 GroupFirstSub galaxies. There doesn't seem to be a strong dependence on group mass, but this definitely shows that the original illustris_python.sublink.numMergers() code significantly underestimates the number of mergers that occur above a given mass ratio.

I'm not sure if this error is also in the IDL or Matlab versions of the code, but I suspect that it might be. I hope some of you find this useful in the future.

Cheers,
Kelly Blumenthal

Updated `sublink.py` file:

``````def maxPastMass(tree, index, partType='stars'):
""" Get maximum past mass (of the given partType) along the main branch of a subhalo
specified by index within this tree. """
ptNum = il.util.partTypeNum(partType)

branchSize = tree['MainLeafProgenitorID'][index] - tree['SubhaloID'][index] + 1
masses = tree['SubhaloMassType'][index: index + branchSize, ptNum]
snaps = tree['SnapNum'][index: index + branchSize]
return np.max(masses), snaps[np.where(masses == np.max(masses))[0]]

def massAtSnap(tree, index, snapnum, partType='stars'):
""" Get mass (of the given partType) along the main branch of a subhalo
at a specific snap shot. """
ptNum = il.util.partTypeNum(partType)

branchSize = tree['MainLeafProgenitorID'][index] - tree['SubhaloID'][index] + 1
masses = tree['SubhaloMassType'][index: index + branchSize, ptNum]
snaps = tree['SnapNum'][index: index + branchSize]

idx = np.where(snaps == snapnum)[0]
return masses[idx]

def numMergers(tree, minMassRatio=1e-10, massPartType='stars', index=0):
""" Calculate the number of mergers in this sub-tree (optionally above some mass ratio threshold). """
reqFields = ['SubhaloID', 'NextProgenitorID', 'MainLeafProgenitorID',
'FirstProgenitorID', 'SubhaloMassType', 'SnapNum']

if not set(reqFields).issubset(tree.keys()):
raise Exception('Error: Input tree needs to have loaded fields: '+', '.join(reqFields))

numMergers   = 0
invMassRatio = 1.0 / minMassRatio

# walk back main progenitor branch
rootID = tree['SubhaloID'][index]
fpID   = tree['FirstProgenitorID'][index]

while fpID != -1:
fpIndex = index + (fpID - rootID)
#removed: fpMass  = maxPastMass(tree, fpIndex, massPartType)

npID = tree['NextProgenitorID'][fpIndex]

while npID != -1:
npIndex = index + (npID - rootID)
#removed: npMass  = maxPastMass(tree, npIndex, massPartType)
npMass, npSnap  = maxPastMass(tree, npIndex, massPartType)

# count if both masses are non-zero, and ratio exceeds threshold
#removed: if fpMass > 0.0 and npMass > 0.0:
if npMass > 0.0:
fpMass = massAtSnap(tree, fpIndex, npSnap, massPartType)
if fpMass > 0.0:
ratio = npMass / fpMass

if ratio >= minMassRatio and ratio <= invMassRatio:
numMergers += 1

npID = tree['NextProgenitorID'][npIndex]

fpID = tree['FirstProgenitorID'][fpIndex]

return numMergers
``````
Hyunmi Song
• 5 Apr '23

Thank you so much for your help! I confirmed that the modified code you provided and the code I modified myself match perfectly.

Thank you again for your kind assistance. It was truly valuable and greatly appreciated

Hyunmi Song
• 1
• 3 May '23

Dear Dylan,
Additionally, I have a question. I would like to create a function similar to numMergers using LhaloTree.
However, it seems that LhaloTree does not have a function like numMergers, is that correct?
If the function doesn't exist, I'm trying to create it myself. But, I'm having difficulty finding the method in LhaloTree to find the next progenitor.

I understand the definition of NextProgenitor in LhaloTree to be as follows :
"The index of the next subhalo from the same snapshot which shares the same descendant, if any (-1 if this is the last). Indexes this TreeX group.".

Does this mean that the value of NextProgenitor refers to the direct index in this dataset when the value is obtained from the entire branch (not just the main branch) through loadTree?
However, when I try to apply this, I encounter the following error :

``````   [1] : T=il.lhalotree.loadTree('../../sims.TNG/TNG50-1/postprocessing/',99,96763,onlyMPB=False,fields=['FirstProgenitor','NextProgenitor','SubhaloNumber'])
[2] : T['NextProgenitor'][:99]
array([       -1, 153231757, 156898234, 156355761, 156356767,  68961306,
156929336, 156929331, 156929328, 156929327, 155974435, 156929308,
156929304, 153325597, 155975385, 152975405, 156929291, 156929281,
156929275, 153322855, 156929262, 156929202,  69389246,  72700207,
156929187,  92540908, 156929169, 156929165, 156929155, 156929131,
156929105, 156929091, 156929085, 156929082, 156929048, 156929023,
153480756, 153047794, 153324314, 156928951, 153319018, 156928931,
156928927, 153048509, 156928918, 156928904, 156927768, 156928812,
156928770, 156928747, 156928736, 156928734,  90392027, 156928724,
156928636, 156928615, 156928608, 156928603, 156928567, 156928515,
156928531, 156928518, 156928462, 156928427, 156928417, 156928411,
156928364, 156928337, 156928308, 156928298,  90669702, 156928253,
156928227,  90678054, 156928199, 156928134, 156928119, 156928108,
156928088, 156928067, 156928053, 156928035, 156928015,  90675068,
156927969, 156927959, 156927950, 156927945, 156927916,  90674139,
156927913, 156927903, 156927894, 156927893,  90674085, 156927885,
156927883, 156927882,  51755329], dtype=int32)
[3] :  T['SubhaloNumber'][T['NextProgenitor'][1]]
--------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-52-ee6788dfa004> in <module>
----> 1 T['SubhaloNumber'][T['NextProgenitor'][1]]

IndexError: index 153231757 is out of bounds for axis 0 with size 94604
``````

The way of finding the nextprogenitor in SubLink was as follows :

``````  < In numMergers function >
rootID=tree['SubhaloID'][index]
fpID   = tree['FirstProgenitorID'][index]
while fpID != -1:
fpIndex = index + (fpID - rootID)
npID = tree['NextProgenitorID'][fpIndex]
while npID != -1:
npID = tree['NextProgenitorID'][npIndex]
fpID = tree['FirstProgenitorID'][fpIndex]
``````

and I'm sure there must be a similar way of finding it in LhaloTree, but I'm not sure what it is.

Please let me know if there is any issue with the approach I have taken! Thank you as always.

Dylan Nelson
• 3 May '23

I am not using LHaloTree much, but I think the issue is the difference between these two statements:

(1) "Indexes this TreeX group."

(2) "Does this mean that the value of NextProgenitor refers to the direct index in this dataset when the value is obtained from the entire branch (not just the main branch) through loadTree?"

I would have agreed with you (it would make sense like #2), but it isn't exactly the case. What is written in #1 is correct.

In particular, if you actually look at the LHaloTree files you see many groups named `Tree{X}`. As the documentation explains, one of these groups can contain subhalos of different halos, and what subhalo is where is a bit confusing. The `loadTree()` function handles this, and returns to you the MPB or the full tree for that subhalo. However, it isn't immediately obvious what `Tree{X}` group it comes from.

In your particular case above, this subhalo is inside `Tree0` of `trees_sf1_099.0.hdf5`, and all the datasets of this tree have a size of `218547998`. So indeed you will find all your NextProgenitor indices inside the datasets of this tree (`153231757` will be in bounds).

In practice, then, you could use the `treeOffsets()` function to get the file and tree number (the first and third returns, respectively). Then, you would probably want to manually load (with h5py) the datasets of this `Tree{X}` group, from the given file, and then follow the NextProgenitor links as you have above.

• Page 1 of 1