| Age | Commit message (Collapse) | Author |
|---|
|
|
|
|
|
The range and energy loss tables have different maximum values for electrons,
muons, and protons so we have to dynamically set the maximum energy of the fit
in order to avoid a GSL interpolation error.
This commit adds {electron,muon,proton}_get_max_energy() functions to return
the maximum energy in the tables and that is then used to set the maximum value
in the fit.
|
|
|
|
Also write out the data cleaning word to the YAML file.
|
|
Small update to the flasher cut to use the non-walk corrected time instead of
just the ECA calibrated time.
Also added some comments to the variables in the event structure.
|
|
|
|
This commit updates the get_event() function to flag any hits which have a
non-zero best charge status word. This essentially gets set when *both* QHS and
QLX are either bad or railed.
I also set this bit if *just* QHS is railed or below 300 since the current
charge model was only set up for QHS. In the future, it would be nice to
instead use the best charge (either QHS or QLX) in the likelihood function. To
do that I need to double check how QLX is normalized and if the current charge
model would work for QLX too.
|
|
|
|
Also increase the maximum kinetic energy to 10^4 GeV which is approximately the
maximum expected energy for cosmic muons at SNO.
|
|
Also, call this function when computing the psi parameter in nll_best().
|
|
This commit adds a data cleaning cut to tag incoming muons by looking for early
OWL hits. It also significantly updates the flasher cut to catch more flashers.
In particular, the flasher cut now does the following:
- loops over *all* paddle cards with at least 4 hits instead of just the paddle
cards with the most hits
- uses QLX to look for charge outliers in the paddle card
- fixes a few bugs (for example, uninitialized values in the charge array)
- adds a check to to see if the given slot is early with respect to all PMTs
within 4 meters to catch the case where the flashing channel is missing from
the event
|
|
This commit adds the following data quality cuts used in SNOMAN:
- neck
- qvnhit
- crate isotropy
- junk
Still need to test these.
|
|
This commit updates get_event() to include OWL, LG, FECD, BUTT, and NECK tubes.
|
|
|
|
|
|
|
|
No reason to use fast_sqrt() in get_theta0_min() since the argument isn't
guaranteed to be between 0 and 1.
|
|
|
|
This commit makes a few small changes to try and reduce the number of divisions
and multiplications done in get_expected_charge() to speed up the likelihood
function.
|
|
Since we already calculate sin(theta) in get_expected_charge() there's no
reason to calculate it again in get_probability(). This *may* already be
optimized out by the compiler.
|
|
|
|
|
|
This commit updates find_peaks() to only return peaks which are at least a
certain number of degrees apart from each other. This is because I found that
for many events the first few peaks would all be essentially the same direction
and so the fit was taking a lot of time fitting essentially the same seed
points. Since I now have to only try 3 peaks in order to get my grid jobs to
run for less than a few hours it's necessary to make sure we aren't just
fitting the same three directions for the "quick" minimization.
I also updated the fit to only use a maximum of 3 seed directions.
|
|
|
|
|
|
|
|
Also, update the step size for the energy during the final minimization to 10%.
|
|
|
|
|
|
This commit changes the format specifier for the values in sprintf_yaml_list()
from %.2g -> %.2f because YAML (at least the python parser) doesn't recognize
values like 1e+03 as floats.
|
|
Also, delete a check on cos_theta_pmt which isn't necessary.
|
|
- set number of shower points to 10 for the main fit
- set step size to 10% of the energy
- set max number of evals during quick minimization phase to 1000
|
|
|
|
I probably need to spend some time to optimize this along with the algorithm
for guessing the peaks, but for now I am just lowering this from 10 -> 5
because with 10 the number of quick minimizations for 3 particles is too big
and so the fits take way too long.
|
|
When plotting the likelihood function I realized that the fast likelihood
calculation was *very* noisy due to the way I calculated the shower and delta
ray charge. Although it works well for single particles, it is not suitable for
distinguishing which seed is the best when doing multi particle fits.
Eventually I may be able to fix this, but for now we just do the normal
likelihood calculation.
I also decreased the number of shower points from 100 -> 10 to speed things up.
|
|
Based on some testing it seems that when fitting muons the likelihood ratio and
angular fits are better without a minimum theta0. I also determined during
testing that the minimum value would cause a discontinuity in the derivative of
the charge as a function of position which could cause the estimate of the
direct charge to be worse.
|
|
|
|
This commit adds a new program called zdab-cat which is kind of like fit, but
just produces the YAML output without actually fitting anything.
|
|
|
|
|
|
Since we now calculate the expected charge from shower photons for muons we
need to initialize the angular distribution and a few other things in
particle_init().
|
|
Similarly to electrons, I fit an analytic form to the ratio of the number of
photons produced via shower particles over the radiative energy loss. In this
case, I chose the functional form:
ratio = a*(1-exp(-T/b))
since the ratio seemed to reach a constant value after a certain energy. I then
simulated a 10 GeV muon and it appears that the ratio might actually decrease
after that, so for higher energies I may have to come up with a different fit
function.
|
|
gsl error
|
|
|
|
|
|
delta rays
This commit introduces a new method for integrating over the particle track to
calculate the number of shower and delta ray photons expected at each PMT. The
reason for introducing a new method was that the previous method of just using
the trapezoidal rule was both inaccurate and not stable. By inaccurate I mean
that the trapezoidal rule was not producing a very good estimate of the true
integral and by not stable I mean that small changes in the fit parameters
(like theta and phi) could produce wildly different results. This meant that
the likelihood function was very noisy and was causing the minimizers to not be
able to find the global minimum.
The new integration method works *much* better than the trapezoidal rule for
the specific functions we are dealing with. The problem is essentially to
integrate the product of two functions over some interval, one of which is very
"peaky", i.e. we want to find:
\int f(x) g(x) dx
where f(x) is peaked around some region and g(x) is relatively smooth. For our
case, f(x) represents the angular distribution of the Cerenkov light and g(x)
represents the factors like solid angle, absorption, etc.
The technique I discovered was that you can approximate this integral via a
discrete sum:
constant \sum_i g(x_i)
where the x_i are chosen to have equal spacing along the range of the integral
of f(x), i.e.
x_i = F^(-1)(i*constant)
This new method produces likelihood functions which are *much* more smooth and
accurate than previously.
In addition, there are a few other fixes in this commit:
- switch from specifying a step size for the shower integration to a number of
points, i.e. dx_shower -> number of shower points
- only integrate to the PSUP
I realized that previously we were integrating to the end of the track even
if the particle left the PSUP, and that there was no code to deal with the
fact that light emitted beyond the PSUP can't make it back to the PMTs.
- only integrate to the Cerenkov threshold
When integrating over the particle track to calculate the expected number
of direct Cerenkov photons, we now only integrate the track up to the point
where the particle's velocity is 1/index. This should hopefully make the
likelihood smoother because previously the estimate would depend on exactly
whether the points we sampled the track were above or below this point.
- add a minimum theta0 value based on the angular width of the PMT
When calculating the expected number of Cerenkov photons we assumed that
the angular distribution was constant over the whole PMT. This is a bad
assumption when the particle is very close to the PMT. Really we should
average the function over all the angles of the PMT, but that would be too
computationally expensive so instead we just calculate a minimum theta0
value which depends on the distance and angle to the PMT. This seems to
make the likelihood much smoother for particles near the PSUP.
- add a factor of sin(theta) when checking if we can skip calculating the
charge in get_expected_charge()
- fix a nan in beta_root() when the momentum is negative
- update PSUP_RADIUS from 800 cm -> 840 cm
|
|
|
|
|
|
This commit updates the test-find-peaks script to plot Cerenkov rings for each
of the peaks. It also updates the script to use quad to find the position
instead of using the MC information. Finally, I added a -n argument to the
script to specify how many peaks to draw.
|
