| Age | Commit message (Collapse) | Author |
|---|
|
Also, call this function when computing the psi parameter in nll_best().
|
|
|
|
|
|
|
|
|
|
Previously I was calculating the expected number of delta ray photons when
integrating over the shower path, but since the delta rays are produced along
the particle path and not further out like the shower photons, this wasn't
correct. The normalization of the probability distribution for the photons
produced along the path was also not handled correctly.
This commit adds a new function called integrate_path_delta_ray() to compute
the expected number of photons from delta rays hitting each PMT. Currently this
means that the likelihood function for muons will be significantly slower than
previously, but hopefully I can speed it up again in the future (for example by
skipping the shower calculation which is negligible for lower energy muons).
|
|
|
|
To enable the fitter to run outside of the src directory, I created a new
function open_file() which works exactly like fopen() except that it searches
for the file in both the current working directory and the path specified by an
environment variable.
|
|
|
|
|
|
This seems to speed things up a little bit.
|
|
Previously, the algorithm used to find peaks was to search for all peaks in the
Hough transform above some constant fraction of the highest peak. This
algorithm could have issues finding smaller peaks away from the highest peak.
The new algorithm instead finds the highest peak in the Hough transform and
then recomputes the Hough transform ignoring all PMT hits within the Cerenkov
cone of the first peak. The next peak is found from this transform and the
process is iteratively repeated until a certain number of peaks are found.
One disadvantage of this new system is that it will *always* find the same
number of peaks and this will usually be greater than the actual number of
rings in the event. This is not a problem though since when fitting the event
we loop over all possible peaks and do a quick fit to determine the starting
point and so false positives are OK because the real peaks will fit better
during this quick fit.
Another potential issue with this new method is that by rejecting all PMT hits
within the Cerenkov cone of the first peak we could miss a second peak very
close to the first peak. This is partially mitigated by the fact that when we
loop over all possible combinations of the particle ids and directions we allow
each peak to be used more than once. For example, when fitting for the
hypothesis that an event is caused by two electrons and one muon and given two
possible directions 1 and 2, we will fit for the following possible direction
combinations:
1 1 1
1 1 2
1 2 1
1 2 2
2 2 1
2 2 2
Therefore if there is a second ring close to the first it is possible to fit it
correctly since we will seed the quick fit with two particles pointing in the
same direction.
This commit also adds a few tests for new functions and changes the energy step
size during the quick fit to 10% of the starting energy value.
|
|
Also, fix a few memory leaks in test.c.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
This commit adds lots of comments to sno_charge.c and makes a couple of other
changes:
- use interp1d() instead of the GSL interpolation routines
- increase MAX_PE to 100
I increased MAX_PE because I determined that it had a rather large impact on
the likelihood function for 500 MeV electrons. This unfortunately slows down
the initialization by a lot. I think I could speed this up by convolving the
single PE charge distribution with a gaussian *before* convolving the charge
distributions to compute the charge distributions for multiple PE.
|
|
This commit adds Rayleigh scattering to the likelihood function. The Rayleigh
scattering lengths come from rsp_rayleigh.dat from SNOMAN which only includes
photons which scattered +/- 10 ns around the prompt peak. The fraction of light
which scatters is treated the same in the likelihood as reflected light, i.e.
it is uniform across all the PMTs in the detector and the time PDF is assumed
to be a constant for a fixed amount of time after the prompt peak.
|
|
|
|
This commit speeds up the likelihood function by about ~20% by using the
precomputed track positions, directions, times, etc. instead of interpolating
them on the fly.
It also switches to computing the number of points to integrate along the track
by dividing the track length by a specified distance, currently set to 1 cm.
This should hopefully speed things up for lower energies and result in more
stable fits at high energies.
|
|
|
|
To characterize the angular distribution of photons from an electromagnetic
shower I came up with the following functional form:
f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
and fit this to data simulated using RAT-PAC at several different energies. I
then fit the alpha and beta coefficients as a function of energy to the
functional form:
alpha = c0 + c1/log(c2*T0 + c3)
beta = c0 + c1/log(c2*T0 + c3).
where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3
are parameters which I fit.
The longitudinal distribution of the photons generated from an electromagnetic
shower is described by a gamma distribution:
f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a).
This parameterization comes from the PDG "Passage of particles through matter"
section 32.5. I also fit the data from my RAT-PAC simulation, but currently I
am not using it, and instead using a simpler form to calculate the coefficients
from the PDG (although I estimated the b parameter from the RAT-PAC data).
I also sped up the calculation of the solid angle by making a lookup table
since it was taking a significant fraction of the time to compute the
likelihood function.
|
|
|
|
|
|
|
|
I noticed when fitting electrons that the cquad integration routine was not
very stable, i.e. it would return different results for *very* small changes in
the fit parameters which would cause the fit to stall.
Since it's very important for the minimizer that the likelihood function not
jump around, I am switching to integrating over the path by just using a fixed
number of points and using the trapezoidal rule. This seems to be a lot more
stable, and as a bonus I was able to combine the three integrals (direct
charge, indirect charge, and time) so that we only have to do a single loop.
This should hopefully make the speed comparable since the cquad routine was
fairly effective at only using as many function evaluations as needed.
Another benefit to this approach is that if needed, it will be easier to port
to a GPU.
|
|
Since we only have the range and dE/dx tables for light water for electrons and
protons it's not correct to use the heavy water density. Also, even though we
have both tables for muons, currently we only load the heavy water table, so we
hardcode the density to that of heavy water.
In the future, it would be nice to load both tables and use the correct one
depending on if we are fitting in the heavy or light water.
|
|
path coefficients
Previously I was adding the log likelihood of the path coefficients instead of
the *negative* log likelihood! When fitting electrons this would sometimes
cause the fit to become unstable and continue increasing the path coefficients
without bound since the gain in the likelihood caused by increasing the
coefficients was more than the loss caused by a worse fit to the PMT data.
Doh!
|
|
|
|
This commit fixes a bug in the calculation of the average rms width of the
angular distribution for a path with a KL expansion. I also made a lot of
updates to the test-path program:
- plot the distribution of the KL expansion coefficients
- plot the standard deviation of the angular distribution as a function of
distance along with the prediction
- plot the simulated and reconstructed path in 3D
|
|
|
|
|
|
This commit adds a function called get_path_length() which computes the path
length inside and outside a sphere for a line segment between two points. This
will be useful for calculating the photon absorption for paths which cross the
AV and for computing the time of flight of photons from a track to a PMT.
|
|
|
|
This commit updates the calculation of the muon kinetic energy as a function of
distance along the track. Previously I was using an approximation from the PDG,
but it doesn't seem to be very accurate and won't generalize to the case of
electrons. The kinetic energy is now calculated using the tabulated values of
dE/dx as a function of energy.
|
|
This commit makes sure that when we conolve the single PE charge distribution
with a gaussian we integrate starting at zero since the PDF is zero for q < 0.
|
|
This commit adds the function ln() to compute log(n) for integer n. It uses a
lookup table for n < 100 to speed things up.
|
|
This commit adds a fast function to calculate the expected number of PE at a
PMT without numerically integrating over the track. This calculation is *much*
faster than integrating over the track (~30 ms compared to several seconds) and
so we use it during the "quick" minimization phase of the fit to quickly find
the best position.
|
|
For some reason the fit seems to have trouble with the kinetic energy.
Basically, it seems to "converge" even though when you run the minimization
again it finds a better minimum with a lower energy. I think this is likely due
to the fact that for muons the kinetic energy only really affects the range of
the muon and this is subject to error in the numerical integration.
I also thought that maybe it could be due to roundoff error in the likelihood
calculation, so I implemented the Kahan summation to try and reduce that. No
idea if it's actually improving things, but I should benchmark it later to see.
|
|
spaced
|
|
|
|
This commit updates the likelihood fit to use the KL path expansion. Currently,
I'm just using one coefficient for the path in both x and y.
|
|
|
|
To fit the path of muons and electrons I use the Karhunen-Loeve expansion of a
random 2D walk in the polar angle in x and y. This allows you to decompose the
path into a sum over sine functions whose coefficients become random variables.
The nice thing about fitting the path in this way is that you can capture
*most* of the variation in the path using a small number of variables by only
summing over the first N terms in the expansion and it is easy to calculate the
probability of the coefficients since they are all uncorrelated.
|
|
|
