aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
AgeCommit message (Collapse)Author
2018-09-13add a function to compute log(n) for integer ntlatorre
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.
2018-09-13speed things up by introducing a minimum ratio between probabilitiestlatorre
Previously to avoid computing P(q,t|n)*P(n|mu) for large n when they were very unlikely I was using a precomputed maximum n value based only on the expected number of PE. However, this didn't take into account P(q|n). This commit updates the likelihood function to dynamically decide when to quit computing these probabilities when the probability for a given n divided by the most likely probability is less than some threshold. This threshold is currently set to 10**(-10) which means we quit calculating these probabilities when the probability is 10 million times less likely than the most probable value.
2018-09-12small updates to speed things uptlatorre
2018-09-11update fast likelihood function to include the pmt response and absorptiontlatorre
2018-09-10add a fast likelihood functiontlatorre
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.
2018-09-04add a function to return the kahan sum of an arraytlatorre
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.
2018-08-31add epsrel argument to likelihood functiontlatorre
2018-08-31switch back to calling cquad just once to speed things uptlatorre
I found when simulating high energy muons that the expected charge for some PMTs which should be getting hit was zero. The reason for this is that the integrand was very sharply peaked at the Cerenkov angle which makes it difficult to integrate for numerical integration routines like cquad. To solve this I split up the integral at the point when the track was at the Cerenkov angle from the PMT to make sure that cquad didn't miss the peak. However, calling cquad twice takes a lot of time so it's not necessarily good to do this for all fits. Also, it's not obvious if it is necessary any more now that the angular distribution calculation was fixed. I think the real reason that cquad was missing those integrals was that for a high energy muon the range is going to be very large (approximately 40 meters for a 10 GeV muon). In this case, I should really only integrate up to the edge of the cavity or PSUP and hopefully cquad picks enough points in there to get a non zero value. I also added a check to only compute tmean when at least one PMT has a valid time. This prevents a divide by zero which causes the likelihood function to return nan.
2018-08-28add path to the likelihood fittlatorre
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.
2018-08-27update code to use get_index_snoman* functions to calculate the index of ↵tlatorre
refraction
2018-08-14fix how the RMS scattering angle is calculatedtlatorre
The RMS scattering angle calculation comes from Equation 33.15 in the PDG article on the passage of particles through matter. It's not entirely obvious if this equation is correct for a long track. It seems like it should be integrated along the track to add up the contributions at different energies, but it's not obvious how to do that with the log term. In any case, the way I was previously calculating it (by using the momentum and velocity at each point along the track) was definitely wrong. I will try this out and perhaps try to integrate it later.
2018-08-14move everything to src directorytlatorre