1st Attempt at QSO-Galaxy Cross Correlation

I've calculated the qso-galaxy correlation functions on the entire stripe, and also the stripe divided into two halves.
The code to do this is below (also in ../qsobias/Correlate/runCorrelation.py):

import numpy as N
from pylab import *
from correlationFunctions import *

#------------------------------------------------------------------------
# Set Angular Correlation Bin Information
#------------------------------------------------------------------------

oversample = 5. # Amount that randoms should be oversampled
corrBins = 15.0 # Number of correlation bins (+1)
mincorr = 0.1 # (Mpc/h comoving distance separation) Must be great than zero if log-binning
maxcorr = 10.0 # (Mphc/h comoving distance separation)
convo = 180./pi # conversion from degrees to radians
tlogbin = 1 # = 0 for uniform spacing, = 1 for log spacing in theta

#------------------------------------------------------------------------
# Create file names (full catalog)
#------------------------------------------------------------------------
workingDir = 'fullrun1'
makeworkingdir(workingDir)
galaxyDataFile, qsoDataFile, randomDataFile, corr2dCodefile, argumentFile, \
runConstantsFile = makeFileNamesFull(workingDir)

#------------------------------------------------------------------------
# Write run constants to a file (full catalog)
#------------------------------------------------------------------------

writeRunConstantsToFile(runConstantsFile, galaxyDataFile, qsoDataFile, \
randomDataFile, corr2dCodefile, argumentFile, oversample, corrBins, \
mincorr, maxcorr, tlogbin)

#------------------------------------------------------------------------
# Compute the Angular Correlation Function (full catalog)
#------------------------------------------------------------------------

runcrossCorrelation(workingDir, argumentFile, corr2dCodefile, galaxyDataFile,\
qsoDataFile, randomDataFile, mincorr, maxcorr, corrBins, tlogbin)

#------------------------------------------------------------------------
# Create file names (Left catalog)
#------------------------------------------------------------------------
workingDir = 'leftrun1'
makeworkingdir(workingDir)
galaxyDataFile, qsoDataFile, randomDataFile, corr2dCodefile, argumentFile, \
runConstantsFile = makeFileNamesLeft(workingDir)

#------------------------------------------------------------------------
# Write run constants to a file (Left catalog)
#------------------------------------------------------------------------

writeRunConstantsToFile(runConstantsFile, galaxyDataFile, qsoDataFile, \
randomDataFile, corr2dCodefile, argumentFile, oversample, corrBins, \
mincorr, maxcorr, tlogbin)

#------------------------------------------------------------------------
# Compute the Angular Correlation Function (Left catalog)
#------------------------------------------------------------------------

runcrossCorrelation(workingDir, argumentFile, corr2dCodefile, galaxyDataFile,\
qsoDataFile, randomDataFile, mincorr, maxcorr, corrBins, tlogbin)

#------------------------------------------------------------------------
# Create file names (Right catalog)
#------------------------------------------------------------------------
workingDir = 'rightrun1'
makeworkingdir(workingDir)
galaxyDataFile, qsoDataFile, randomDataFile, corr2dCodefile, argumentFile, \
runConstantsFile = makeFileNamesRight(workingDir)

#------------------------------------------------------------------------
# Write run constants to a file (Right catalog)
#------------------------------------------------------------------------

writeRunConstantsToFile(runConstantsFile, galaxyDataFile, qsoDataFile, \
randomDataFile, corr2dCodefile, argumentFile, oversample, corrBins, \
mincorr, maxcorr, tlogbin)

#------------------------------------------------------------------------
# Compute the Angular Correlation Function (Right catalog)
#------------------------------------------------------------------------

runcrossCorrelation(workingDir, argumentFile, corr2dCodefile, galaxyDataFile,\
qsoDataFile, randomDataFile, mincorr, maxcorr, corrBins, tlogbin)

#------------------------------------------------------------------------
# Must wait for jobs to finish running & output files to be created
#------------------------------------------------------------------------

#------------------------------------------------------------------------
# Plot Correlation Functions
#------------------------------------------------------------------------

#------------------------------------------------------------------------
# Plot Full Stripe 82
#------------------------------------------------------------------------

workingDir = "fullrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta,omega=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta, omega,'r-', label = "Full Stripe")
# Number of dd pairs in each bin
ddpairs = [3029, 5404, 9578, 17175, 31653, 58113, 106518, 194566, 358573, 656623, 1211740, 2234135, 4121513, 7602541, 13931607]
errors = 1/N.sqrt(ddpairs)
errorbar(theta, omega, yerr=errors,fmt='ro')

xlabel('Separation Distance (Mpc/h)')
ylabel('omega')
title('2D QSO-Galaxy Cross-correlation Function')
pylab.legend(loc=3)

(plots on wiki)

#------------------------------------------------------------------------
# Plot Full Stripe 82, Left-Half and Right-Half of Stripe 82
#------------------------------------------------------------------------

workingDir = "fullrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta,omega=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta, omega, label = "Full Stripe")

workingDir = "leftrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta1,omega1=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta1, omega1, label = "Left Stripe")

workingDir = "rightrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta2,omega2=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta2, omega2, label = "Right Stripe")

xlabel('Separation Distance (Mpc/h)')
ylabel('omega')
title('2D QSO-Galaxy Cross-correlation Function')
pylab.legend(loc=3)

(plots on wiki)

#------------------------------------------------------------------------
# Plot Full, Left-Half, Right-Half of Stripe 82, and Mean of all three + 3-sigma errors
#------------------------------------------------------------------------

workingDir = "fullrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta,omega=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta, omega, label = "Full Stripe")

