Event list read in mode has a small bug

[clark2668]

The event list read-in mode has a small "bug." AraSim displaces the vertices of all events that are read-in by 10km. So it does x_new = x_old + displacement, and y_new = y_old + displacement. You have to make sure to correct for this when you do the event list generation. If you don't, you end up with about a 1% under-estimation of the fiducial volume. Here's a little toy MC (written by Claude) which does the calculation and comparison.

So you either need to account for the shift at generation, or AraSim should really rotate the Earth, rather than do a literal lateral translation without also fixing the depth points.

#!/usr/bin/env python3
"""
Monte Carlo calculation of the volume of a cylinder that lies inside the Earth,
where the top of the cylinder is tangent to the Earth's surface.
Batched to avoid memory issues.
"""

import numpy as np

# Constants
R_EARTH = 6371.0  # km

# Cylinder parameters
r_cyl = 15.0  # km
h_cyl = 3.0   # km

# Monte Carlo - batched
n_total = 100_000_000
batch_size = 10_000_000
n_batches = n_total // batch_size

rng = np.random.default_rng(seed=42)

n_inside_total = 0
n_still_inside_total = 0

for i in range(n_batches):
    # Uniform sampling in cylinder
    r = np.sqrt(rng.uniform(0, r_cyl**2, batch_size))
    theta = rng.uniform(0, 2*np.pi, batch_size)
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    z = rng.uniform(R_EARTH - h_cyl, R_EARTH, batch_size)

    # Check which points are inside Earth
    inside_earth = (x**2 + y**2 + z**2) <= R_EARTH**2
    n_inside_total += np.sum(inside_earth)

    # Student mistake: shift x,y for points inside Earth
    x_good = x[inside_earth]
    y_good = y[inside_earth]
    z_good = z[inside_earth]

    shift = np.sqrt(2) * 10.0  # km
    x_shifted = x_good + shift
    y_shifted = y_good + shift

    still_inside = (x_shifted**2 + y_shifted**2 + z_good**2) <= R_EARTH**2
    n_still_inside_total += np.sum(still_inside)

# Results
V_cylinder = np.pi * r_cyl**2 * h_cyl
fraction_inside = n_inside_total / n_total
V_inside = V_cylinder * fraction_inside

fraction_still = n_still_inside_total / n_inside_total
V_still = V_inside * fraction_still

print(f"Cylinder volume:        {V_cylinder:.2f} km³")
print(f"Volume inside Earth:    {V_inside:.2f} km³")
print(f"Fraction inside:        {fraction_inside:.6f}")

print(f"\n--- After student mistake (shift x,y by sqrt(2)*10 km) ---")
print(f"Fraction still inside:  {fraction_still:.6f}")
print(f"Volume still inside:    {V_still:.2f} km³")