Add more mask

mytools/contour_pixels.cpp

@@ -7,8 +7,17 @@ using namespace std;
 
 static const bool DEBUG = true;
 
+void computeFilledPixels(Healpix_Map<float>& m, vector<int>& filled)
+{
+  filled.clear();
+  for (int p = 0; p < m.Npix(); p++)
+    if (m[p] > 0)
+      filled.push_back(p);
+}
+
 void computeContourPixels(Healpix_Map<float>& m, vector<int>& contour)
 {
+  contour.clear();
   for (int p = 0; p < m.Npix(); p++)
     {
       fix_arr<int, 8> result;

mytools/contour_pixels.hpp

@@ -4,6 +4,8 @@
 #include <vector>
 #include <healpix_map.h>
 
+
 void computeContourPixels(Healpix_Map<float>& map, std::vector<int>& contour);
+void computeFilledPixels(Healpix_Map<float>& map, std::vector<int>& contour);
 
 #endif

mytools/generateFromCatalog.cpp

@@ -92,6 +92,65 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data)
   
 }
 
+static double cube(double x)
+{
+  return x*x*x;
+}
+
+void generateBoxMask(generateFromCatalog_info& args ,
+                         Healpix_Map<float>& mask,
+                         vector<int>& pixel_list,
+                         NYU_VData& data,
+                         ParticleData& output_data)
+{
+  int idx = -1;
+  int insertion = 0;
+  double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI;
+  int numToInsert;
+
+  idx = output_data.id_mask;
+
+  cout << "Generate box mask..." << endl;
+
+  double Rmax = output_data.Lmax, t = args.box_thickness_arg;
+  output_data.Lmax += args.box_thickness_arg;
+
+
+  volume *= Rmax*Rmax/1e6 * args.box_thickness_arg;
+  numToInsert = (int)floor(volume*args.box_density_arg);
+
+  for (int i = 0; i < numToInsert; i++)
+    {
+      Position p;
+      bool stop_here;
+
+      do
+        {
+
+          int p0 = (int)floor(drand48()*pixel_list.size());
+          vec3 v = mask.pix2vec(pixel_list[p0]);
+          double r = Rmax*pow(drand48()*(cube(1+t/Rmax)-1) + 1,1./3);
+
+          p.xyz[0] = v.x * r;
+          p.xyz[1] = v.y * r;
+          p.xyz[2] = v.z * r;
+
+          stop_here = true;
+          for (int j = 0; j < 3; j++)
+            {
+              if (p.xyz[j] > output_data.box[j][0] ||
+                  p.xyz[j] < output_data.box[j][1])
+                stop_here = false;
+            }
+        }
+      while (!stop_here);
+
+      output_data.pos.push_back(p);
+      output_data.id_gal.push_back(idx);
+      insertion++;
+    }
+}
+
 void generateSurfaceMask(generateFromCatalog_info& args ,
 			 Healpix_Map<float>& mask, 
 			 vector<int>& pixel_list, 
@@ -296,6 +355,8 @@ int main(int argc, char **argv)
 
   generateGalaxiesInCube(data, output_data);
   generateSurfaceMask(args_info, mask, pixel_list, data, output_data);
+  computeFilledPixels(mask,pixel_list);
+  generateBoxMask(args_info, mask, pixel_list, data, output_data);
   
   saveForZobov(output_data, args_info.output_arg, args_info.params_arg);
   //  saveData(output_data);

mytools/generateFromCatalog.ggo

@@ -8,4 +8,7 @@ option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" st
 option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)"  double optional default="1"
 
 option "output"	- "Filename of particle datafile" string required
-option "params" - "Output parameters of the datacube" string required
+option "params" - "Output parameters of the datacube" string required
+
+option "box_thickness"  - "Thickness of the limiting mask" double optional default="1"
+option  "box_density" - "Fake Galaxy density on the box" double optional default="1"

voztie.c

@@ -94,9 +94,9 @@ int main(int argc, char *argv[]) {
 	fread(&volstemp,1,sizeof(float),part);
 	if (vols[orig[p]] > -1.)
 	  if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) {
-	    printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
-		   orig[p],vols[orig[p]],volstemp);
-	    exit(0);
+//	    printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
+//		   orig[p],vols[orig[p]],volstemp);
+//	    exit(0);
 	  }
 	vols[orig[p]] = volstemp;
       }