bccp / nbodykit

Analysis kit for large-scale structure datasets, the massively parallel way
http://nbodykit.rtfd.io
GNU General Public License v3.0
111 stars 60 forks source link

Data points being painted to mesh in a funny way. #659

Open danpryer opened 3 years ago

danpryer commented 3 years ago

Hey,

I'm trying to do a convolved FFT power spectrum measurement on some simulated galaxy data I have. The data is distributed in -x_min < x < x_max (i.e. the cartesian coordinates can be negative, but i'm still having an issue if i shift all of the coordinates so that they are in the range 0 < x < x_max instead). The data is also not isotropic (there is an angular survey mask that has already been applied to the data).

I also have a random catalogue that i created myself, which matches the data very closely.

The issue I have is that when i paint the data and randoms to grid in nbodykit, it seems like some sort of coordinate shift is done on them as the geometry is totally distorted when i do a "plt.imshow(mesh.preview(axes=[0,1]))":

Screenshot 2021-09-13 at 11 16 25

Instead of seeing circular arks around the centre of the box where the observer should be, you can see the arcs have been broken up, and flipped such that they arc around the corners of the box.

The resulting power spectrum is, as you can imagine, totally wrong.

I do not specify the box size or the box centre, but you can see from the mesh.attrs that it seems to be detecting these correctly (i've calculated them myself, and they match what i think they should be).

When i paint the data and randoms to their own grids using my own custom algorithms they look like this (what i would expect):

Screenshot 2021-09-13 at 11 40 51

Any pointers on where i could be going wrong?

Heres the main bits of the code:

data = fitsio.read('test_data.fits') rand = fitsio.read('test_rand.fits')

cat_data = ArrayCatalog(data) cat_rand = ArrayCatalog(rand)

cat_data['Position'] = cat_data['xpos'][:, None] [1, 0, 0] + cat_data['ypos'][:, None] [0, 1, 0] + cat_data['zpos'][:, None] [0, 0, 1] cat_rand['Position'] = cat_rand['xpos'][:, None] [1, 0, 0] + cat_rand['ypos'][:, None] [0, 1, 0] + cat_rand['zpos'][:, None] [0, 0, 1]

FSKY = 0.18 zhist = RedshiftHistogram(cat_rand, FSKY, cosmo_nbody, redshift='Z') alpha = 1.0 cat_data.csize / cat_rand.csize nofz = InterpolatedUnivariateSpline(zhist.bin_centers, alphazhist.nbar) cat_data['NZ'] = nofz(cat_data['REDSHIFT_ESTIMATE']) cat_rand['NZ'] = nofz(cat_rand['Z'])

fkp = FKPCatalog(cat_data, cat_rand, BoxPad=0.02)

fkp['data/NZ'] = nofz(cat_data['REDSHIFT_ESTIMATE']) fkp['randoms/NZ'] = nofz(cat_rand['Z'])

fkp['data/FKPWeight'] = 1.0 / (1.0 + fkp['data/NZ'] pk_fid) #pk_fid set to 10000 fkp['randoms/FKPWeight'] = 1.0 / (1.0 + fkp['randoms/NZ'] pk_fid)

mesh = fkp.to_mesh(Nmesh=Nmesh, window=window, compensated=True, interlaced=True, nbar='NZ', fkp_weight='FKPWeight')

rainwoodman commented 3 years ago

Did this affect the estimation of the power spectrum?

Keep in mind that the box is periodic -- it appears (at a glance) that if you shift the periodic box by half a box in x and y the two results (nbodykit and yours) match.

Could you mind also sharing a subsample of the randoms (e.g. a uniform subsample of 1000 points spanning the area)?

danpryer commented 3 years ago

Ok.. excuse my stupidity.. its because I had window set to 'TSC' instead of 'tsc' ... apologies ☠️

rainwoodman commented 3 years ago

Glad you sorted this out.

I suppose the reason for the 'weird' geometry in nbodykit's internal mesh is that some part of the FKP assumes the observer sits at 0, 0, 0. If the center is moved to center of box then the line of sight direction for computing the poles are incorrect.

If that's the case, maybe adding an observer=(ox, oy, oz) argument to FKP catalog?

danpryer commented 2 years ago

Thanks! And I guess it's not a big issue really, at the end of the day the results it produces are correct :)