Age | Commit message (Collapse) | Author |
|
|
|
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.
|
|
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 fixes a potential heap overflow in interp1d() which could occur if
x was very close to the last value in the xp array. The bounds check is now
performed on the index rather than the x values.
|
|
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 fixes a bug in interp1d() which would cause the function to read
past the end of the arrays if the interpolation point was equal to the last
point in the xp array.
|
|
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
|
|
|
|
|