Happy Halloween Everyone!

## Sunday, October 31, 2010

## Monday, October 25, 2010

### 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

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)

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 =

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:

So their convention is the opposite of mine above so we should use the following:

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

theta =

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:

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

Labels:
code,
coordinate,
mask,
randoms,
spherical

## Friday, October 22, 2010

### 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.

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.

## Wednesday, October 20, 2010

### 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:

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

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:

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

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

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

stripe82targets = stripe82targets[wphot]

stripe82quasars = stripe82quasars[wspec]

stripe82data = stripe82data[wdata]

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

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

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> 20.03

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_

(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, d12totaltarget = n_elements(where(

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> stripe82quasars.dec, 2./3600, i1, i2, d12

> print, n_elements(where(

> 0.543214))/45.93

> 10.47

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, d12totaltarget = n_elements(where(

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_

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,

; Low Redshift Quasars targeted per square degree

lqpsd = n_elements(where(

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-

;targets = mrdfits(erinfile, 2)

adamfile = '/home/jessica/boss/chunk1_

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_

; Read in the truthtable file

truthfile = '/home/jessica/boss/BOSS_

collateinfo = mrdfits(truthfile, 1)

;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.15confcut = where(stripe82data.z_conf_

stripe82data.Z_PERSON GE 2.15)

lowqsocut = where(stripe82data.z_conf_

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)).areawspec = where(stripe82quasars.

wdata = where(stripe82data.

;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/

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 thresholdtpsd = round(20.0*area)

sortRatio = targets[reverse(sort(targets.

threshold = sortRatio[tpsd-1]

print, compthresh

;Area

print, thisArea

;Likelihood Targets per square degree

thistpsd = n_elements(where(

threshold))/thisArea

print, thistpsd

; Total objects targeted

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

stripe82data.dec, 2./3600, i1, i2, d12totaltarget = n_elements(where(

threshold))/thisArea

print, totaltarget

spherematch, stripe82targets.ra, stripe82targets.dec,

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

qpsd = n_elements(where(

threshold))/thisArea

print, qpsd

spherematch, stripe82targets.ra, stripe82targets.dec,

; Likelihood Low Redshift Quasars targeted per square degree

lqpsd = n_elements(where(

threshold))/thisArea

print, lqpsd

;Likelihood targeted not qsos

print, totaltarget - qpsd - lqpsd

;Total z > 2.15 qsos

totalqsos = n_elements(stripe82quasars)/

print, totalqsos

;Total 0.5 < z < 2.15 qsos

totallowqsos = n_elements(stripe82lowquasars)

print, totallowqsos

Labels:
Adam Myers,
completeness,
efficiency,
likelihood

## Tuesday, October 19, 2010

### 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:

## Monday, October 18, 2010

### 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:

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:

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

Labels:
Alexia Schulz,
d,
David Schlegel,
Shirley Ho

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

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:

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:

- How do I assign a random place in a HEALPix pixel and convert this to a ra/dec position?
- 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?
- 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)?

Labels:
coordinate,
geometry,
HEALPix,
mask,
randoms,
Shirley Ho

Subscribe to:
Posts (Atom)