Happy Halloween!

Happy Halloween Everyone!

Jessica as 'Hit Girl'

Code for Random Catalog with Mask

From Wolfram Math World:

To pick a random point on the surface of a sphere, it is incorrect to select spherical coordinates θ and φ from uniform distributions θ ∈ [0,2π) and φ ∈ [0,π], since the area element dΩ=sinφdθdφ is a function of φ, and hence points picked in this way will be "bunched" near the poles (left figure above).

To obtain points such that any small area on the sphere is expected to contain the same number of points (right figure above), choose u and v to be random variates on (0,1). Then

θ = 2πu

φ = cos^-1(2v - 1)

gives the spherical coordinates for a set of points which are uniformly distributed over the surface of the sphere. This works since the differential element of solid angle is given by

dΩ = sinφdθdφ = -dθd(cosφ)

So for us to pick a random ra and dec on the sphere we want to use the following in idl:

theta = 2*!pi*RANDOMU(S, 1)
phi = acos(2*RANDOMU(S, 1) - 1)

This should give us a random theta and phi where θ ∈ [0,2π) and φ ∈ [0,π].

Unfortunately, just to be difficult HEALPix expects the following as inputs for ang2pix_nest:

theta - angle (along meridian), in [0,Pi], theta=0 : north pole
phi - angle (along parallel), in [0,2*Pi]

So their convention is the opposite of mine above so we should use the following:
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

Here is the code:

file = "maskfile.fits"
; nside_in is the number of pixels n the mask
read_fits_map, file, weight, nside=nside_in, ordering=order

; make random phi and theta
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

; Find pixel ipnest
ang2pix_nest, nside_in, theta, phi, ipnest

; Find weight of pixel
thisweight = weight[ipnest]

; check to make sure pixel center (thistheta, thisphi) is close to phi and theta
pix2ang_nest, nside_in, ipnest, thistheta, thisphi

(Checked this for several theta/phi and it looked good)


So here is the code to make the random catalog of size randomsize:

randomsize = 100
randomtheta = findgen(randomsize)*0.0
randomphi = findgen(randomsize)*0
index = 0

while index lt randomsize do begin
    thisindex = index 
    while index eq thisindex do begin
       phi = 2*!pi*RANDOMU(S, 1)
       theta = acos(2*RANDOMU(S, 1) - 1)
       ang2pix_nest, nside_in, theta, phi, ipnest
       thisweight = weight[ipnest]
       if thisweight gt 2.0 then begin
          randomtheta[index] = theta
          randomphi[index] = phi
          index = index+1
       endif  
    endwhile
endwhile

Randoms with a Mask

I wrote code to generate the randoms using Shirley's mask. Ran it, and now need to compare it to the actual data. The actual data isn't in a fits file.. it's an ASCII file with no header, so I am not exactly sure which column is what. Shirley left the country, so she is out of communication right now. Soon I should be able to check that the randoms fall over the data. Once I do that, I can run the correlation functions and see how the reconstruction looks on this data.

Exciting!

Alexia: Tomorrow is Cosmology in Northern California, so I will be busy with that all day. But it would be good for us to touch base on Monday again.

Doing Efficiencies Correctly

Adam and I have been going back and forth about the best way to compute the efficiencies and completeness for the likelihood paper.

Below is an email I sent him tonight:

> print, n_elements(where(stripe82targets.like_ratio_core gt 0.543214))/45.93
> 20.03

I get the same number as you here: 20.0305

But I am confused as to why you are calling this number to total
fibers. Wouldn't the fibers be the match between stripe82targets gt
threshold and all those in the BOSS_Quasars_3PCplus_v1.1_

wcomp.fits?
(I call this stripe82data)

Because what you calculate above is all potential targets that pass
the threshold, but a large fraction of them were not observed when you
look at the spherematch between stripe82targets and stripe82data:

; Total objects targeted

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82data.ra,

stripe82data.dec, 2./3600, i1, i2, d12
totaltarget = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, totaltarget
11.5174

Which means that of the 20 objects per square degree we expected to
target, only 12.7 were actually targeted... this makes since when we
are using a spec completeness of 0.5... and should go up as we
increase go down in the table correct?

I agree that the area = 45.93 and the threshold = 0.543214

However when I do the below, I get:

> spherematch, stripe82targets.ra, stripe82targets.dec, stripe82quasars.ra,
> stripe82quasars.dec, 2./3600, i1, i2, d12
> print, n_elements(where(stripe82targets[i1].like_ratio_core gt
> 0.543214))/45.93
> 10.47

10.2330

Not sure where the 0.74 difference is coming from....

When I look at the total targets number of objects per square degree
that passed the likelihood threshold but were targeted, this was:

; Total targeted

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82data.ra,

stripe82data.dec, 2./3600, i1, i2, d12
totaltarget = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, totaltarget
11.5174

Which means that we have only 1.284 targets per square degree that
were targeted by likelihood and not quasars correct?

The number that are lower redshift quasars is:

lowqsocut = where(stripe82data.z_conf_person GE 2 and
stripe82data.Z_PERSON GE 0.5 and stripe82data.Z_PERSON LT 2.15)
stripe82lowquasars = stripe82data[lowqsocut] ; with 0.5 < redshift < 2.15

spherematch, stripe82targets.ra, stripe82targets.dec,

stripe82lowquasars.ra, stripe82lowquasars.dec, 2./3600, i1, i2, d12

; Low Redshift Quasars targeted per square degree
lqpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, lqpsd
0.326581



In terms of filling out your table I get the following numbers:
Spec Area Total # Like z > 2.15 other not ALL > 2.15 ALL other
Comp fibers targets QSOs QSOs QSOs QSOs QSOs

