Fix theta uncertainty bug

[napoliion]

Before you submit this PR: make sure to put all operations-related information in a wiki-note, a PR should be about code and is publicly accessible

What does the code in this PR do / what does it improve?

This fixes a bug found by @sebvetter (figure below) in how theta uncertainties were handled due to contours that passed over the -pi and pi boundary. These contours would have an uncertainty that is an order of magnitude smaller than the contours that do not cross this boundary.

Can you briefly describe how it works?

The theta uncertainty is calculated as the difference between the min and max $\theta$ of a contour. The correction is to rotate the contour around (0,0) by the average $\theta$ within the contour so that the contour is now around $\theta = 0$ to avoid the seam, then add $\pi$ and modulo $2\pi$ to handle large angular contours.

# Correction for circular nature of angle 
avg_theta = np.arctan2(
    np.mean(events[f"{type_}_position_contour_cnf_naive"][..., 1], axis=1),
    np.mean(events[f"{type_}_position_contour_cnf_naive"][..., 0], axis=1),
)

avg_theta = np.reshape(avg_theta, (avg_theta.shape[0],1))
theta_array_shift = (np.subtract(theta_array, avg_theta) + np.pi) % (2 * np.pi)
theta_min = np.min(theta_array_shift, axis=1)
theta_max = np.max(theta_array_shift, axis=1)

theta_diff = theta_max - theta_min

Can you give a minimal working example (or illustrate with a figure)?

Given a generated fake contour, we can show that this method works for the three cases:

• edge case described above at the $-\pi$ and $+\pi$ boundary
• contour that overlaps (0,0)
• small angle contour with $\theta < \pi$

Below is a quick minimal modified example of the code above on a fake generated contour located at the seam we want to test this method on:

# Generate fake contour located on -pi and pi seam
angles = np.linspace(0, 2 * np.pi, 16, endpoint=False)
radius_uncertainty = 2 # width of contour
xdist_from_origin = 30 # x dist of contour from (0,0)

# Calculate x and y coordinates for each angle
circle_contour = np.array([
    [radius_uncertainty * np.cos(angle), radius_uncertainty * np.sin(angle)] for angle in angles])

circle_contour = np.vstack([circle_contour, circle_contour[0]])

# Calculate expected theta diff at theta = pi/2 to calculate
check_contour = circle_contour.copy()
circle_contour[:,1] += np.abs(xdist_from_origin) # shift up y-axis
test_theta_array = np.arctan2(
    circle_contour[..., 1],
    circle_contour[..., 0],
)
theta_min = np.min(test_theta_array)
theta_max = np.max(test_theta_array)

theta_diff = theta_max - theta_min

expected_theta_diff = theta_diff/2

print(f'Theta Diff: {expected_theta_diff:.12f}')

which prints:

Theta Diff: 0.066568163776

After generating the above contour, the contour can be moved to show that this still returns the expected values at angles on and off of the $-\pi$ and $\pi$ seam.

# Now shift to -pi and pi seam for testing correction
# Define cases
test_case1 = check_contour.copy() # over -pi and pi seam
test_case2 = check_contour.copy() # at theta = pi/2

test_case1[:,0] -= xdist_from_origin # shift left on x-axis
test_case2[:,1] += np.abs(xdist_from_origin) # shift up on y-axis

# Calculate if theta_diff matches expected theta_diff
for contour in [test_case1, test_case2]:
    test_theta_array = np.arctan2(
        contour[..., 1],
        contour[..., 0],
    )
    
    avg_theta = np.arctan2(
        np.mean(contour[..., 1]),
        np.mean(contour[..., 0]),
    )

    # Correction for circular nature of angle
    theta_array_shift = (np.subtract(test_theta_array, avg_theta) + np.pi) % (2 * np.pi)
    theta_min = np.min(theta_array_shift)
    theta_max = np.max(theta_array_shift)
    
    theta_diff = theta_max - theta_min
    print(f'Theta Diff: {expected_theta_diff:.12f}')

which also prints:

Theta Diff: 0.066568163776
Theta Diff: 0.066568163776

Because the contour is not a perfect circle in the above example (only 17 points are selected, as per the CNF contour array shape), rotating the contour will cause slight differences in theta_diff depending on the orientation of this contour.

To check that this works on data, I have loaded a set of events from a Kr83m run to test. The above figure was reproduced and after rerunning with the code change, the bug is no longer visible.

For event contours:

For peak contours:
(The sharp dip in the center is a result of very small radii close to 0.)

[coveralls]

coverage: 89.797% (+0.03%) from 89.771%
when pulling e4c73c5 on theta_uncertainty_fix
into 6838d38 on master.

[dachengx]

Thanks @napoliion .

In general, I think we should consider a corner case where theta_array_shift still crosses the $-\pi$ / $\pi$ boundary.

[juehang]

In addition to the specific comments, the shared code should be refractored into a static method as @dachengx suggested, and the inconsistency in naming should be fixed.

[dachengx]

In addition to the specific comments, the shared code should be refractored into a static method as @dachengx suggested, and the inconsistency in naming should be fixed.

@juehang is referring to position_contours_cnf -> position_contour_cnf (slack thread)