From 8281f0bed612245f1243a669f1a992644c865f08 Mon Sep 17 00:00:00 2001
From: Oliver Hahn <oliverjhahn@gmail.com>
Date: Fri, 25 Mar 2022 18:17:10 +0100
Subject: [PATCH] added customization to particle masking

---
 example.conf                  |  3 +++
 include/particle_generator.hh | 30 ++++++++++++++++++++++--------
 2 files changed, 25 insertions(+), 8 deletions(-)

diff --git a/example.conf b/example.conf
index 9781316..76f3a1d 100644
--- a/example.conf
+++ b/example.conf
@@ -23,6 +23,9 @@ ParticleLoad    = sc       # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc
                            # (increases number of particles by given factor!), 
                            # or 'glass' or 'masked'
 
+## if `ParticleLoad = masked' then you can specify here how masking should take place
+# ParticleMaskType = 3     # bit mask for particle mask (0=center,1=center+edges,2=center+faces,3=center+edges+faces)
+
 ## if `ParticleLoad = glass' then specify here where to load the glass distribution from
 # GlassFileName   = glass128.hdf5
 # GlassTiles      = 1
diff --git a/include/particle_generator.hh b/include/particle_generator.hh
index 37bb946..b0c1a1b 100644
--- a/include/particle_generator.hh
+++ b/include/particle_generator.hh
@@ -38,11 +38,11 @@ namespace particle
         lattice_rsc = 3, // RSC: refined simple cubic
     };
 
-    const std::vector<std::vector<bool>> lattice_masks =
-    {
-        // mask from Richings et al. https://arxiv.org/pdf/2005.14495.pdf
-        {true,true,true,true,true,true,true,false},
-    };
+    // const std::vector<std::vector<bool>> lattice_masks =
+    // {
+    //     // mask from Richings et al. https://arxiv.org/pdf/2005.14495.pdf
+    //     {true,true,true,true,true,true,true,false},
+    // };
 
     const std::vector<std::vector<vec3_t<real_t>>> lattice_shifts =
         {
@@ -266,11 +266,25 @@ namespace particle
                 // set the particle mask 
                 if( cf.get_value_safe("setup","DoBaryons",false) )
                 {
+                    unsigned mask_type = cf.get_value_safe("setup","ParticleMaskType",3);
                     // mask everywehere 0, except the last element
-                    for( auto& m : particle_type_mask_) m = 0;
-                    particle_type_mask_[7] = 1;
+                    for( auto& m : particle_type_mask_) m = -1;
+                    particle_type_mask_[0] = 0; // CDM at corner of unit cube
+                    if( mask_type & 1<<0 ) {
+                        // edge centers
+                        particle_type_mask_[1] = 0; // CDM
+                        particle_type_mask_[2] = 0; // CDM
+                        particle_type_mask_[4] = 0; // CDM
+                    }
+                    if( mask_type & 1<<1 ){
+                        // face centers
+                        particle_type_mask_[3] = 0; // CDM
+                        particle_type_mask_[5] = 0; // CDM
+                        particle_type_mask_[6] = 0; // CDM
+                    }
+                    particle_type_mask_[7] = 1; // baryon at cell center (aka opposite corner)
                 }else{
-                    // mask everywehere 0
+                    // mask everywhere 0, all particle locations are occupied by CDM
                     for( auto& m : particle_type_mask_) m = 0;
                 }