aboutsummaryrefslogtreecommitdiff
path: root/src/misc.h
AgeCommit message (Collapse)Author
2018-11-11update likelihood function to fit electrons!tlatorre
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.
2018-10-19add interp2d() for fast bilinear 2D interpolationtlatorre
2018-10-19update path integral to use a fixed number of pointstlatorre
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.
2018-10-18fix the likelihood function to return the *negative* log likelihood of the ↵tlatorre
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!
2018-10-17fix a bug in the theta0 calculation for a pathtlatorre
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
2018-09-17add get_path_length()tlatorre
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.
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-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 interp1d function to do fast interpolation when the x values are evenly ↵tlatorre
spaced
2018-08-31rotate and translate the path in path_init to speed things uptlatorre
2018-08-14move everything to src directorytlatorre