Add synfast galaxy map for large scale structure background test

dguevel · Aug 20, 2024 · 931c51d · 931c51d
1 parent 7ba3dfa
commit 931c51d
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 2 deletions.
diff --git a/KIPAC/nuXgal/GalaxySample.py b/KIPAC/nuXgal/GalaxySample.py
@@ -453,6 +453,33 @@ def __init__(self):
         GalaxySample.__init__(self, "WISE", self.mask())
 
 
+class Synfast(GalaxySample):
+    """WISE Galaxy sample
+
+    Synfast instance of WISE-2MASS galaxy sample map based on ~5M galaixes
+    """
+    @staticmethod
+    def mask():
+        """Contstruct and return the mask for this sample"""
+        c_icrs = SkyCoord(ra=Defaults.exposuremap_phi * u.radian,
+                          dec=(np.pi/2 - Defaults.exposuremap_theta)*u.radian, frame='icrs')
+
+        # planck dust map
+        planck_dustmap_file = os.path.join(Defaults.NUXGAL_ANCIL_DIR, 'HFI_Mask_GalPlane-apo0_2048_R2.00_nside128.fits')
+        planck_dustmap = hp.read_map(planck_dustmap_file)
+
+        # large and small magellanic clouds
+        lmc = hp.query_disc(Defaults.NSIDE, hp.ang2vec(80.894200, -69.756100, lonlat=True), 5*np.pi/180)
+        smc = hp.query_disc(Defaults.NSIDE, hp.ang2vec(13.158300, -72.800300, lonlat=True), 5*np.pi/180)
+        mc = np.zeros(Defaults.NPIXEL)
+        mc[lmc] = 1
+        mc[smc] = 1
+        return np.where((np.abs(c_icrs.galactic.b.degree) < 10) | (planck_dustmap == 0) | (mc == 1))
+
+    def __init__(self):
+        """C'tor"""
+        GalaxySample.__init__(self, "synfast", self.mask())
+
 class GalaxySample_unWise_z04(GalaxySample):
     """WISE Galaxy sample
 
@@ -596,6 +623,7 @@ class GalaxySampleLibrary:
         'unWISE_z=0.6': GalaxySample_unWise_z06,
         'unWISE_z=1.0': GalaxySample_unWise_z10,
         'unWISE_z=1.5': GalaxySample_unWise_z15,
+        'synfast': Synfast,
         'Atmospheric': GalaxySample_Atmospheric,
         'SDSS_z0.0_0.1': GalaxySample_SDSS_00_01,
         'SDSS_z0.1_0.2': GalaxySample_SDSS_01_02,

diff --git a/KIPAC/nuXgal/data/ancil/synfast_galaxymap.fits b/KIPAC/nuXgal/data/ancil/synfast_galaxymap.fits
diff --git a/scripts/nuXgal-generate-data.py b/scripts/nuXgal-generate-data.py
@@ -107,6 +107,7 @@ def main():
                         help='Path to MC signal file if not using default')
     parser.add_argument('--isotropic', type=float, default=0,
                         help='Add an isotropic component by fraction of measured')
+    parser.add_argument('--inj-template', default='')
 
     args = parser.parse_args()
 
@@ -158,13 +159,30 @@ def main():
         )
         isotropic_counts = isotropic_event_generator.trial_runner.to_ns(
             (1e5)**2 * 1.44e-18 * llh.f_sky * 4 * np.pi, E0=1e5, unit=1)
+
+    if args.inj_template:
+        inj_galaxy_sample = GALAXY_LIBRARY.get_sample(args.inj_template)
+        inj_event_generator = CskyEventGenerator(
+            args.N_yr,
+            inj_galaxy_sample,
+            gamma=2.28,
+            Ebinmin=args.ebinmin,
+            Ebinmax=args.ebinmax,
+            idx_mask=llh.idx_mask,
+            mc_background=args.mcbg,
+            path_sig=args.path_sig
+        )
 
     result_list = []
 
     for ninj in args.inject:
         for i in tqdm(np.arange(args.n_trials)):
-            trial = llh.event_generator.SyntheticTrial(
-                ninj, keep_total_constant=False)[0]
+            if args.inj_template:
+                trial = inj_event_generator.SyntheticTrial(
+                    ninj, keep_total_constant=False)[0]
+            else:
+                trial = llh.event_generator.SyntheticTrial(
+                    ninj, keep_total_constant=False)[0]
             if args.pnt_src:
                 pt_trial = pnt_trial_runner.get_one_trial(ninj_pt)[0]
                 trial = append_signal_trials(trial, pt_trial)