0.5 45.93 20.03 11.52 10.23 0.33 0.96 23.32 6.79
0.6 37.80 20.45 12.43 11.03 0.37 1.0 24.71 8.25
0.7 26.12 21.11 12.74 10.78 0.53 1.42 24.99 11.73
0.8 14.02 22.40 13.27 10.42 0.85 1.99 23.19 22.26
0.9 2.12 18.35 10.82 8.47 0.94 1.41 20.24 146.86


Below is the code I used to generate the table.
Let me know if you have thoughts.
Jessica

; Read in the targeting file
;erinfile = '/home/jessica/boss/boss-qso-stripe82median.fits'
;targets = mrdfits(erinfile, 2)
adamfile = '/home/jessica/boss/chunk1_wcomp.fits'
targets = mrdfits(adamfile, 1)

;stripe82cut = where((targets.dec GE -1.25) and (targets.dec LE 1.25)
and ((targets.ra GT 317) or (targets.ra LT 45)))
;stripe82targets = targets[stripe82cut]
stripe82targets = targets

newratio = stripe82targets.LIKE_RATIO_CORE

; Read in the truthtable file
truthfile = '/home/jessica/boss/BOSS_Quasars_3PCplus_v1.1_wcomp.fits'
collateinfo = mrdfits(truthfile, 1)

stripe82data = collateinfo

;High human confidence quasars (with z > 2.15)
confcut = where(stripe82data.z_conf_person GE 2 and
stripe82data.Z_PERSON GE 2.15)

stripe82quasars = stripe82data[confcut] ; with redshift > 2.15
lowqsocut = where(stripe82data.z_conf_person GE 2 and
stripe82data.Z_PERSON GE 0.5 and stripe82data.Z_PERSON LT 2.15)
stripe82lowquasars = stripe82data[lowqsocut] ; with 0.5 < redshift < 2.15


compthresh = 0.9

wphot = where(stripe82targets.completeness ge compthresh)
wspec = where(stripe82quasars.completeness ge compthresh)
wdata = where(stripe82data.completeness ge compthresh)

;plot, stripe82targets.ra, stripe82targets.dec, psym=3
;oplot, stripe82targets[wphot].ra, stripe82targets[wphot].dec, psym=3,
color=FSC_COLOR('blue')
;oplot, stripe82quasars[wspec].ra, stripe82quasars[wspec].dec, psym=7,
color=FSC_COLOR('red')

areas = mrdfits('/home/jessica/boss/area_wcomp.fits',1)

thisArea = areas(where(areas.comp_thresh eq compthresh)).area

stripe82targets = stripe82targets[wphot]
stripe82quasars = stripe82quasars[wspec]
stripe82data = stripe82data[wdata]

area = 219.9 ; we know the exact area of the chunk1 file
tpsd = round(20.0*area)
sortRatio = targets[reverse(sort(targets.LIKE_RATIO_CORE))].LIKE_RATIO_CORE
threshold = sortRatio[tpsd-1]

;Comp threshold
print, compthresh

;Area
print, thisArea


;Likelihood Targets per square degree
thistpsd = n_elements(where(stripe82targets.like_ratio_core gt
threshold))/thisArea
print, thistpsd

; Total objects targeted

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82data.ra,

stripe82data.dec, 2./3600, i1, i2, d12
totaltarget = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, totaltarget

spherematch, stripe82targets.ra, stripe82targets.dec,
stripe82quasars.ra, stripe82quasars.dec, 2./3600, i1, i2, d12

; Likelihood Quasars targeted per square degree
qpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, qpsd

spherematch, stripe82targets.ra, stripe82targets.dec,

stripe82lowquasars.ra, stripe82lowquasars.dec, 2./3600, i1, i2, d12

; Likelihood Low Redshift Quasars targeted per square degree
lqpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt
threshold))/thisArea
print, lqpsd


;Likelihood targeted not qsos
print, totaltarget - qpsd - lqpsd


;Total z > 2.15 qsos
totalqsos = n_elements(stripe82quasars)/thisArea
print, totalqsos

;Total 0.5 < z < 2.15 qsos
totallowqsos = n_elements(stripe82lowquasars)/thisArea
print, totallowqsos

Cosmology in Northern California

In the small chance that anyone (other than Alexia) reads this blog, I feel like I should promote the conference I am organizing considering it takes away so much of my research time these days. Check it out:

http://panic.berkeley.edu/~jessica/cinc/

Next steps...

I talked to Shirley and Alexia and they both agree that the procedure outlined in my last post makes sense. Shirley told me that the weights in her file don't mean anything, so I should just treat them is a step function. I need to talk to David to make sure he agrees with all this.

Thinking I should post-pone my trip to Los Alamos until December. Feeling like I'd like to have the data in better shape before I get there, and the flights are cheaper.

Steps:

1. Talk to David and make sure he agrees with these steps. Also talk about funding for trip.
2. Generate randoms as laid out in last post.
3. Book tickets for December.
4. Keep looking for temp housing in Los Alamos.

Another idea...

Another idea I have about generating this random catalog. First I generate the random ra and random dec, then I use the function ang2pix_nest (converts angle to pixel) to figure out which pixel it occurs in. Once we have that, depending on the weight we can keep the random value or throw it away. However, this might take a really long time if most of our objects are outside of the mask.

I also am not sure how to use the weight value. Should it be a step function, if weight > 0.2 then keep, if not then throw away, or something more subtle like the weight represents the expectation value of the probability distribution of if the random is in the pixel or not. In terms of the redshifts, should I do what I did before and match the random redshift distribution to the data?

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)?