| Age | Commit message (Collapse) | Author |
|---|
|
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 adds the probability that a channel is miscalibrated and/or doesn't
make it into the event to the likelihood. This was added because I noticed when
looking at the likelihood for one very high energy event that there was a
single PMT that should have been hit that wasn't in the event and which was not
marked as bad in DQXX.
I did some testing and the addition of this term does not seem to significantly
affect that atmospheric MC or the psi values for flashers. One unexpected
improvement is that it seems that external muons are more likely to correctly
reconstruct at the PSUP with this change. I haven't determined the exact cause
but I suspect it's because there is some mismodelling of the likelihood for
muons near the edge of the detector when they exit and that adding this term
allows the likelihood to ignore these PMT hits.
|
|
This commit adds four scripts:
1. calculate-atmospheric-oscillations
This script uses an independent python package called nucraft to calculate the
neutrino oscillation probabilities for atmospheric neutrinos. These
probabilities are calculated as a function of energy and cosine of the zenith
angle and stored in text files.
2. apply-atmospheric-oscillations
This script combines the high energy 2D atmospheric neutrino flux from Barr and
the low energy 1D flux from Battistoni and then applies neutrino oscillations
to them. The results are then stored in new flux files that can be used with a
modified version of GENIE to simulate the oscillated atmospheric neutrino flux
from 10 MeV to 10 GeV.
3. plot-atmospheric-fluxes
This is a simple script to plot the atmospheric flux files produced by
apply-atmospheric-oscillations.
4. plot-atmospheric-oscillations
This is a simple script to plot the 2D neutrino oscillation probabilities.
|
|
electron showers
Also update the a parameter based on a simple 0 degree polynomial fit to the
shower profiles above 100 MeV.
|
|
|
|
|
|
This commit adds the --save command line argument to plot-energy to save either
the corner plots or the energy distribution plots. It also updates the code to
make plots similar to plot-fit-results.
In addition there are a bunch of other small changes:
- plot the theoretical Michel spectrum for free muons
- energy plots now assume there are only a max of 2 particles fit for each event
- display particle IDs as letters instead of numbers, i.e. 2022 -> eu
|
|
|
|
|
|
the database
Also add a -r or --reprocess command line option to reprocess runs which are
already in the database.
|
|
results
|
|
|
|
To select stopping muons we simply look for the muons before a Michel event.
The muon distance is calculated by first projecting the muon fit back to the
PSUP along the fitted direction and then taking the distance between this point
and the fitted position of the Michel event. I then calculate the expected
kinetic energy of the muon by using the muon lookup tables of the CSDA range to
convert the distance to an energy.
I also changed a few other things like changing as_index=False ->
group_keys=False when grouping events. The reason for this is just that if we
do this we don't have to reset the index and drop the new index after calling
apply().
I also fixed a small bug I had introduced recently where I selected only prompt
events before finding the michel events and atmospheric events.
This commit updates the plot-energy script to plot the energy bias and resolution for stopping muons by computing the expected
|
|
This commit updates the dc script to calculate the instrumental contamination
to now treat all 4 high level variables as correlated for muons. Previously I
had assumed that the reconstructed radius was independent from udotr, z, and
psi, but based on the corner plots it seems like the radius is strongly
correlated with udotr.
I also updated the plotting code when using the save command line argument to
be similar to plot-fit-results.
|
|
This commit adds a script to calculate the background contamination using a
method inspired by the bifurcated analysis method used in SNO. The method works
by looking at the distribution of several high level variables (radius, udotr,
psi, and reconstructed z position) for events tagged by the different data
cleaning cuts and assuming that any background events which sneak past the data
cleaning cuts will have a similar distribution (for certain backgrounds this is
assumed and for others I will actually test this assumption. For more details
see the unidoc). Then, by looking at the distribution of these high level
variables for all the untagged events we can use a maximum likelihood fit to
determine the residual contamination.
There are also a few other updates to the plot-energy script:
- add a --dc command line argument to plot corner plots for the high level
variables used in the contamination analysis
- add a fudge factor to the Ockham factor of 100 per extra particle
- fix a bug by correctly setting the final kinetic energy to the sum of the
individual kinetic energies instead of just the first particle
- fix calculation of prompt events by applying at the run level
|
|
|
|
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.
|
|
|
|
|
|
|
|
|
|
|
|
is large
Previously to achieve a large speedup in the likelihood calculation I added a
line to skip calculating the charge if:
abs((cos(theta)-cos_theta_cerenkov)/(sin_theta*theta0)) > 5
However I noticed that this was causing discontinuities in the likelihood
function when fitting low energy muons so I'm putting it behind a compile time
flag for now.
|
|
This commit updates plot-fit-results to use the median when plotting the energy
and position bias and the interquartile range (times 1.35) when plotting the
energy and position resolution. The reason is that single large outliers for
higher energy muons were causing the energy bias and resolution to no longer
represent the central part of the distribution well.
|
|
This commit updates how we handle PMTs whose type is different in the
snoman.ratdb file and the SNOMAN bank again. In particular, we now trust the
snoman.ratdb type *only* for the NCD runs and mark the PMT as invalid for the
D2O and salt phases.
This was spurred by noticing that with the current code GTID 9228 in run 10,000
was being marked as a neck event even though it was clearly a muon and XSNOED
only showed one neck hit. It was marked as a neck event because there were 2
neck PMT hits in the event: 3/15/9 and 13/15/0. After reading Stan's email more
carefully I realized that 3/15/9 was only installed as a neck PMT in the NCD
phase. I don't really know what type of PMT it was in the D2O and salt phases
(maybe an OWL), but in any case since I don't know the PMT position I don't
think we can use this PMT for these phases.
|
|
This commit fixes a small bug in cat-grid-jobs which was causing it to print
the wrong filename when there was no git_sha1 attrs in the HDF5 file.
|
|
I noticed that many of my jobs were failing with the following error:
module: command not found
My submit description files *should* only be selecting nodes with modules because of this line:
requirements = (HAS_MODULES =?= true) && (OSGVO_OS_STRING == "RHEL 7") && (OpSys == "LINUX")
which I think I got from
https://support.opensciencegrid.org/support/solutions/articles/12000048518-accessing-software-using-distributed-environment-modules.
I looked up what the =?= operator does and it's a case sensitive search. I also
found another site
(https://support.opensciencegrid.org/support/solutions/articles/5000633467-steer-your-jobs-with-htcondor-job-requirements)
which uses the normal == operator. Therefore, I'm going to switch to the ==
operator and hope that fixes the issue.
|
|
This commit updates get_expected_charge() to always use the index of refraction
for d2o instead of choosing the index of d2o or h2o based on the position of
the particle. The reason for this is that selecting the index based on the
position was causing discontinuities in the likelihood function for muon tracks
which crossed the AV.
|
|
This commit fixes a bug in gen-dark matter which was causing the positions to
all be generated with positive x, y, and z values. Doh!
|
|
This commit adds a new test to test the quad fitter when the t0 quantile
argument is less than 1.
|
|
This commit updates get_event() to clear any PMT flags except for PMT_FLAG_DQXX
from all PMT hits before loading the event. Although I *was* previously
clearing the other flags for hit PMTs, I was not clearing flags for PMTs which
were *not* hit. This was causing non deterministic behaviour, i.e. I was
getting different results depending on if I ran the fitter over a whole file or
just a single event.
|
|
This commit updates the likelihood function to use the PMT hit time without the
time walk correction applied (when the charge is greater than 1.5 PE) instead
of the multiphoton PCA time. The reason is that after talking with Chris Kyba I
realized that the multiphoton PCA time was calibrated to give the mean PMT hit
time when mulitiple photons hit at the same time instead of the time when the
first photon hits which is what I assume in my likelihood function.
Therefore I now use the regular PMT hit time without time walk correction
applied which should be closer to the first order statistic.
|
|
|
|
This commit updates the likelihood function to initialize mu_indirect to 0.0
since it's a static array. This can have an impact when the fit position is
outside of the PSUP and we skip calculating the charges.
|
|
|
|
|
|
This commit updates the crate, card, and channel variables in is_slot_early()
to be ints instead of size_ts. The reason is I saw a warning when building with
clang and realized that the abs(card - flasher_card) == 1 check wasn't working
if flasher_card was 1 greater than card because I was using unsigned ints.
|
|
|
|
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.
|
|
|
|
|
|
|
|
|
|
|
