Showing posts with label HEALPix. Show all posts
Showing posts with label HEALPix. Show all posts

Tuesday, November 2, 2010

Masking both Data and Randoms

I've masked both the data and the randoms now:


The randoms (magenta) now fall exactly on top of the data (cyan). The code to do this is as follows:

;Read in data
file1 = "/home/shirleyho/research/SDSS_PK/powerspec_dir/DATA/ra_dec_z_run_rerun_camcol.z_gt_0.01"

readcol, file1, RA, DEC, z, RUN, rerun, camcol, x, x, x, x, x, x, x, x, format='(F,F,F,D,D,D,D,D,F,F,F,F,F,F)'

datasize = n_elements(ra)

;Read in mask
file = "/home/shirleyho/research/SDSS_PK/powerspec_dir/MASK_SDSS/dir_fits/tmp_unified_score_ebv_psf_bsm_nside1024.fits"

read_fits_map, file, weight, nside=nside_in, ordering=order


index = 0L
;Convert ra/dec to phi/theta expected in HEALPix
dataphi = ra*!pi/180.
datatheta = (90.0-dec)*!pi/180.
;Will be filled in with weight corresponding to ra/dec of data
dataweight = findgen(datasize)*0.0

;Function I wrote to get the weight for the data
junk = maskdata(datatheta, dataphi, index, datasize, nside_in, weight, dataweight)
;Only keeps data with a weight greater than 0.2 (inside mask)
datacut = where(dataweight GE 0.2)

;Make the randoms
randomsize = datasize*5L
randomtheta = findgen(randomsize)*0.0
randomphi = findgen(randomsize)*0.0
index = 0L
;Function I wrote to fill in randomtheta and randomphi with points inside the mask
junk = makemaskrandoms(randomtheta, randomphi, index, randomsize, nside_in, weight)

;Plot the data and the randoms
plot, ra[sub], dec[sub], psym=3, symsize=2, xr=[0,360],yr=[0,90], XTITLE = "RA", YTITLE = "DEC", TITLE = "RA vs Dec of Data and Randoms", charsize = 2, charthick = 1, thick = 2, xthick=2, ythick=2
oplot, ra[datacut], dec[datacut], psym=3, color=fsc_color('cyan')
oplot, randomra[0:100000], randomdec[0:100000], psy = 3, color=fsc_color('magenta')

Now to compute some correlation functions!

Monday, October 18, 2010

Learning About HEALPix

I've been playing with Shirley's LRG catalog and mask. The mask uses something called HEALPix which I've never used before. It basically breaks up the spherical sky into equal area pixels:



The map is in the form of weights for each pixel. Shirley says that if the weight is greater than 0.2 she considers that pixel in the mask.

So to create a random catalog, I just go through the weights, and if it is greater than 0.2 I place objects in that pixel.

The key is to now convert from a Healpix mask to ra/dec.
There are two routines that do this that I have found:

DPF_PIX2ANG_NEST
pix2ang_nest

The problem is that these routines return the ra/dec of the center of the pixel. And to avoid putting extra power at the separations of the pixels, we want to select random points within the pixel. However, as you can see below, the ra range of these pixels changes dramatically as a function of cos (theta) (i.e. large at the poles and small at the equator) and so it isn't a simple procedure to fill in the area of a pixel evenly with random points:


I've searched to see if there is an existing package that does this, and I can't seem to find anything. In the HEALPix documentation there is a comment on random number generation, it doesn't seem to be what I am interested in here. I also have a email out to Shirley about this.

There are many functions that deal with HEALPix on their web page and well as in IDLUTILS. I'm currently combing through them to see if there is a function that gives me the ra/dec ranges for a given pixel, or even better, automatically picks a random point given a pixel.

Some questions:
  1. How do I assign a random place in a HEALPix pixel and convert this to a ra/dec position?
  2. The mask contains a "weight" for each HEALPix pixel, how to do factor in the weight into the random catalog. My intuition is that pixel's with higher weight should have more objects in them, but will this add clustering signal to the random catalog that we don't want?
  3. How do we distinguish what is signal because of observation effects (we only looked at that patch of the sky once, or we didn't target as many objects there) and what is signal due to clustering of the data (there is a huge galaxy cluster in that part of the sky, and that is why there is more data there)?