clusters/main: Implement _identify_halos() to determine genuine halos

New configuration options 'extragalactic/halos/genuine_{emfacc,index}_th'
are added to set the thresholds.

2 files changed, 30 insertions, 0 deletions

diff --git a/fg21sim/configs/config.spec b/fg21sim/configs/config.spec
index fd45949..5a0731d 100644
--- a/fg21sim/configs/config.spec
+++ b/fg21sim/configs/config.spec
@@ -456,3 +456,8 @@ stream = option("stderr", "stdout", "", default="stderr")
   f_rc = float(default=0.1)
   # The slope parameter (i.e., beta).
   beta = float(default=0.6667)
+
+  # The emissivity acceleration factor and spectral index thresholds for
+  # determining whether the halo is genuine/formed.
+  genuine_emfacc_th = float(default=100)
+  genuine_index_th = float(default=3.5)
diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py
index 4e1f134..cfeaf80 100644
--- a/fg21sim/extragalactic/clusters/main.py
+++ b/fg21sim/extragalactic/clusters/main.py
@@ -93,6 +93,8 @@ class GalaxyClusters:
 
         sec = "extragalactic/halos"
         self.eta_b = configs.getn(sec+"/x_cr")
+        self.genuine_emfacc_th = configs.getn(sec+"/genuine_emfacc_th")
+        self.genuine_index_th = configs.getn(sec+"/genuine_index_th")
 
         if self.use_dump_halos_data and (not self.use_dump_catalog_data):
             self.use_dump_catalog_data = True
@@ -387,6 +389,28 @@ class GalaxyClusters:
 
         logger.info("Calculated halo emissions.")
 
+    def _identify_halos(self):
+        """
+        Determine the formation/genuineness of each radio halo.
+
+        A halo is genuine if both its emissivity acceleration factor >=
+        genuine_emfacc_th and its spectral index <= genuine_index_th
+        at any frequency.
+        """
+        logger.info("Identify the genuineness of radio halos ...")
+        emfacc_th = self.genuine_emfacc_th
+        index_th = self.genuine_index_th
+        n = 0
+        for hdict in self.halos:
+            emfacc = hdict["emissivity_facc"]
+            index = hdict["spec_index"]
+            genuine = np.any((emfacc >= emfacc_th) & (index <= index_th))
+            hdict["genuine"] = genuine
+            n += genuine
+
+        logger.info("Identified %d (%.1f%%) genuine halos." %
+                    (n, n*100/len(self.halos)))
+
     def _dropout_halos(self):
         """
         Considering that the (very) massive galaxy clusters are very rare,
@@ -588,6 +612,7 @@ class GalaxyClusters:
             self._simulate_halos()
 
         self._calc_halos_emission()
+        self._identify_halos()
         self._dropout_halos()
         self._draw_halos()