From 3eb3a0cf5d724eda3cef3e88899553e791c550d7 Mon Sep 17 00:00:00 2001 From: Marie Aubert Date: Thu, 11 Oct 2018 17:52:01 +0200 Subject: [PATCH 1/4] useComoving passed to mock_sphere particles + modified lower redshift bound to be consistent with pruneVoids --- c_tools/mock/generateFromCatalog.cpp | 71 ++++++++++++++++++++++------ 1 file changed, 57 insertions(+), 14 deletions(-) diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index 7c046e1..7598bc8 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -104,11 +104,17 @@ double expanFun (double z, void * p) { void loadData(const string& fname, NYU_VData & data) { ifstream f(fname.c_str()); - + int lastidx = -1 ; while (!f.eof()) { NYU_Data d; f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z; + // Maubert - avoid double counting of last particle in array if EOF is after a blank line + if (d.index == lastidx ){ + continue; + } + lastidx = d.index; + // End Maubert d.uniqueID = d.index; data.push_back(d); } @@ -165,7 +171,8 @@ void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data, if (useComoving) { //double pos = gsl_interp_eval(interp, redshifts, dL, data[i].cz, acc); - gsl_integration_qng(&expanF, 1.e-6, data[i].cz/LIGHT_SPEED, + // Maubert - Lower boundary in redshift set to 0 to be consistent with pruneVoids (was 1.e-6 before). + gsl_integration_qng(&expanF, 0.0, data[i].cz/LIGHT_SPEED, 1.e-6, 1.e-6, &result, &error, &nEval); double Dc = result*LIGHT_SPEED; @@ -220,8 +227,24 @@ void generateSurfaceMask(generateFromCatalog_info& args , vector& pixel_list, vector& full_mask_list, NYU_VData& data, - ParticleData& output_data) + ParticleData& output_data, + bool useComoving, + double omegaM) { + + //Maubert - Needed for comobile distance in mock_sphere + gsl_function expanF; + expanF.function = &expanFun; + struct my_expan_params expanParams; + expanParams.Om = omegaM; + expanParams.w0 = -1.0; + expanParams.wa = 0.0; + expanF.params = &expanParams; + + double result, error ; + size_t nEval; + //End Maubert - Needed for comobile distance in mock_sphere + // Find the first free index int idx = -1; int insertion = 0; @@ -367,17 +390,36 @@ void generateSurfaceMask(generateFromCatalog_info& args , // PMS // TEST - insert mock galaxies along spheres of survey redshift boundaries fp = fopen("mock_sphere.txt", "w"); - + //Maubert - insert mock galaxies according to useComoving specification + // Added & compute rmin & rmax out of loop + double rmin ; + double rmax ; + if (useComoving) { + gsl_integration_qng(&expanF, 0.0, args.zMin_arg, + 1.e-6, + 1.e-6, &result, &error, &nEval); + rmin = result* LIGHT_SPEED; + + gsl_integration_qng(&expanF, 0.0, args.zMax_arg, + 1.e-6, + 1.e-6, &result, &error, &nEval); + rmax = result* LIGHT_SPEED; + } else { + rmin = args.zMin_arg * LIGHT_SPEED; + rmax = args.zMax_arg * LIGHT_SPEED; + } + //for (int q = 0; q < 0; q++) { for (int q = 0; q < full_mask_list.size(); q++) { vec3 v = mask.pix2vec(full_mask_list[q]); Position p; - double r = args.zMin_arg * LIGHT_SPEED; - if (r > 0.) { - p.xyz[0] = v.x * r; - p.xyz[1] = v.y * r; - p.xyz[2] = v.z * r; + //double r = args.zMin_arg * LIGHT_SPEED; + + if (rmin > 0.) { + p.xyz[0] = v.x * rmin; + p.xyz[1] = v.y * rmin; + p.xyz[2] = v.z * rmin; output_data.pos.push_back(p); output_data.id_gal.push_back(idx); output_data.ra.push_back(-1); @@ -391,10 +433,10 @@ void generateSurfaceMask(generateFromCatalog_info& args , (p.xyz[2])); } - r = args.zMax_arg * LIGHT_SPEED; - p.xyz[0] = v.x * r; - p.xyz[1] = v.y * r; - p.xyz[2] = v.z * r; + //r = args.zMax_arg * LIGHT_SPEED; + p.xyz[0] = v.x * rmax; + p.xyz[1] = v.y * rmax; + p.xyz[2] = v.z * rmax; output_data.pos.push_back(p); output_data.id_gal.push_back(idx); output_data.ra.push_back(-1); @@ -590,7 +632,8 @@ int main(int argc, char **argv) generateGalaxiesInCube(data, output_data, args_info.useComoving_flag, args_info.omegaM_arg); generateSurfaceMask(args_info, mask, pixel_list, full_mask_list, - data, output_data); + data, output_data,args_info.useComoving_flag, + args_info.omegaM_arg); saveForZobov(output_data, args_info.output_arg, args_info.params_arg); // saveData(output_data); From 7bf85c07c86dd09d612b05c80c87ee5f01a339aa Mon Sep 17 00:00:00 2001 From: Marie Aubert Date: Thu, 11 Oct 2018 18:01:23 +0200 Subject: [PATCH 2/4] useComoving definition moved in the Sample instance to be correctly passed --- pipeline/datasets/example_observation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pipeline/datasets/example_observation.py b/pipeline/datasets/example_observation.py index 091e4ca..b8a7fd6 100644 --- a/pipeline/datasets/example_observation.py +++ b/pipeline/datasets/example_observation.py @@ -49,9 +49,6 @@ figDir = os.getenv("PWD")+"/../figs/example_observation/" ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" -# if true, convert to comoving space using LCDM cosmology -useComoving = True - # optimization: maximum number of parallel threads to use numZobovThreads = 2 @@ -97,6 +94,9 @@ newSample = Sample( # density of mock particles in cubic Mpc/h # (make this as high as you can afford) fakeDensity = 0.05, + + # if true, convert to comoving space using LCDM cosmology + useComoving = True ) dataSampleList.append(newSample) From 531a2c30f3770afc4ae6e833027e1ea89edd0e74 Mon Sep 17 00:00:00 2001 From: Marie Aubert Date: Thu, 11 Oct 2018 18:05:59 +0200 Subject: [PATCH 3/4] Uncommented previous part of the code making sure that upper redshift boundary of zRange is respected --- c_tools/stacking/pruneVoids.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 1a61cc8..f59377b 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -921,7 +921,7 @@ int main(int argc, char **argv) { } voids.resize(iGood); -/* + //Maubert - Uncommented this part : to be sure that voids to not cross maximum redshift asked for in zrange iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { // just in case @@ -933,7 +933,8 @@ int main(int argc, char **argv) { } } voids.resize(iGood); -*/ + // Maubert - End of Uncommented part + printf(" 4th filter: rejected %d outside redshift boundaries\n", numNearZ); From 6a5aa85f2f864f737630e77673a068a314a2aba9 Mon Sep 17 00:00:00 2001 From: Marie Aubert Date: Fri, 12 Oct 2018 11:11:49 +0200 Subject: [PATCH 4/4] Implemented corrections from reviews --- c_tools/mock/generateFromCatalog.cpp | 14 +++++++------- c_tools/stacking/pruneVoids.cpp | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index 7598bc8..d95311d 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -104,16 +104,16 @@ double expanFun (double z, void * p) { void loadData(const string& fname, NYU_VData & data) { ifstream f(fname.c_str()); - int lastidx = -1 ; + while (!f.eof()) { NYU_Data d; f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z; // Maubert - avoid double counting of last particle in array if EOF is after a blank line - if (d.index == lastidx ){ - continue; - } - lastidx = d.index; + if (!f) + { + continue; + } // End Maubert d.uniqueID = d.index; data.push_back(d); @@ -414,7 +414,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , vec3 v = mask.pix2vec(full_mask_list[q]); Position p; - //double r = args.zMin_arg * LIGHT_SPEED; + if (rmin > 0.) { p.xyz[0] = v.x * rmin; @@ -433,7 +433,7 @@ void generateSurfaceMask(generateFromCatalog_info& args , (p.xyz[2])); } - //r = args.zMax_arg * LIGHT_SPEED; + p.xyz[0] = v.x * rmax; p.xyz[1] = v.y * rmax; p.xyz[2] = v.z * rmax; diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index f59377b..ba5cddf 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -921,7 +921,7 @@ int main(int argc, char **argv) { } voids.resize(iGood); - //Maubert - Uncommented this part : to be sure that voids to not cross maximum redshift asked for in zrange + //Maubert - Uncommented this part : to be sure that voids do not cross maximum redshift asked for in zrange iGood = 0; for (iVoid = 0; iVoid < voids.size(); iVoid++) { // just in case