Bug in two_point_clustering.tpcf_multipole

[johannesulf]

@duncandc As mentioned in #649, based on some science tests I believe that the multipoles calculated by halotools are a factor of 2 too low. Mathematically, the multipoles are the coefficients for expanding xi(s, mu) into a series of Legendre polynomials P_l (mu), where l is the order. Take the following code.

import numpy as np
from halotools.mock_observables import tpcf_multipole

xi_mu_s_fake = np.ones((49, 49))
mu_bins = np.linspace(0, 1, 50)
print tpcf_multipole(xi_mu_s_fake, mu_bins, order=0)

We calculate the monopole (order = 0) of a 2PCF that is 1 at all radii and LOS angles. Clearly, the monopole should be 1 at all radii because the 0th order Legendre polynomial is also 1. However, I get an array of 0.5. A factor of 2 too low. Let's look at the source code of tpcf_multipole.

# process inputs
s_mu_tcpf_result = np.atleast_1d(s_mu_tcpf_result)
mu_bins = np.atleast_1d(mu_bins)
order = int(order)

# calculate the center of each mu bin
mu_bin_centers = (mu_bins[:-1]+mu_bins[1:])/(2.0)

# get the Legendre polynomial of the desired order.
Ln = legendre(order)

 # numerically integrate over mu
 result = (2.0*order + 1.0)/2.0 *\
        np.sum(s_mu_tcpf_result * np.diff(mu_bins) * Ln(mu_bin_centers), axis=1)

This code calculates the coefficients/multipoles by integrating over mu. The problem is that mu is only defined in [0, 1] in halotools. The source code of tpcf_multipole explicitly states that and s_mu_tpcf also only accepts [0, 1]. However, the integration should be going from -1 to 1 in mu, not 0 to 1. That means all multipoles of even order (monopole, quadropole, hexadecapole...) where the Legendre polynomial is even are currently underestimated by a factor of 2. On the other hand, all multipoles of odd order now give large non-zero values. However, they should be zero because their Legendre polynomials are odd and the 2PCF even as a function of mu. Does that make sense?

An easy fix would be to change the last line to

 # numerically integrate over mu
 result = (2.0*order + 1.0)/2.0 *\
        np.sum(s_mu_tcpf_result * np.diff(mu_bins) * (Ln(mu_bin_centers) + Ln(-mu_bin_centers)), axis=1)

[duncandc]

I think @johannesulf's fix here is exactly the right thing to do. I can implement this, or Johannes can since its his good eye...

@aphearin

[duncandc]

This is also a good test. The odd multipoles should be zero to within numerical precision.

[aphearin]

After some discussion with @duncandc, we decided to include a brute force comparison of npairs_s_mu to a pure python implementation in order to consider this resolved. See #711.

At some point, it would be useful to find someone who has an independently written multipole calculation that uses the distant observer approximation, but since we do not know of such a code, we'll stick with improving npairs_s_mu testing for the time being.