search

rmjarvis / TreeCorr

Code for efficiently computing 2-point and 3-point correlation functions. For documentation, go to
http://rmjarvis.github.io/TreeCorr/
Other
98 stars 37 forks source link

feature request: halo ellipticity estimators #36

Closed clampitt closed 8 years ago

clampitt commented 8 years ago

Similar to running an NG (tangential shear) correlation, but the user passes in an extra parameter, call it phi, for each object in the lens catalog. Treecorr returns the result of two optimally weighted estimators of the elliptical lensing shear, and two corresponding cross-components. All four results are arrays of length N, where N is the number of bins in angle (or perhaps projected scale in Mpc/h).

This extra input parameter phi is used to rotate the shear field for each lens into a common Cartesian coordinate system (gamma_1 and gamma2 instead of gamma+ and gamma_x). Having performed this rotation, the desired estimators to calculate are gamma_4theta = Sum_i (w_i * gamma_1i * cos(4theta) + w_i * gamma_2i * sin(4theta)) / Sum_i (w_i) gamma_const = Sum_i (w_i * gamma_1i) / Sum_i (w_i)

gamma_4theta_x = Sum_i (w_i * gamma_1i * cos(4theta) - w_i * gamma_2i * sin(4theta)) / Sum_i (w_i) gamma_const_x = Sum_i (w_i * gamma_2i) / Sum_i (w_i)

where the index i runs over all lens-source pairs, and w_i is the usual weight for gamma (for example, sum of squares of measurement error and shape noise).

rmjarvis commented 8 years ago

On the plane back from Madrid, I realized that these estimators are actually just the normal xi+ and xi- for the cross correlation between the halo ellipticities and the sources, where the halo ellipticity is normalized to have |g| = 1 (i.e., it is just the phase g1 + i g2 = exp(i phi)).

I made a unit test to confirm this, which is now in the code here (on master, but it works with the current v3.2.1 release -- no extra source code required).

The basic code flow is along the lines of the following:

lens_absg = numpy.sqrt(lens_g1**2 + lens_g2**2)
lens_cat = treecorr.Catalog(ra=lens_ra, dec=lens_dec, g1=lens_g1/lens_absg, g2=lens_g2/lens_absg)
source_cat = treecorr.Catalog(ra=source_ra, dec=source_dec, g1=source_g1, g2=source_g2, w=source_w)
gg = treecorr.GGCorrelation(min_sep=1, max_sep=30, bin_size=0.1)
gg.process(lens_cat, source_cat)
gamma_4theta = gg.xim
gamma_const = gg.xip
gamma_4theta_x = gg.xim_im
gamma_const_x = gg.xip_im

The test I linked to doesn't actually use weights, but those are already implemented. And it uses x,y, just because it is easier to construct the test in Euclidean coordinates. But ra, dec would normally be used on real data.

rmjarvis commented 8 years ago

So, given the above, do you think there should be a helper function in TreeCorr for this? It's pretty straightforward, despite being non-obvious, so I'm not sure how helpful a helper function would actually be.

Perhaps just some additional documentation to point this out? Maybe in the GGCorrelation class doc string.

Or maybe just in your paper about this, you lay out the math to prove this point, in which case it would be clear to other people making this measurement what functions in TreeCorr to use.

rmjarvis commented 8 years ago

I think there is nothing further to do on this issue. Closing now.