| Age | Commit message (Collapse) | Author |
|---|
|
|
|
This commit updates the code to calculate the number of Cerenkov photons from
secondary particles produced in an electromagnetic shower from electrons to use
an energy dependent formula I fit to data simulated with RAT-PAC.
|
|
This commit updates the charge likelihood calculation to calculate:
P(hit,q|n) = P(q|hit,n)*P(hit|n)
This has almost no effect on the fit results, but is technically correct.
|
|
This commit updates the optics code to calculate the rayleigh scattering length
using the Einstein-Smoluchowski formula instead of using the effective rayleigh
scattering lengths from the RSPR bank.
|
|
Thanks clang!
|
|
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).
|
|
This commit speeds up the likelihood function by integrating the charge along
the track inline instead of creating an array and then calling trapz(). It also
introduces two global variables avg_index_d2o and avg_index_h2o which are the
average indices of refraction for D2O and H2O weighted by the PMT quantum
efficiency and the Cerenkov spectrum.
|
|
This commit speeds up the likelihood calculation by eliminating most calls to
acos(). This is done by updating the PMT response lookup tables to be as a
function of the cosine of the angle between the photon and the PMT normal
instead of the angle itself.
|
|
Previously I was computing the fraction of light absorbed and scattered by
calculating an average absorption and scattering length weighted by the
Cerenkov spectrum and the PMT quantum efficiency, which isn't correct since we
should be averaging the absorption and scattering probabilities, not the
absorption and scattering lengths.
This commit fixes this by instead computing the average probability that a
photon is absorbed or scattered as a function of the distance travelled by
integrating the absorption and scattering probabilities over all wavelengths
weighted by the PMT quantum efficiency and the Cerenkov spectrum.
|
|
|
|
|
|
|
|
|
|
This is so that in the future if we only integrate over the path in the PSUP we
don't overestimate the Cerenkov light from delta rays.
|
|
This commit updates the likelihood function to take into account Cerenkov light
produced from delta rays produced by muons. The angular distribution of this
light is currently assumed to be constant along the track and parameterized in
the same way as the Cerenkov light from an electromagnetic shower. Currently I
assume the light is produced uniformly along the track which isn't exactly
correct, but should be good enough.
|
|
This commit adds a new function fit_event2() to fit multiple vertices. To seed
the fit, fit_event2() does the following:
- use the QUAD fitter to find the position and initial time of the event
- call find_peaks() to find possible directions for the particles
- loop over all possible unique combinations of the particles and direction
vectors and do a "fast" minimization
The best minimum found from the "fast" minimizations is then used to start the fit.
This commit has a few other updates:
- adds a hit_only parameter to the nll() function. This was necessary since
previously PMTs which weren't hit were always skipped for the fast
minimization, but when fitting for multiple vertices we need to include PMTs
which aren't hit since we float the energy.
- add the function guess_energy() to guess the energy of a particle given a
position and direction. This function estimates the energy by summing up the
QHS for all PMTs hit within the Cerenkov cone and dividing by 6.
- fixed a bug which caused the fit to freeze when hitting ctrl-c during the
fast minimization phase.
|
|
|
|
|
|
|
|
|
|
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.
|
|
|
|
integral
|
|
This commit speeds up the fast likelihood calculation by only computing the
time PDF for a single photon. Since the majority of the time in the fast
likelihood calculation is spent computing the time PDF this should speed things
up by quite a bit. I suspect this won't have a big effect on the likelihood
value, but I should do some more testing.
|
|
|
|
|
|
|
|
|
|
This function is only used when the expected number of photons reaching a PMT
is *very* small. In this case, we still need to estimate the PMT hit time PDF
for indirect light which is modelled as a flat distribution starting at the
time where the PMT is most likely to be hit from direct light. Since we compute
the most likely time for a PMT to be hit from direct light by computing the
integral of the expected charge times the time and then dividing by the total
charge, when the total charge is very small this can introduce large errors.
Note that this code already existed but it was computed in the likelihood
function. This commit just moves it to its own function to make things look
nicer.
|
|
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.
|
|
|
|
|
|
This commit speeds up the likelihood calculation by returning zero early if the
angle between the PMT and the track is far from the Cerenkov angle.
Specifically we check to see that the angle is 5 "standard deviations" away.
Where the standard deviation is taken to be the RMS width of the angular
distribution.
|
|
|
|
|
|
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.
|
|
|
|
|
|
|
|
|
|
things up
|
|
|
|
|
|
|
|
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.
|
|
This commit fixes a bug which was double counting the pmt response when
computing the direct charge and incorrectly multiplying the reflected charge by
the pmt response. I think this was just a typo left in when I added the
reflected charge.
|
|
Occasionally when fitting electrons the kinetic energy at the last step would
be high enough that the electron never crossed the BETA_MIN threshold which
would cause the gsl routine to throw an error.
This commit updates particle_init() to set the kinetic energy at the last
step to zero to make sure that we can bisect the point along the track where
the speed drops to BETA_MIN.
|
|
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.
|
