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

#134 Fix error in ParallelTransport when multiple objects are all at the center. #135

Closed rmjarvis closed 3 years ago

rmjarvis commented 3 years ago

This fixes the bug reported by @SandraUnruh in issue #134.

Sandra noticed that catalogs with duplicated objects could result in xi+ being biased low. It only shows up in spherical coordinates and is more pronounced when there are non-unit weights.

The root problem was in the step of parallel transporting shears to the center of a cell to take the average. When all the shears were at the same place, some calculations could sometimes be numerically very tiny but not zero, when perfect math would have them be 0.0. This then shows up in a denominator. The numerator is also essentially zero, so you get 0/0, which meant the phase factor exp(i phi) was some random value, not equal to 1, which is what you should have when you don't actually need to parallel transport at all (since all the points are at the center of the cell).

I had a check for exactly 0.0 to skip the phase factor, which worked when there was only one object in a cell. But in your case with the duplicated objects, rounding errors prevented it from being exact. Essentially sin(phi) was calculated as something of order 3.e-19/1.e-17 (it varies for each cell -- probably some were worse than this), which was not exactly 0, and so the shears were rotated a little bit. The rotation was random, which had the systematic effect of diluting the signal, thus explaining the decrease you saw for the duplicated catalog.

The solution was to change the check for normA + normB from == 0 to < 1.e-12. This sum is the denominator in question and is nominally (sin(theta) cos(dec))^2 where theta is the angular distance of the point to the cell center and dec is the declination of the cell center. If this is < 1.e-12, then the object is less than an arcsec from the cell center, so we can skip the parallel transport rotation.