| Age | Commit message (Collapse) | Author |
|---|
|
This commit updates both the PMT response and optics code to be able to
load constants from the salt phase.
|
|
|
|
This commit updates fit.c to start with 5 peaks for the direction seeds. I
chose this number because I did some testing with the test-find-peaks program
on the atmospheric MC and it looks like 5 peaks were necessary to capture the
majority of the peaks.
|
|
This commit updates the fit program to fit each event and particle hypothesis
twice, once using the normal quad implementation and the other by cutting on
the 10% quantile of times. The first way is much much better when the event is
fully contained since quad will return a really good starting point, and the
second is much better for muons where we want to seed the fit near the entry
point of the muon.
Ideally we would only need a single call and I have an idea of how to update
QUAD to maybe return reasonable guesses in both cases. The idea is to take the
cloud of quad points and find the position and time that has the smallest time
such that it is only a certain Mahalabonis distance from the distribution. This
(I think) corresponds roughly to what I would do by eye where you look at the
distribution of quad points in the cloud and see that it forms a track, and
pick a point at the start of the track.
I started working on this second idea but haven't successfully tested it out
yet.
|
|
This commit updates the find peaks algorithm with several improvements which
together drastically improve its ability to find Cerenkov rings:
- when computing the Hough transform, instead of charge we weight each PMT hit
by the probability that it is a multi-photon PMT hit
- we don't subtract off previously found rings (this makes the code simpler and
I don't think it previously had a huge effect)
- ignore PMT hits who are within approximately 5 degrees of any previously
found ring (previously we ignored all hits within the center of previously
found rings)
- ignore PMT hits which have a time residual of more than 10 nanoseconds to
hopefully ignore more reflected and/or scattered light
- switch from weighting the Hough transform by exp(-fabs(cos(theta)-1/n)/0.1)
-> exp(-pow(cos(theta)-1/n,2)/0.01). I'm still not sure if this has a huge
effect, but the reason I switched is that the PDF for Cerenkov light looks
closer to the second form.
- switch to calling quad with f = 1.0 in test-find-peaks (I still need to add
this update to fit.c but will do that in a later commit).
|
|
|
|
This commit updates the ./fit program to add a ctrl-z handler to allow you to
skip events. This is really handy when testing nwe things. Currently if you
press ctrl-z and it has already done at least one of the initial fits, it will
skip to move on to the final minimization stage. If you press ctrl-z during the
final minimization, it will skip fitting the event. Currently this will *not*
save the result to the file but I may change that in the future.
|
|
This commit updates get_expected_photons() to check if there are any shower
photons or delta ray photons before adding them since if there aren't any
shower photons or delta ray photons the PDF isn't constructed.
|
|
PSUP
|
|
This commit updates submit-grid-jobs so that it keeps a database of jobs. This
allows the script to make sure that we only have a certain number of jobs in
the job queue at a single time and automatically resubmitting failed jobs. The
idea is that it can now be run once to add jobs to the database:
$ submit-grid-jobs ~/zdabs/SNOCR_0000010000_000_p4_reduced.xzdab.gz
and then be run periodically via crontab:
PATH=/usr/bin:$HOME/local/bin
SDDM_DATA=$HOME/sddm/src
DQXX_DIR=$HOME/dqxx
0 * * * * submit-grid-jobs --auto --logfile ~/submit.log
Similarly I updated cat-grid-jobs so that it uses the same database and can
also be run via a cron job:
PATH=/usr/bin:$HOME/local/bin
SDDM_DATA=$HOME/sddm/src
DQXX_DIR=$HOME/dqxx
0 * * * * cat-grid-jobs --logfile cat.log --output-dir $HOME/fit_results
I also updated fit so that it keeps track of the total time elapsed including
the initial fits instead of just counting the final fits.
|
|
|
|
|
|
get_hough_transform()
This commit adds two improvements to the quad fitter:
1. I updated quad to weight the random PMT hit selection by the probability
that the PMT hit is a multiphoton hit. The idea here is that we really only
want to sample direct light and for high energy events the reflected and
scattered light is usually single photon.
2. I added an option to quad to only use points in the quad cloud which are
below a given quantile of t0. The idea here is that for particles like muons
which travel more than a few centimeters in the detector the quad cloud usually
looks like the whole track. Since we want the QUAD fitter to find the position
of the *start* of the track we select only those quad cloud points with an
early time so the position is closer to the position of the start of the track.
Also, I fixed a major bug in get_hough_transform() in which I was using the
wrong index variable when checking if a PMT was not flagged, a normal PMT, and
was hit. This was causing the algorithm to completely miss finding more than
one ring while I was testing it.
|
|
This commit updates guess_energy() which is used to seed the energy for the
likelihood fit. Previously we estimated the energy by summing up the charge in
a 42 degree cone around the proposed direction and then dividing that by 6
(since electrons in SNO and SNO+ produce approximately 6 hits/MeV). Now,
guess_energy() estimates the energy by calculating the expected number of
photons produced from Cerenkov light, EM showers, and delta rays for a given
particle at a given energy. The most likely energy is found by bisecting the
difference between the expected number of photons and the observed charge to
find when they are equal.
This improves things dramatically for cosmic muons which have energies of ~200
GeV. Previously the initial guess was always very low (~1 GeV) and the fit
could take > 1 hour to increase the energy.
|
|
|
|
This commit updates the zebra code to store a pointer to the first MAST bank in
the zebraFile struct so that we can jump to it when iterating over the logical
records. I had naively assumed based on the documenation in the SNOMAN
companion that the first bank in a logical record was guaranteed to be a MAST
bank, but that doesn't seem to be the case. This also explains why I was
sometimes seeing RHDR and ZDAB banks as the first bank in a logical record.
|
|
This commit adds the sub_run variable to the ev array in the HDF5 output file
and updates plot-energy to order the events using the run and sub_run
variables. This fixes a potential issue where I was sorting by GTID before, but
the GTID can wrap around and so isn't guaranteed to put the events in the right
order.
|
|
|
|
|
|
This commit updates the zebra code to properly handle all the errors from
get_bytes(). I also updated fit and cat-zdab to not display the errors about
the FTX banks by default unless you run them with the -v command line option.
|
|
|
|
This commit updates the fit program to accept a particle combo from the command
line so you can fit for a single particle combination hypothesis. For example
running:
$ ./fit ~/zdabs/mu_minus_700_1000.hdf5 -p 2020
would just fit for the 2 electron hypothesis.
The reason for adding this ability is that my grid jobs were getting evicted
when fitting muons in run 10,000 since it takes 10s of hours to fit for all the
particle hypothesis. With this change, and a small update to the
submit-grid-jobs script we now submit a single grid job per particle
combination hypothesis which should make each grid job run approximately 4
times faster.
|
|
|
|
|
|
This commit fixes a bug introduced in a previous commit when I moved all the
code to add the run, gtid, etc. to the event object to get_event(). We actually
need the run number before then to load the DQXX file for the right run.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
The first logical record in the SNOCR files don't have an EV bank which was
causing the output file to have an emtpy list element. This commit fixes the
issue by checking for an empty EV bank before printing the list delimiter.
|
|
|
|
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.
|
|
Also increase the maximum kinetic energy to 10^4 GeV which is approximately the
maximum expected energy for cosmic muons at SNO.
|
|
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 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.
|
|
- 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.
|
|
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.
|
|
|
|
|
|
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
|
