Key: - To do p Partially done ^ Done ? Not sure if it's done - write cross section section in paper ? explain integral in paper - not sure what integral I was talking about here. I just checked the git history and this was before I had written the fitter, so it must be talking about the integral to compute the expected event rate. p develop signal generator - made signal generator for two particle decays p develop high energy and multi-ring fitter ^ write talk for Monday ^ write a section on atmospheric backgrounds in paper ^ check index of refraction calculation against tables in the paper ^ write more tests for the solid angle calculation - use table from paper ^ make sure normal direction is normalized ^ speed up likelihood calculation - rotate and translate path points in path_init - create fast interp1d for equally spaced points - optimize with -O2 - only calculate average time if mu_direct[i] is not too small and PMT is hit - add likelihood function which calculates total expected charge and only loops over hit PMTs - add fast likelihood which doesn't use a 2D lookup table - add function to compute the Hessian via finite differences - add function to calculate covariance matrix - add function to compute the approximate volume in the likelihood space where the likelihood is greater than some value ^ add output to likelihood ^ add a set of seed points to start the likelihood ^ minimize the likelihood in a do while loop until it converges to a single point ^ once the first position is found, compute the Hough transform to look for other rings p add rayleigh scattering to likelihood function - Done, but it's pretty hacky p add muon shower PDF? delta rays? - Done for delta rays, but need to add shower photons - calculate path length through acrylic? - did some tests on this and it seems to fix the radial energy bias for electrons, so it's definitely an issue. However, currently the small 5% energy bias is not really a big deal, although I would like to fix this eventually. p update estar table with energies above 1 GeV - added values up to 10 GeV ^ speed up normalize() and interp1d() ^ add shower PDF to get_total_charge_approx() x compute mean of starting points ^ speed up time PDF calculation for fast likelihood p calculate psi parameter as a goodness of fit test - Still need to tweak this. It seems to have a strong nhit dependence ^ add a maximum length for the range - now only fit to the end of the PSUP ^ add a time limit for the fit ^ multi-particle fit p optimize CHARGE_FRACTION ? find out why likelihood is sometimes returning nan ^ add pdf for cerenkov light from delta rays ^ tweak find peaks algorithm ^ update zebra to allow use of links p minimum energy for each particle type - I do already have a minimum energy which is the Cerenkov threshold, however I think what I meant here was to have something significantly higher than that since a particle just at the Cerenkov threshold won't produce many PMT hits. I think I was mostly concerned about this issue because many fits were always fitting better by including a low energy electron, however I have solved that problem by simply dropping multi particle fits which have an electron with an energy less than 20 MeV. ^ check all zebra links ^ fix zero logical record size bug in zebra.c p speed up fit - have optimized most of the likelihood function. I think further optimization will require a different technique for evaluating the integrals ^ nhit cut p add flasher cut p default max particles - submit-grid-jobs fits up to 3 particles by default ^ add environment variable to control where to look for pmt.txt and titles files p update charge code to handle highest values for qhs and qlx - don't have this coded into the charge probability code, but I now use QHL if QHS is railed and flag any PMT with a railed QHL ^ use qhs and qlx based on variable in PMT bank - Note: I don't use the best charge variable in the PMT bank, but I do select QHL when QHS is railed ^ fix bug in path_init() ^ load DQXX file based on run number ^ update rayleigh scattering lengths ^ fix delta ray calculation p figure out why electron energy bias is +10% - updated code to use a more accurate approximation for the number of shower photons which seems to bring the bias down to 5% p update find peaks algorithm to do single particle quick fits - I did update the find peaks algorithm to work *much* better, but it still doesn't do single particle quick fits. I think the next improvement would be to add the ability to actually determine the number of rings. ^ figure out how to combine SNO fitter data with my fitter for final analysis. For example, how to apply neutron follower cut? - double check that zdab-reprocess is correct for the D2O and salt phases since it appears to be from the NCD phase ? add code to compute expected deviation from nll_best to normalize psi - tried several different versions of this and nothing seemed to perform as well as psi/nhit. ^ add term to likelihood for probability that a channel is misccalibrated or didn't make it into the event - when calculating the first order statistic for a gaussian is alpha = pi/8 a better choice than 0.375? - speed up charge initialization - determine *real* mean and standard deviation of single PE charge distribution. TELLIE SNO+ data? ^ extend electron range table - extended the electron range table up to 1 TeV ? fix time PDF. Currently take the first order statistic of dark noise + indirect light + direct light all together, but this isn't correct. - thought more about this and it's not actually obvious if what I'm doing is wrong. Need to think more about this. olor: #008800; font-weight: bold } /* Keyword.Declaration */ .highlight .kn { color: #008800; font-weight: bold } /* Keyword.Namespace */ .highlight .kp { color: #008800 } /* Keyword.Pseudo */ .highlight .kr { color: #008800; font-weight: bold } /* Keyword.Reserved */ .highlight .kt { color: #888888; font-weight: bold } /* Keyword.Type */ .highlight .m { color: #0000DD; font-weight: bold } /* Literal.Number */ .highlight .s { color: #dd2200; background-color: #fff0f0 } /* Literal.String */ .highlight .na { color: #336699 } /* Name.Attribute */ .highlight .nb { color: #003388 } /* Name.Builtin */ .highlight .nc { color: #bb0066; font-weight: bold } /* Name.Class */ .highlight .no { color: #003366; font-weight: bold } /* Name.Constant */ .highlight .nd { color: #555555 } /* Name.Decorator */ .highlight .ne { color: #bb0066; font-weight: bold } /* Name.Exception */ .highlight .nf { color: #0066bb; font-weight: bold } /* Name.Function */ .highlight .nl { color: #336699; font-style: italic } /* Name.Label */ .highlight .nn { color: #bb0066; font-weight: bold } /* Name.Namespace */ .highlight .py { color: #336699; font-weight: bold } /* Name.Property */ .highlight .nt { color: #bb0066; font-weight: bold } /* Name.Tag */ .highlight .nv { color: #336699 } /* Name.Variable */ .highlight .ow { color: #008800 } /* Operator.Word */ .highlight .w { color: #bbbbbb } /* Text.Whitespace */ .highlight .mb { color: #0000DD; font-weight: bold } /* Literal.Number.Bin */ .highlight .mf { color: #0000DD; font-weight: bold } /* Literal.Number.Float */ .highlight .mh { color: #0000DD; font-weight: bold } /* Literal.Number.Hex */ .highlight .mi { color: #0000DD; font-weight: bold } /* Literal.Number.Integer */ .highlight .mo { color: #0000DD; font-weight: bold } /* Literal.Number.Oct */ .highlight .sa { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Affix */ .highlight .sb { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Backtick */ .highlight .sc { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Char */ .highlight .dl { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Delimiter */ .highlight .sd { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Doc */ .highlight .s2 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Double */ .highlight .se { color: #0044dd; background-color: #fff0f0 } /* Literal.String.Escape */ .highlight .sh { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Heredoc */ .highlight .si { color: #3333bb; background-color: #fff0f0 } /* Literal.String.Interpol */ .highlight .sx { color: #22bb22; background-color: #f0fff0 } /* Literal.String.Other */ .highlight .sr { color: #008800; background-color: #fff0ff } /* Literal.String.Regex */ .highlight .s1 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Single */ .highlight .ss { color: #aa6600; background-color: #fff0f0 } /* Literal.String.Symbol */ .highlight .bp { color: #003388 } /* Name.Builtin.Pseudo */ .highlight .fm { color: #0066bb; font-weight: bold } /* Name.Function.Magic */ .highlight .vc { color: #336699 } /* Name.Variable.Class */ .highlight .vg { color: #dd7700 } /* Name.Variable.Global */ .highlight .vi { color: #3333bb } /* Name.Variable.Instance */ .highlight .vm { color: #336699 } /* Name.Variable.Magic */ .highlight .il { color: #0000DD; font-weight: bold } /* Literal.Number.Integer.Long */
#!/usr/bin/env python
from __future__ import print_function, division
import subprocess
from os.path import splitext, split

def run_fit(filename):
    head, tail = split(filename)
    root, ext = splitext(tail)
    output = root + '.hdf5'
    cmd = ["./fit", filename, "-o", output]
    subprocess.call(cmd)

if __name__ == '__main__':
    import argparse
    from multiprocessing import Pool, cpu_count
    import signal
    import os

    parser = argparse.ArgumentParser("fit multiple zdab files")
    parser.add_argument("-j", "--jobs", type=int, default=None,
                        help="number of jobs")
    parser.add_argument("filenames", nargs="+",
                        help="zdab files")
    args = parser.parse_args()

    jobs = args.jobs

    if jobs is None:
        jobs = cpu_count()

    # see https://stackoverflow.com/questions/11312525/catch-ctrlc-sigint-and-exit-multiprocesses-gracefully-in-python
    handler = signal.signal(signal.SIGINT, signal.SIG_IGN)
    p = Pool(jobs)
    signal.signal(signal.SIGINT, handler)
    try:
        result = p.map_async(run_fit,args.filenames)
        result.get()
    except KeyboardInterrupt:
        print("ctrl-c caught. quitting...")
        p.terminate()
        os.killpg(os.getpgid(0),9)
    else:
        p.close()

    p.join()