workingDir = "leftrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta1,omega1=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta1, omega1, label = "Left Stripe")

workingDir = "rightrun1"
outfile = "./"+workingDir+"/wpsOutput.dat"
theta2,omega2=readCorrOutfile(outfile) #Reads correlation output file
loglog(theta2, omega2, label = "Right Stripe")

datamaxtrix = [omega,omega1,omega2]
mean = N.mean(datamatrix,axis=0)
errors = 3*N.std(datamatrix,axis=0)/sqrt(3)
loglog(theta, mean, 'k-',label = "Mean")
errorbar(theta, mean, yerr=errors,fmt='ko')

xlabel('Separation Distance (Mpc/h)')
ylabel('omega')
title('2D QSO-Galaxy Cross-correlation Function')
pylab.legend(loc=3)

(plots on wiki)

I sent the above plots to Martin, Alexie, Nic and David. Martin thought they looked about right. The slope of the mean correlation function above is -0.74. Martin and Nikhil found that the 3D correlation function was about -1.7, and we would expect the projected 2D correlation function to be about one less than this.... so that is good.

1) Try dividing the QSO samples into different redshift/brightness sets and see if we get different correlation results.
2) Implement jackknifing ability into the code.

Sheldon's Correlaton Function

Erin Sheldon is helping me with my 2D correlation function. He has one that runs a lot faster than the one I have been using. He has been amazing in helping me to get it working with my data set.

After some incompatibility problems with various packages we needed for running our respective code it seems to be able working now (well it compiles and passes the checks). However, when I try to run it I get the following error:

RuntimeError: found log_binsize < 0

I've pinged Erin about this, so we'll see if he has ideas. The code I'm using to test this is in the following log file: .../logs/100405log.py

7:15 Update: Erin Rocks! He just got back to me with a fix. It is running now! I'm going to do some tests to see how long it takes and compare it's answer with my correlation function's answer.

Outputs:
Generating HTM ids
Generating reverse indices
rmin: 0.001
rmax: 10
degrees?: True
nbin: 25
logrmin: -3
logrmax: 1
log binsize: 0.16
len(scale_array) = 0

Each dot is 500 points
......................................................................
35000/53071 pair count: 176696860
....................................

In [15]: rlower
Out[15]:
array([ 1.00000000e-03, 1.44543977e-03, 2.08929613e-03,
3.01995172e-03, 4.36515832e-03, 6.30957344e-03,
9.12010839e-03, 1.31825674e-02, 1.90546072e-02,
2.75422870e-02, 3.98107171e-02, 5.75439937e-02,
8.31763771e-02, 1.20226443e-01, 1.73780083e-01,
2.51188643e-01, 3.63078055e-01, 5.24807460e-01,
7.58577575e-01, 1.09647820e+00, 1.58489319e+00,
2.29086765e+00, 3.31131121e+00, 4.78630092e+00,
6.91830971e+00])

In [16]: rupper
Out[16]:
array([ 1.44543977e-03, 2.08929613e-03, 3.01995172e-03,
4.36515832e-03, 6.30957344e-03, 9.12010839e-03,
1.31825674e-02, 1.90546072e-02, 2.75422870e-02,
3.98107171e-02, 5.75439937e-02, 8.31763771e-02,
1.20226443e-01, 1.73780083e-01, 2.51188643e-01,
3.63078055e-01, 5.24807460e-01, 7.58577575e-01,
1.09647820e+00, 1.58489319e+00, 2.29086765e+00,
3.31131121e+00, 4.78630092e+00, 6.91830971e+00,
1.00000000e+01])

In [17]: counts
Out[17]:
array([ 468, 512, 758, 1216, 2176, 4266,
6962, 13332, 25388, 46232, 83912, 146620,
246088, 377354, 523490, 674092, 972164, 1683416,
3167032, 6144742, 12083396, 22744048, 41867790, 71782726,
105281816])

Likelihood Optimatization

I've been playing with various parameters in the likelihood method trying to find the most efficient cuts. The three things I have been playing with are:

added errors (what to use as the input adderr into likelihood_compute)
everything - variability (what happens when we take away variable objects from the L_everything file)
QSO + variability (what happens when we add the variable objects to the L_QSO likelihoods)

Here are my findings (these are on the co-added fluxes, next step is to re-run with single epoch):

Normal Errors (adderr = [0.014, 0.01, 0.01, 0.01, 0.014])
Percent of those targeted are quasars (based on 40/degree^2 targeting)
No variability: 0.649914
Variability everything: 0.656196
Variability everything + QSO = 0.651057

5X Errors (adderr = 5*[0.014, 0.01, 0.01, 0.01, 0.014])
Percent of those targeted are quasars (based on 40/degree^2 targeting)
No variability: 0.641348
Variability everything: 0.645346
Variability everything + QSO = 0.603655

7X Errors (adderr = 7*[0.014, 0.01, 0.01, 0.01, 0.014])
Percent of those targeted are quasars (based on 40/degree^2 targeting)
No variability: 0.627641
Variability everything: 0.624786
Variability everything + QSO = 0.572244

No Errors (adderr = 0.0*[0.014, 0.01, 0.01, 0.01, 0.014])
Percent of those targeted are quasars (based on 40/degree^2 targeting)
No variability: 0.644203
Variability everything: 0.627070
Variability everything + QSO = 0.619075

It looks like the errors we were running have the best numbers, and using the variability everything, but not the variability QSO.

I am going to play more with the definitions of variable everything and variable qso to see if I can get these to work better. I also want to play with not cutting on a L_ratio = 0.01, but perhaps changing L_ratio as a function of L_QSO (it seems we might be able to get a few more objects if we have L_ratio cut decrease as L_QSO gets large.