Auger points out in their paper that the binomial probability formula gives only the probability P of coming up with a greater alignment for the given set of parameters (redshift, angular window, energy cutoff). It does not account for the fact that these parameters were not chosen at random - they were chosen to maximize the correlation. In order to take this into account, they generate randomized data sets and repeat the parameter scans on each simulated data set. The fraction of data sets with a smaller P value gives the true probability of the data being drawn from an isotropic distribution.
I decided to spend some time optimizing my code so that I could run it on a large number of data sets in a reasonable amount of time. The new code is here: agn macro. It can now search 1 data set per minute (on deneb).
The random data sets should be isotropic, ie, drawn from a distribution proportional to the sky exposure. The energy distribution of each simulated data set is the same as for the real data set - in other words, for each shower in the real data set, I assigned a new galactic latitude/longitude from a shower drawn at random from a large set of simulated showers. (The pattern of sky exposure of CHICOS does not vary significantly with the shower energy, so I did not worry about matching the shower energy with the energy of the simulated shower.) Simulated showers had to pass a cut on angle, theta < 1 radian.
First, I have checked that this process gave a good distribution of simulated showers. Below is the sky distribution of a few simulated data sets, and a histogram of the angular separation between the simulated showers and the set of AGN. It is what you would expect from a random distribution of points on a sphere, which is good.
The first correlation I did, with z=0.018 and angle/energy allowed to vary, resulted in a best value of P=0.032. However, when I did the same analysis on a set of 1000 simulated data sets, I found that nearly 50% had a smaller value of P. So this result is actually average for an isotropic distribution.
When I expanded the scan on the data to include z, I found a new smallest value of P=.02. However, when this expanded search is done on randomized data, ~85% of the data sets come up with a better correlation for some combination of z, angle, energy.
I am forced to conclude that we do not see a real correlation after all.
For comparison, Auger's best P value was 5x10^-9, with an overall probabilty (from the random data sets) of <10^-5.