From dc6d330eef60b7ca10e029d9a9af434454575daa Mon Sep 17 00:00:00 2001
From: Claudia Lagos <claudia.lagos@icrar.org>
Date: Tue, 21 Dec 2021 16:10:06 +0800
Subject: [PATCH] Adding most bound BH ID and position to output

Added the position and ID of the BH that is closest to the most bound dark matter particle.
This was tested in COLIBRE and works well.
---
 src/allvars.cxx                |  3 +++
 src/allvars.h                  |  5 +++++
 src/io.cxx                     |  9 +++++++++
 src/substructureproperties.cxx | 16 ++++++++++++++++
 4 files changed, 33 insertions(+)

diff --git a/src/allvars.cxx b/src/allvars.cxx
index 5de1a9bd..c8f78a8e 100644
--- a/src/allvars.cxx
+++ b/src/allvars.cxx
@@ -1811,6 +1811,9 @@ void PropDataHeader::declare_all_datasets(const Options &opt)
 #ifdef BHON
     declare_dataset<unsigned long>("n_bh");
     declare_dataset("Mass_bh", MASS);
+    //save ID most bound BH and its position
+    declare_dataset<long>("ID_mbp_bh");
+    declare_XYZ_datasets("cmbp_bh", LENGTH);
 #endif // BHON
 
 #ifdef HIGHRES
diff --git a/src/allvars.h b/src/allvars.h
index 1ed6d87f..415bdf48 100644
--- a/src/allvars.h
+++ b/src/allvars.h
@@ -1621,6 +1621,10 @@ struct PropData
     ///mean accretion rate, metallicty
     Double_t acc_bh, acc_bh_mostmassive;
 
+    ///ID of most bound BH and its position
+    Int_t ibound_bh = 0;
+    Coordinate gposmbp_bh;
+
     ///\name blackhole aperture/radial profiles
     //@{
     vector<int> aperture_npart_bh;
@@ -2478,6 +2482,7 @@ struct PropData
 #endif
 #ifdef BHON
         M_bh*=opt.h;
+        gposmbp_bh*=opt.h/opt.a;
 #endif
 #ifdef HIGHRES
         M_interloper*=opt.h;
diff --git a/src/io.cxx b/src/io.cxx
index 3cb3d92a..02c47917 100644
--- a/src/io.cxx
+++ b/src/io.cxx
@@ -2177,6 +2177,15 @@ void WriteProperties(Options &opt, const Int_t ngroups, PropData *pdata){
 
         for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].M_bh;
         Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (Int_t i=0;i<ngroups;i++) ((long long*)data)[i]=pdata[i+1].ibound_bh;
+        Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (int k=0;k<3;k++){
+            for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].gposmbp_bh[k];
+            Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+        }
+
 #endif
 #ifdef HIGHRES
         for (Int_t i=0;i<ngroups;i++) ((unsigned long*)data)[i]=pdata[i+1].n_interloper;
diff --git a/src/substructureproperties.cxx b/src/substructureproperties.cxx
index 90489ec2..f5da189a 100644
--- a/src/substructureproperties.cxx
+++ b/src/substructureproperties.cxx
@@ -1136,6 +1136,14 @@ private(EncMassSF,EncMassNSF,Krot_sf,Krot_nsf,Ekin_sf,Ekin_nsf)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+		//save particle information if the most bound id has not been assigned.
+		//note that this works because particles has been previously ordered by distance to defined centre.
+		if(pdata[i].ibound_bh == 0){
+			pdata[i].ibound_bh = Pval->GetPID();
+			for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+		}
+
             }
         }
 #endif
@@ -2073,6 +2081,14 @@ private(j,Pval,x,y,z,vx,vy,vz,jval,jzval,zdist,Rdist)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+                //save particle information if the most bound id has not been assigned.
+                //note that this works because particles has been previously ordered by distance to defined centre.
+                if(pdata[i].ibound_bh == 0){
+                        pdata[i].ibound_bh = Pval->GetPID();
+                        for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+                }
+
             }
         }
 #endif