mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-04 23:31:12 +00:00
Slight re-organization of C/C++ tools. Significant modifications to support observational data. Python and pipeline scripts added
This commit is contained in:
P.M. Sutter
parent
15496df4ff
commit
14abbc2018
42 changed files with 16252 additions and 557 deletions
531
c_tools/mock/generateFromCatalog.cpp
Normal file
531
c_tools/mock/generateFromCatalog.cpp
Normal file
|
|
@ -0,0 +1,531 @@
|
|||
#include <healpix_map.h>
|
||||
#include <healpix_map_fitsio.h>
|
||||
#include <boost/format.hpp>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include "generateFromCatalog_conf.h"
|
||||
#include "contour_pixels.hpp"
|
||||
#include <netcdfcpp.h>
|
||||
#include <CosmoTool/fortran.hpp>
|
||||
#include <gsl/gsl_interp.h>
|
||||
#include <gsl/gsl_integration.h>
|
||||
|
||||
#define LIGHT_SPEED 299792.458
|
||||
|
||||
using namespace std;
|
||||
using boost::format;
|
||||
using namespace CosmoTool;
|
||||
|
||||
struct NYU_Data
|
||||
{
|
||||
int index;
|
||||
int sector;
|
||||
int region;
|
||||
double ra, dec;
|
||||
double cz;
|
||||
double fgotten;
|
||||
double phi_z;
|
||||
};
|
||||
|
||||
struct Position
|
||||
{
|
||||
double xyz[3];
|
||||
};
|
||||
|
||||
struct ParticleData
|
||||
{
|
||||
vector<int> id_gal;
|
||||
vector<double> ra;
|
||||
vector<double> dec;
|
||||
vector<double> redshift;
|
||||
int id_mask;
|
||||
// PMS
|
||||
int mask_index;
|
||||
// END PMS
|
||||
vector<Position> pos;
|
||||
double box[3][2];
|
||||
double Lmax;
|
||||
};
|
||||
|
||||
typedef vector<NYU_Data> NYU_VData;
|
||||
|
||||
// this defines the expansion function that we will integrate
|
||||
// Laveaux & Wandelt (2012) Eq. 24
|
||||
struct my_expan_params { double Om; double w0; double wa; };
|
||||
double expanFun (double z, void * p) {
|
||||
struct my_expan_params * params = (struct my_expan_params *)p;
|
||||
double Om = (params->Om);
|
||||
double w0 = (params->w0);
|
||||
double wa = (params->wa);
|
||||
|
||||
//const double h0 = 1.0;
|
||||
const double h0 = 0.71;
|
||||
double ez;
|
||||
|
||||
double wz = w0 + wa*z/(1+z);
|
||||
|
||||
ez = Om*pow(1+z,3) + (1.-Om);
|
||||
//ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz);
|
||||
|
||||
ez = sqrt(ez);
|
||||
//ez = sqrt(ez)/h0;
|
||||
|
||||
ez = 1./ez;
|
||||
|
||||
return ez;
|
||||
}
|
||||
|
||||
void loadData(const string& fname, NYU_VData & data)
|
||||
{
|
||||
ifstream f(fname.c_str());
|
||||
|
||||
while (!f.eof())
|
||||
{
|
||||
NYU_Data d;
|
||||
f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z;
|
||||
data.push_back(d);
|
||||
}
|
||||
}
|
||||
|
||||
void generateGalaxiesInCube(NYU_VData& data, ParticleData& output_data,
|
||||
bool useLCDM)
|
||||
{
|
||||
double d2r = M_PI/180;
|
||||
|
||||
gsl_function expanF;
|
||||
expanF.function = &expanFun;
|
||||
struct my_expan_params expanParams;
|
||||
double maxZ = 2.0, z, result, error, *dL, *redshifts;
|
||||
int numZ = 1000, iZ;
|
||||
size_t nEval;
|
||||
|
||||
expanParams.Om = 0.27;
|
||||
expanParams.w0 = -1.0;
|
||||
expanParams.wa = 0.0;
|
||||
expanF.params = &expanParams;
|
||||
|
||||
dL = (double *) malloc(numZ * sizeof(double));
|
||||
redshifts = (double *) malloc(numZ * sizeof(double));
|
||||
|
||||
for (iZ = 0; iZ < numZ; iZ++) {
|
||||
z = iZ * maxZ/numZ;
|
||||
//gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval);
|
||||
dL[iZ] = result;
|
||||
redshifts[iZ] = z;
|
||||
}
|
||||
|
||||
gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ);
|
||||
gsl_interp_init(interp, redshifts, dL, numZ);
|
||||
gsl_interp_accel *acc = gsl_interp_accel_alloc();
|
||||
|
||||
output_data.pos.resize(data.size());
|
||||
output_data.id_gal.resize(data.size());
|
||||
output_data.ra.resize(data.size());
|
||||
output_data.dec.resize(data.size());
|
||||
output_data.redshift.resize(data.size());
|
||||
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
output_data.box[j][0] = -INFINITY;
|
||||
output_data.box[j][1] = INFINITY;
|
||||
}
|
||||
|
||||
for (int i = 0; i < data.size(); i++)
|
||||
{
|
||||
double ra = data[i].ra*d2r, dec = data[i].dec*d2r;
|
||||
Position& p = output_data.pos[i];
|
||||
|
||||
if (useLCDM) {
|
||||
//double pos = gsl_interp_eval(interp, redshifts, dL, data[i].cz, acc);
|
||||
gsl_integration_qng(&expanF, 1.e-6, data[i].cz/LIGHT_SPEED,
|
||||
1.e-6,
|
||||
1.e-6, &result, &error, &nEval);
|
||||
double Dc = result*LIGHT_SPEED;
|
||||
cout << "HELLO " << data[i].cz << " " << Dc << endl;
|
||||
p.xyz[0] = Dc*cos(ra)*cos(dec);
|
||||
p.xyz[1] = Dc*sin(ra)*cos(dec);
|
||||
p.xyz[2] = Dc*sin(dec);
|
||||
} else {
|
||||
p.xyz[0] = data[i].cz*cos(ra)*cos(dec);
|
||||
p.xyz[1] = data[i].cz*sin(ra)*cos(dec);
|
||||
p.xyz[2] = data[i].cz*sin(dec);
|
||||
}
|
||||
//printf("CREATE %e %e\n", data[i].cz, sqrt(p.xyz[0]*p.xyz[0] + p.xyz[1]*p.xyz[1] + p.xyz[2]*p.xyz[2]));
|
||||
output_data.id_gal[i] = data[i].index;
|
||||
output_data.ra[i] = ra;
|
||||
output_data.dec[i] = dec;
|
||||
output_data.redshift[i] = data[i].cz;
|
||||
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
if (p.xyz[j] > output_data.box[j][0])
|
||||
output_data.box[j][0] = p.xyz[j];
|
||||
if (p.xyz[j] < output_data.box[j][1])
|
||||
output_data.box[j][1] = p.xyz[j];
|
||||
}
|
||||
//printf("INSERT GAL %d %e %e %e\n", output_data.id_gal[i], p.xyz[0], p.xyz[1], p.xyz[2]);
|
||||
}
|
||||
cout << format("Galaxy position generated: %d galaxies") % output_data.pos.size() << endl;
|
||||
cout << format("box is %g < x < %g; %g < y < %g; %g < z < %g")
|
||||
% (1e-2*output_data.box[0][1]) % (1e-2*output_data.box[0][0])
|
||||
% (1e-2*output_data.box[1][1]) % (1e-2*output_data.box[1][0])
|
||||
% (1e-2*output_data.box[2][1]) % (1e-2*output_data.box[2][0]) << endl;
|
||||
|
||||
gsl_interp_free(interp);
|
||||
}
|
||||
|
||||
void generateSurfaceMask(generateFromCatalog_info& args ,
|
||||
Healpix_Map<float>& mask,
|
||||
vector<int>& pixel_list,
|
||||
vector<int>& full_mask_list,
|
||||
NYU_VData& data,
|
||||
ParticleData& output_data)
|
||||
{
|
||||
// Find the first free index
|
||||
int idx = -1;
|
||||
int insertion = 0;
|
||||
double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI;
|
||||
int numToInsert;
|
||||
|
||||
for (int i = 0; i < output_data.id_gal.size(); i++)
|
||||
{
|
||||
if (idx < output_data.id_gal[i])
|
||||
idx = output_data.id_gal[i]+1;
|
||||
}
|
||||
|
||||
output_data.id_mask = idx;
|
||||
|
||||
// PMS
|
||||
output_data.mask_index = output_data.id_gal.size();
|
||||
// END PMS
|
||||
cout << "Generate surface mask..." << endl;
|
||||
|
||||
double Rmax = -1;
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
Rmax = max(Rmax, max(output_data.box[j][0], -output_data.box[j][1]));
|
||||
}
|
||||
|
||||
output_data.Lmax = Rmax;
|
||||
|
||||
// PMS - write a small text file with galaxy position
|
||||
FILE *fp;
|
||||
fp = fopen("galaxies.txt", "w");
|
||||
for (int i = 0; i < data.size(); i++) {
|
||||
Position& p = output_data.pos[i];
|
||||
fprintf(fp, "%e %e %e\n",
|
||||
(p.xyz[0]),
|
||||
(p.xyz[1]),
|
||||
(p.xyz[2]));
|
||||
}
|
||||
fclose(fp);
|
||||
// END PMS
|
||||
|
||||
cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl;
|
||||
volume *= Rmax*Rmax*Rmax/3/1e6;
|
||||
numToInsert = (int)floor(volume*args.density_fake_arg);
|
||||
cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl;
|
||||
|
||||
cout << format("Will insert %d particles") % numToInsert << endl;
|
||||
|
||||
fp = fopen("mock_galaxies.txt", "w");
|
||||
|
||||
double pct = 0;
|
||||
for (int i = 0; i < numToInsert; i++) {
|
||||
double new_pct = i*100./numToInsert;
|
||||
|
||||
if (new_pct-pct > 5.) {
|
||||
pct = new_pct;
|
||||
cout << format(" .. %3.0f %%") % pct << endl;
|
||||
}
|
||||
|
||||
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(),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);
|
||||
|
||||
// PMS : write mock galaxies to a small file for plotting
|
||||
fprintf(fp, "%e %e %e\n",
|
||||
(p.xyz[0]),
|
||||
(p.xyz[1]),
|
||||
(p.xyz[2]));
|
||||
// END PMS
|
||||
output_data.pos.push_back(p);
|
||||
output_data.id_gal.push_back(idx);
|
||||
output_data.ra.push_back(-1);
|
||||
output_data.dec.push_back(-1);
|
||||
output_data.redshift.push_back(-1);
|
||||
//printf("INSERT MOCK %d %e %e %e\n", idx, p.xyz[0], p.xyz[1], p.xyz[2]);
|
||||
insertion++;
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
|
||||
// PMS
|
||||
// TEST - insert mock galaxies along box edge
|
||||
fp = fopen("mock_boundary.txt", "w");
|
||||
double dx[3];
|
||||
dx[0] = output_data.box[0][1] - output_data.box[0][0];
|
||||
dx[1] = output_data.box[1][1] - output_data.box[1][0];
|
||||
dx[2] = output_data.box[2][1] - output_data.box[2][0];
|
||||
|
||||
int nPart = 100;
|
||||
// TEST
|
||||
for (int iDir = 0; iDir < 0; iDir++) {
|
||||
for (int iFace = 0; iFace < 0; iFace++) {
|
||||
//for (int iDir = 0; iDir < 3; iDir++) {
|
||||
//for (int iFace = 0; iFace < 2; iFace++) {
|
||||
|
||||
int iy = (iDir + 1) % 3;
|
||||
int iz = (iDir + 2) % 3;
|
||||
|
||||
for (int i = 0; i < nPart; i++) {
|
||||
for (int j = 0; j < nPart; j++) {
|
||||
Position p;
|
||||
|
||||
p.xyz[iDir] = output_data.box[iDir][iFace];
|
||||
p.xyz[iy] = i * dx[iy]/nPart + output_data.box[iy][0];
|
||||
p.xyz[iz] = j * dx[iz]/nPart + output_data.box[iz][0];
|
||||
|
||||
output_data.pos.push_back(p);
|
||||
output_data.id_gal.push_back(idx);
|
||||
output_data.ra.push_back(-1);
|
||||
output_data.dec.push_back(-1);
|
||||
output_data.redshift.push_back(-1);
|
||||
insertion++;
|
||||
|
||||
fprintf(fp, "%e %e %e\n",
|
||||
(p.xyz[0]),
|
||||
(p.xyz[1]),
|
||||
(p.xyz[2]));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
fclose(fp);
|
||||
// END PMS
|
||||
|
||||
// PMS
|
||||
// TEST - insert mock galaxies along spheres of survey redshift boundaries
|
||||
fp = fopen("mock_sphere.txt", "w");
|
||||
|
||||
for (int p = 0; p < full_mask_list.size(); p++) {
|
||||
vec3 v = mask.pix2vec(full_mask_list[p]);
|
||||
|
||||
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;
|
||||
output_data.pos.push_back(p);
|
||||
output_data.id_gal.push_back(idx);
|
||||
output_data.ra.push_back(-1);
|
||||
output_data.dec.push_back(-1);
|
||||
output_data.redshift.push_back(-1);
|
||||
insertion++;
|
||||
fprintf(fp, "%e %e %e\n",
|
||||
(p.xyz[0]),
|
||||
(p.xyz[1]),
|
||||
(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;
|
||||
output_data.pos.push_back(p);
|
||||
output_data.id_gal.push_back(idx);
|
||||
output_data.ra.push_back(-1);
|
||||
output_data.dec.push_back(-1);
|
||||
output_data.redshift.push_back(-1);
|
||||
insertion++;
|
||||
fprintf(fp, "%e %e %e\n",
|
||||
(p.xyz[0]),
|
||||
(p.xyz[1]),
|
||||
(p.xyz[2]));
|
||||
}
|
||||
fclose(fp);
|
||||
// END PMS
|
||||
|
||||
cout << format("Done. Inserted %d particles.") % insertion << endl;
|
||||
}
|
||||
|
||||
void saveData(ParticleData& pdata)
|
||||
{
|
||||
NcFile f("particles.nc", NcFile::Replace);
|
||||
|
||||
assert(f.is_valid());
|
||||
|
||||
NcDim *d = f.add_dim("space", 3);
|
||||
NcDim *p = f.add_dim("Np", pdata.pos.size());
|
||||
NcVar *v = f.add_var("particles", ncDouble, d, p);
|
||||
double *x = new double[pdata.pos.size()];
|
||||
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
|
||||
for (int i = 0; i < pdata.pos.size(); i++)
|
||||
x[i] = pdata.pos[i].xyz[j];
|
||||
|
||||
v->put_rec(d, x, j);
|
||||
}
|
||||
|
||||
v = f.add_var("id_gal", ncInt, p);
|
||||
v->put(&pdata.id_gal[0], pdata.id_gal.size());
|
||||
|
||||
delete[] x;
|
||||
|
||||
}
|
||||
|
||||
void saveForZobov(ParticleData& pdata, const string& fname, const string& paramname)
|
||||
{
|
||||
UnformattedWrite f(fname);
|
||||
static const char axis[] = { 'X', 'Y', 'Z' };
|
||||
double Lmax = pdata.Lmax;
|
||||
|
||||
f.beginCheckpoint();
|
||||
f.writeInt32(pdata.pos.size());
|
||||
f.endCheckpoint();
|
||||
|
||||
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
cout << format("Writing %c components...") % axis[j] << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < pdata.pos.size(); i++)
|
||||
{
|
||||
f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax));
|
||||
}
|
||||
f.endCheckpoint();
|
||||
}
|
||||
|
||||
NcFile fp(paramname.c_str(), NcFile::Replace);
|
||||
|
||||
fp.add_att("range_x_min", -Lmax/100.);
|
||||
fp.add_att("range_x_max", Lmax/100.);
|
||||
fp.add_att("range_y_min", -Lmax/100.);
|
||||
fp.add_att("range_y_max", Lmax/100.);
|
||||
fp.add_att("range_z_min", -Lmax/100.);
|
||||
fp.add_att("range_z_max", Lmax/100.);
|
||||
fp.add_att("mask_index", pdata.mask_index); // PMS
|
||||
|
||||
int nOutputPart = pdata.mask_index;
|
||||
//int nOutputPart = pdata.pos.size();
|
||||
|
||||
NcDim *NumPart_dim = fp.add_dim("numpart_dim", nOutputPart);
|
||||
NcVar *v = fp.add_var("particle_ids", ncInt, NumPart_dim);
|
||||
NcVar *vra = fp.add_var("RA", ncInt, NumPart_dim);
|
||||
NcVar *vdec = fp.add_var("DEC", ncInt, NumPart_dim);
|
||||
NcVar *vredshift = fp.add_var("z", ncInt, NumPart_dim);
|
||||
//NcVar *v2 = fp.add_var("expansion", ncDouble, NumPart_dim);
|
||||
|
||||
//double *expansion_fac = new double[pdata.pos.size()];
|
||||
|
||||
//for (int i = 0; i < pdata.pos.size(); i++)
|
||||
// expansion_fac[i] = 1.0;
|
||||
|
||||
v->put(&pdata.id_gal[0], nOutputPart);
|
||||
vra->put(&pdata.ra[0], nOutputPart);
|
||||
vdec->put(&pdata.dec[0], nOutputPart);
|
||||
vredshift->put(&pdata.redshift[0], nOutputPart);
|
||||
//v2->put(expansion_fac, pdata.pos.size());
|
||||
|
||||
//delete[] expansion_fac;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
generateFromCatalog_info args_info;
|
||||
generateFromCatalog_conf_params args_params;
|
||||
|
||||
generateFromCatalog_conf_init(&args_info);
|
||||
generateFromCatalog_conf_params_init(&args_params);
|
||||
|
||||
args_params.check_required = 0;
|
||||
if (generateFromCatalog_conf_ext (argc, argv, &args_info, &args_params))
|
||||
return 1;
|
||||
|
||||
if (!args_info.configFile_given)
|
||||
{
|
||||
if (generateFromCatalog_conf_required (&args_info, GENERATEFROMCATALOG_CONF_PACKAGE))
|
||||
return 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
args_params.check_required = 1;
|
||||
args_params.initialize = 0;
|
||||
if (generateFromCatalog_conf_config_file (args_info.configFile_arg,
|
||||
&args_info,
|
||||
&args_params))
|
||||
return 1;
|
||||
}
|
||||
|
||||
generateFromCatalog_conf_print_version();
|
||||
|
||||
cout << "Loading NYU data..." << endl;
|
||||
vector<NYU_Data> data;
|
||||
Healpix_Map<float> o_mask;
|
||||
vector<int> pixel_list;
|
||||
vector<int> full_mask_list;
|
||||
ParticleData output_data;
|
||||
|
||||
loadData(args_info.catalog_arg, data);
|
||||
|
||||
|
||||
cout << "Loading mask..." << endl;
|
||||
read_Healpix_map_from_fits(args_info.mask_arg, o_mask);
|
||||
|
||||
Healpix_Map<float> mask;
|
||||
|
||||
mask.SetNside(128, RING);
|
||||
mask.Import(o_mask);
|
||||
|
||||
computeContourPixels(mask,pixel_list);
|
||||
|
||||
// We compute a cube holding all the galaxies + the survey surface mask
|
||||
|
||||
cout << "Placing galaxies..." << endl;
|
||||
generateGalaxiesInCube(data, output_data, args_info.useLCDM_flag);
|
||||
generateSurfaceMask(args_info, mask, pixel_list, full_mask_list,
|
||||
data, output_data);
|
||||
|
||||
saveForZobov(output_data, args_info.output_arg, args_info.params_arg);
|
||||
// saveData(output_data);
|
||||
|
||||
// PMS
|
||||
FILE *fp = fopen("mask_index.txt", "w");
|
||||
fprintf(fp, "%d", output_data.mask_index);
|
||||
fclose(fp);
|
||||
|
||||
fp = fopen("sample_info.txt", "w");
|
||||
fprintf(fp, "Lmax = %f\n", output_data.Lmax);
|
||||
fprintf(fp, "mask_index = %d\n", output_data.mask_index);
|
||||
fprintf(fp, "total_particles = %d\n", output_data.pos.size());
|
||||
fclose(fp);
|
||||
|
||||
fp = fopen("total_particles.txt", "w");
|
||||
fprintf(fp, "%d", output_data.pos.size());
|
||||
fclose(fp);
|
||||
printf("Done!\n");
|
||||
// END PMS
|
||||
return 0;
|
||||
}
|
||||
17
c_tools/mock/generateFromCatalog.ggo
Normal file
17
c_tools/mock/generateFromCatalog.ggo
Normal file
|
|
@ -0,0 +1,17 @@
|
|||
package "generateFromCatalog"
|
||||
version "alpha"
|
||||
|
||||
option "configFile" - "Configuration filename" string optional
|
||||
|
||||
option "catalog" - "Input NYU-VAGC catalog" string required
|
||||
option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" string required
|
||||
option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)" double optional default="1"
|
||||
|
||||
option "zMin" - "Minimum redshift of data" double required
|
||||
|
||||
option "zMax" - "Maximum redshift of data" double required
|
||||
|
||||
option "output" - "Filename of particle datafile" string required
|
||||
option "params" - "Output parameters of the datacube" string required
|
||||
|
||||
option "useLCDM" - "Convert to real space using LCDM cosmology" flag off
|
||||
508
c_tools/mock/generateMock.cpp
Normal file
508
c_tools/mock/generateMock.cpp
Normal file
|
|
@ -0,0 +1,508 @@
|
|||
#include <cmath>
|
||||
#include <cassert>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <CosmoTool/loadSimu.hpp>
|
||||
#include <CosmoTool/loadRamses.hpp>
|
||||
#include <CosmoTool/loadGadget.hpp>
|
||||
#include <CosmoTool/loadFlash.hpp>
|
||||
#include <CosmoTool/interpolate.hpp>
|
||||
#include <CosmoTool/fortran.hpp>
|
||||
#include "generateMock_conf.h"
|
||||
#include "gslIntegrate.hpp"
|
||||
#include <netcdfcpp.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace CosmoTool;
|
||||
|
||||
#define LIGHT_SPEED 299792.458
|
||||
|
||||
SimuData *myLoadGadget(const char *fname, int id, int flags)
|
||||
{
|
||||
return loadGadgetMulti(fname, id, flags);
|
||||
}
|
||||
|
||||
SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, SimuData *(*loadFunction)(const char *fname, int id, int flags))
|
||||
{
|
||||
SimuData *d, *outd;
|
||||
bool singleFile = false;
|
||||
|
||||
try
|
||||
{
|
||||
d = loadFunction(gadgetname, -1, 0);
|
||||
singleFile = true;
|
||||
}
|
||||
catch (const NoSuchFileException& e)
|
||||
{
|
||||
try
|
||||
{
|
||||
d = loadFunction(gadgetname, 0, 0);
|
||||
}
|
||||
catch(const NoSuchFileException& e)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
outd = new SimuData;
|
||||
|
||||
outd->NumPart = d->TotalNumPart;
|
||||
outd->BoxSize = d->BoxSize/1000;
|
||||
outd->TotalNumPart = outd->NumPart;
|
||||
outd->Hubble = d->Hubble;
|
||||
outd->Omega_Lambda = d->Omega_Lambda;
|
||||
outd->Omega_M = d->Omega_M;
|
||||
outd->time = d->time;
|
||||
|
||||
for (int k = 0; k < 3; k++)
|
||||
outd->Pos[k] = new float[outd->NumPart];
|
||||
outd->Vel[2] = new float[outd->NumPart];
|
||||
delete d;
|
||||
|
||||
int curCpu = singleFile ? -1 : 0;
|
||||
cout << "loading file 0 " << endl;
|
||||
try
|
||||
{
|
||||
while (1)
|
||||
{
|
||||
d = loadFunction(gadgetname, curCpu, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID);
|
||||
for (int k = 0; k < 3; k++)
|
||||
for (int i = 0; i < d->NumPart; i++)
|
||||
{
|
||||
assert(d->Id[i] >= 1);
|
||||
assert(d->Id[i] <= outd->TotalNumPart);
|
||||
outd->Pos[k][d->Id[i]-1] = d->Pos[k][i]/1000;
|
||||
outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i];
|
||||
}
|
||||
|
||||
if (goRedshift)
|
||||
for (int i = 0; i < d->NumPart; i++)
|
||||
outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.;
|
||||
|
||||
delete d;
|
||||
if (singleFile)
|
||||
break;
|
||||
curCpu++;
|
||||
cout << "loading file " << curCpu << endl;
|
||||
}
|
||||
}
|
||||
catch (const NoSuchFileException& e)
|
||||
{
|
||||
}
|
||||
|
||||
return outd;
|
||||
}
|
||||
|
||||
|
||||
SimuData *doLoadMultidark(const char *multidarkname)
|
||||
{
|
||||
SimuData *outd;
|
||||
FILE *fp;
|
||||
int actualNumPart;
|
||||
|
||||
outd = new SimuData;
|
||||
cout << "opening multidark file " << multidarkname << endl;
|
||||
fp = fopen(multidarkname, "r");
|
||||
if (fp == NULL) {
|
||||
cout << "could not open file!" << endl;
|
||||
return 0;
|
||||
}
|
||||
fscanf(fp, "%f\n", &outd->BoxSize);
|
||||
fscanf(fp, "%f\n", &outd->Omega_M);
|
||||
fscanf(fp, "%f\n", &outd->Hubble);
|
||||
fscanf(fp, "%f\n", &outd->time);
|
||||
fscanf(fp, "%d\n", &outd->NumPart);
|
||||
|
||||
outd->time = 1./(1.+outd->time); // convert to scale factor
|
||||
outd->TotalNumPart = outd->NumPart;
|
||||
outd->Omega_Lambda = 1.0 - outd->Omega_M;
|
||||
|
||||
for (int k = 0; k < 3; k++)
|
||||
outd->Pos[k] = new float[outd->NumPart];
|
||||
outd->Vel[2] = new float[outd->NumPart];
|
||||
|
||||
cout << "loading multidark particles" << endl;
|
||||
actualNumPart = 0;
|
||||
for (int i = 0; i < outd->NumPart; i++) {
|
||||
fscanf(fp, "%f %f %f %f\n", &outd->Pos[0][i], &outd->Pos[1][i],
|
||||
&outd->Pos[2][i], &outd->Vel[2][i]);
|
||||
|
||||
if (outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 &&
|
||||
outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) {
|
||||
break;
|
||||
} else {
|
||||
actualNumPart++;
|
||||
}
|
||||
}
|
||||
fclose(fp);
|
||||
|
||||
outd->NumPart = actualNumPart;
|
||||
outd->TotalNumPart = actualNumPart;
|
||||
return outd;
|
||||
}
|
||||
|
||||
|
||||
static double cubic(double a)
|
||||
{
|
||||
return a*a*a;
|
||||
}
|
||||
|
||||
struct TotalExpansion
|
||||
{
|
||||
double Omega_M, Omega_L;
|
||||
|
||||
double operator()(double z)
|
||||
{
|
||||
return 1/sqrt(Omega_M*cubic(1+z) + Omega_L);
|
||||
}
|
||||
};
|
||||
|
||||
Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z1, int N = 5000)
|
||||
{
|
||||
TotalExpansion e_computer;
|
||||
double D_tilde, Q, Qprime;
|
||||
InterpolatePairs pairs;
|
||||
|
||||
e_computer.Omega_M = OM;
|
||||
e_computer.Omega_L = OL;
|
||||
|
||||
pairs.resize(N);
|
||||
ofstream f("comoving_distance.txt");
|
||||
|
||||
for (int i = 0; i < N; i++)
|
||||
{
|
||||
double z = z0 + (z1-z0)/N*i;
|
||||
|
||||
pairs[i].second = z;
|
||||
pairs[i].first = gslIntegrate(e_computer, 0, z, 1e-3);
|
||||
f << z << " " << pairs[i].first << endl;
|
||||
}
|
||||
|
||||
return buildFromVector(pairs);
|
||||
}
|
||||
|
||||
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact)
|
||||
{
|
||||
int x0, x1, x2;
|
||||
|
||||
switch (axis) {
|
||||
case 0:
|
||||
x0 = 1; x1 = 2; x2 = 0;
|
||||
break;
|
||||
case 1:
|
||||
x0 = 0; x1 = 2; x2 = 1;
|
||||
break;
|
||||
case 2:
|
||||
x0 = 0; x1 = 1; x2 = 2;
|
||||
break;
|
||||
default:
|
||||
abort();
|
||||
}
|
||||
|
||||
double z0 = 1/data->time - 1;
|
||||
Interpolate z_vs_D = make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 0., z0+8*data->BoxSize*100/LIGHT_SPEED); // Redshift 2*z0 should be sufficient ?
|
||||
double z_base = reshift ? z0 : 0;
|
||||
TotalExpansion e_computer;
|
||||
double baseComovingDistance;
|
||||
|
||||
expfact = new double[data->NumPart];
|
||||
|
||||
cout << "Using base redshift z=" << z0 << " " << z0+8*data->BoxSize*100/LIGHT_SPEED << endl;
|
||||
|
||||
e_computer.Omega_M = data->Omega_M;
|
||||
e_computer.Omega_L = data->Omega_Lambda;
|
||||
baseComovingDistance = LIGHT_SPEED/100.* gslIntegrate(e_computer, 0, z0, 1e-3);
|
||||
cout << "Comoving distance = " << baseComovingDistance << " Mpc/h" << endl;
|
||||
|
||||
for (uint32_t i = 0; i < data->NumPart; i++)
|
||||
{
|
||||
float& x = data->Pos[x0][i];
|
||||
float& y = data->Pos[x1][i];
|
||||
float& z = data->Pos[x2][i];
|
||||
float& v = data->Vel[2][i];
|
||||
float z_old = z;
|
||||
|
||||
double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED;
|
||||
try
|
||||
{
|
||||
|
||||
// Distorted redshift
|
||||
if (reduced_red == 0)
|
||||
z = 0;
|
||||
else {
|
||||
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
|
||||
}
|
||||
expfact[i] = z / z_old;
|
||||
// Add peculiar velocity
|
||||
if (pecvel)
|
||||
z += v/100;
|
||||
}
|
||||
catch(const InvalidRangeException& e) {
|
||||
cout << "Trying to interpolate out of the tabulated range." << endl;
|
||||
cout << "The offending value is z=" << reduced_red << endl;
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void generateOutput(SimuData *data, int axis,
|
||||
const std::string& fname)
|
||||
{
|
||||
UnformattedWrite f(fname);
|
||||
|
||||
cout << "Generating output particles to " << fname << endl;
|
||||
|
||||
int x0, x1, x2;
|
||||
|
||||
switch (axis) {
|
||||
case 0:
|
||||
x0 = 1; x1 = 2; x2 = 0;
|
||||
break;
|
||||
case 1:
|
||||
x0 = 0; x1 = 2; x2 = 1;
|
||||
break;
|
||||
case 2:
|
||||
x0 = 0; x1 = 1; x2 = 2;
|
||||
break;
|
||||
default:
|
||||
abort();
|
||||
}
|
||||
|
||||
f.beginCheckpoint();
|
||||
f.writeInt32(data->NumPart);
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing X components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < data->NumPart; i++)
|
||||
{
|
||||
f.writeReal32(data->Pos[x0][i]);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing Y components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < data->NumPart; i++)
|
||||
{
|
||||
f.writeReal32(data->Pos[x1][i]);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing Z components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < data->NumPart; i++)
|
||||
{
|
||||
f.writeReal32(data->Pos[x2][i]);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
}
|
||||
|
||||
void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info)
|
||||
{
|
||||
float subsample = args_info.subsample_given ? args_info.subsample_arg : 1.0;
|
||||
uint32_t goodParticles = 0;
|
||||
double ranges[3][2] = {
|
||||
{ args_info.rangeX_min_arg, args_info.rangeX_max_arg },
|
||||
{ args_info.rangeY_min_arg, args_info.rangeY_max_arg },
|
||||
{ args_info.rangeZ_min_arg, args_info.rangeZ_max_arg }
|
||||
};
|
||||
double mul[3];
|
||||
float minmax[2][3];
|
||||
int *particle_id;
|
||||
bool *random_acceptance = 0;
|
||||
|
||||
boxed = new SimuData;
|
||||
boxed->Hubble = simu->Hubble;
|
||||
boxed->Omega_M = simu->Omega_M;
|
||||
boxed->Omega_Lambda = simu->Omega_Lambda;
|
||||
boxed->time = simu->time;
|
||||
boxed->BoxSize = simu->BoxSize;
|
||||
|
||||
random_acceptance = new bool[simu->NumPart];
|
||||
|
||||
for (int j = 0; j < 3; j++) minmax[1][j] = minmax[0][j] = simu->Pos[j][0];
|
||||
|
||||
for (uint32_t i = 0; i < simu->NumPart; i++)
|
||||
{
|
||||
bool acceptance = true;
|
||||
|
||||
for (int j = 0; j < 3; j++) {
|
||||
acceptance =
|
||||
acceptance &&
|
||||
(simu->Pos[j][i] > ranges[j][0]) &&
|
||||
(simu->Pos[j][i] < ranges[j][1]);
|
||||
minmax[0][j] = min(simu->Pos[j][i], minmax[0][j]);
|
||||
minmax[1][j] = max(simu->Pos[j][i], minmax[1][j]);
|
||||
}
|
||||
random_acceptance[i] = acceptance && (drand48() <= subsample);
|
||||
if (random_acceptance[i])
|
||||
goodParticles++;
|
||||
}
|
||||
|
||||
cout << "Subsample fraction: " << subsample << endl;
|
||||
cout << "Min range = " << ranges[0][0] << " " << ranges[1][0] << " " << ranges[2][0] << endl;
|
||||
cout << "Max range = " << ranges[0][1] << " " << ranges[1][1] << " " << ranges[2][1] << endl;
|
||||
|
||||
cout << "Min position = " << minmax[0][0] << " " << minmax[0][1] << " " << minmax[0][2] << endl;
|
||||
cout << "Max position = " << minmax[1][0] << " " << minmax[1][1] << " " << minmax[1][2] << endl;
|
||||
|
||||
cout << "Number of accepted particles: " << goodParticles << endl;
|
||||
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
boxed->Pos[j] = new float[goodParticles];
|
||||
boxed->Vel[j] = 0;
|
||||
mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
|
||||
}
|
||||
|
||||
cout << "Rescaling factors = " << mul[0] << " " << mul[1] << " " << mul[2] << endl;
|
||||
boxed->NumPart = goodParticles;
|
||||
|
||||
particle_id = new int[goodParticles];
|
||||
double *expansion_fac = new double[goodParticles];
|
||||
|
||||
uint32_t k = 0;
|
||||
for (uint32_t i = 0; i < simu->NumPart; i++)
|
||||
{
|
||||
bool acceptance = random_acceptance[i];
|
||||
|
||||
if (acceptance)
|
||||
{
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
boxed->Pos[j][k] = (simu->Pos[j][i]-ranges[j][0])*mul[j];
|
||||
assert(boxed->Pos[j][k] > 0);
|
||||
assert(boxed->Pos[j][k] < 1);
|
||||
}
|
||||
particle_id[k] = i;
|
||||
expansion_fac[k] = efac[i];
|
||||
k++;
|
||||
}
|
||||
}
|
||||
|
||||
delete[] random_acceptance;
|
||||
|
||||
NcFile f(args_info.outputParameter_arg, NcFile::Replace);
|
||||
|
||||
f.add_att("range_x_min", ranges[0][0]);
|
||||
f.add_att("range_x_max", ranges[0][1]);
|
||||
f.add_att("range_y_min", ranges[1][0]);
|
||||
f.add_att("range_y_max", ranges[1][1]);
|
||||
f.add_att("range_z_min", ranges[2][0]);
|
||||
f.add_att("range_z_max", ranges[2][1]);
|
||||
|
||||
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
|
||||
NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);
|
||||
NcVar *v2 = f.add_var("expansion", ncDouble, NumPart_dim);
|
||||
|
||||
v->put(particle_id, boxed->NumPart);
|
||||
v2->put(expansion_fac, boxed->NumPart);
|
||||
|
||||
delete[] particle_id;
|
||||
delete[] expansion_fac;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
generateMock_info args_info;
|
||||
generateMock_conf_params args_params;
|
||||
SimuData *simu, *simuOut;
|
||||
|
||||
generateMock_conf_init(&args_info);
|
||||
generateMock_conf_params_init(&args_params);
|
||||
|
||||
args_params.check_required = 0;
|
||||
if (generateMock_conf_ext (argc, argv, &args_info, &args_params))
|
||||
return 1;
|
||||
|
||||
if (!args_info.configFile_given)
|
||||
{
|
||||
if (generateMock_conf_required (&args_info, GENERATEMOCK_CONF_PACKAGE))
|
||||
return 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
args_params.check_required = 1;
|
||||
args_params.initialize = 0;
|
||||
if (generateMock_conf_config_file (args_info.configFile_arg,
|
||||
&args_info,
|
||||
&args_params))
|
||||
return 1;
|
||||
}
|
||||
|
||||
generateMock_conf_print_version();
|
||||
|
||||
if (args_info.ramsesBase_given || args_info.ramsesId_given)
|
||||
{
|
||||
if (args_info.ramsesBase_given && args_info.ramsesId_given) {
|
||||
//simu = doLoadRamses(args_info.ramsesBase_arg,
|
||||
// args_info.ramsesId_arg,
|
||||
// args_info.axis_arg, false);
|
||||
}
|
||||
else
|
||||
{
|
||||
cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl;
|
||||
return 1;
|
||||
}
|
||||
|
||||
if (simu == 0)
|
||||
{
|
||||
cerr << "Error while loading" << endl;
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
else if (args_info.gadget_given || args_info.flash_given || args_info.multidark_given)
|
||||
{
|
||||
if (args_info.gadget_given && args_info.flash_given)
|
||||
{
|
||||
cerr << "Do not know which file to use: Gadget or Flash ?" << endl;
|
||||
return 1;
|
||||
}
|
||||
|
||||
if (args_info.multidark_given) {
|
||||
simu = doLoadMultidark(args_info.multidark_arg);
|
||||
}
|
||||
|
||||
if (args_info.gadget_given) {
|
||||
simu = doLoadSimulation(args_info.gadget_arg, args_info.axis_arg, false, myLoadGadget);
|
||||
}
|
||||
if (args_info.flash_given) {
|
||||
simu = doLoadSimulation(args_info.flash_arg, args_info.axis_arg, false, loadFlashMulti);
|
||||
}
|
||||
|
||||
if (simu == 0)
|
||||
{
|
||||
cerr << "Error while loading " << endl;
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
cerr << "Either a ramses snapshot or a gadget snapshot is required." << endl;
|
||||
return 1;
|
||||
}
|
||||
cout << "Hubble = " << simu->Hubble << endl;
|
||||
cout << "Boxsize = " << simu->BoxSize << endl;
|
||||
cout << "Omega_M = " << simu->Omega_M << endl;
|
||||
cout << "Omega_Lambda = " << simu->Omega_Lambda << endl;
|
||||
|
||||
double *expfact;
|
||||
|
||||
if (args_info.cosmo_flag){
|
||||
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag, expfact);
|
||||
} else {
|
||||
expfact = new double[simu->NumPart];
|
||||
for (int j = 0; j < simu->NumPart; j++) expfact[j] = 1.0;
|
||||
}
|
||||
|
||||
makeBox(simu, expfact, simuOut, args_info);
|
||||
delete simu;
|
||||
|
||||
generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);
|
||||
|
||||
delete simuOut;
|
||||
|
||||
printf("Done!\n");
|
||||
return 0;
|
||||
}
|
||||
31
c_tools/mock/generateMock.ggo
Normal file
31
c_tools/mock/generateMock.ggo
Normal file
|
|
@ -0,0 +1,31 @@
|
|||
package "generateMock"
|
||||
version "0"
|
||||
|
||||
option "configFile" - "Configuration file path" string optional
|
||||
|
||||
# Ramses data
|
||||
option "ramsesBase" - "Base directory for ramses" string optional
|
||||
option "ramsesId" - "Ramses snapshot id" int optional
|
||||
|
||||
option "gadget" - "Base name of gadget snapshot (without parallel writing extension)" string optional
|
||||
option "flash" - "Base name for FLASH snapshot" string optional
|
||||
option "multidark" - "Base name for multidark snapshot" string optional
|
||||
|
||||
option "axis" - "Redshift axis (X=0, Y=1, Z=2)" int optional default="2"
|
||||
|
||||
option "output" - "Output filename for particles" string required
|
||||
option "outputParameter" - "Output geometry parameter file for postprocessing" string required
|
||||
|
||||
option "rangeX_min" - "Minimum range in X for making the box" double required
|
||||
option "rangeX_max" - "Maximum range in X for making the box" double required
|
||||
option "rangeY_min" - "Minimum range in Y for making the box" double required
|
||||
option "rangeY_max" - "Maximum range in Y for making the box" double required
|
||||
option "rangeZ_min" - "Minimum range in Z for making the box (after distortion)" double required
|
||||
option "rangeZ_max" - "Maximum range in Z for making the box (after distortion)" double required
|
||||
|
||||
option "preReShift" - "Reshift the zero of the Z axis" flag off
|
||||
option "peculiarVelocities" - "Added peculiar velocities distortion" flag off
|
||||
|
||||
option "cosmo" - "Apply cosmological redshift" flag on
|
||||
|
||||
option "subsample" - "Subsample the input simulation by the specified amount" double optional
|
||||
47
c_tools/mock/generateTestMock.cpp
Normal file
47
c_tools/mock/generateTestMock.cpp
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
#include <CosmoTool/fortran.hpp>
|
||||
|
||||
using namespace CosmoTool;
|
||||
using namespace std;
|
||||
|
||||
#define LX 1.0
|
||||
#define LY 1.0
|
||||
#define LZ 1.0
|
||||
|
||||
#define NUMPART (16*16*16)
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
UnformattedWrite f("particles.bin");
|
||||
|
||||
f.beginCheckpoint();
|
||||
f.writeInt32(NUMPART);
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing X components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < NUMPART; i++)
|
||||
{
|
||||
f.writeReal32(drand48()*LX);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing Y components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < NUMPART; i++)
|
||||
{
|
||||
f.writeReal32(drand48()*LY);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
|
||||
cout << "Writing Z components..." << endl;
|
||||
f.beginCheckpoint();
|
||||
for (uint32_t i = 0; i < NUMPART; i++)
|
||||
{
|
||||
f.writeReal32(drand48()*LZ);
|
||||
}
|
||||
f.endCheckpoint();
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue