#!/usr/bin/env python # Copyright (c) 2019, Anthony Latorre # # This program is free software: you can redistribute it and/or modify it # under the terms of the GNU General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) # any later version. # # This program is distributed in the hope that it will be useful, but WITHOUT # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for # more details. # # You should have received a copy of the GNU General Public License along with # this program. If not, see . """ Script to do final null hypothesis test that the events in the 20 MeV - 10 GeV range are consistent with atmospheric neutrino events. To run it just run: $ ./chi2 [list of data fit results] --mc [list of atmospheric MC files] --muon-mc [list of muon MC files] --steps [steps] After running you will get a plot showing the best fit results of the MC and data along with p-values for each of the possible particle combinations (single electron, single muon, double electron, etc.). """ from __future__ import print_function, division import numpy as np from scipy.stats import iqr, poisson from matplotlib.lines import Line2D from scipy.stats import iqr, norm, beta, percentileofscore from scipy.special import spence from itertools import izip_longest from sddm.stats import * from sddm.dc import estimate_errors, EPSILON, truncnorm_scaled import emcee from sddm import printoptions from sddm.utils import fast_cdf, correct_energy_bias # Likelihood Fit Parameters # 0 - Atmospheric Neutrino Flux Scale # 1 - Electron energy bias # 2 - Electron energy resolution # 3 - Muon energy bias # 4 - Muon energy resolution # 5 - External Muon scale FIT_PARS = [ 'Atmospheric Neutrino Flux Scale', 'Electron energy bias', 'Electron energy resolution', 'Muon energy bias', 'Muon energy resolution', 'External Muon scale'] # Uncertainty on the energy scale # # - the muon energy scale and resolution terms come directly from measurements # on stopping muons, so those are known well. # - for electrons, we only have Michel electrons at the low end of our energy # range, and therefore we don't really have any good way of constraining the # energy scale or resolution. However, if we assume that the ~7% energy bias # in the muons is from the single PE distribution (it seems likely to me that # that is a major part of the bias), then the energy scale should be roughly # the same. Since the Michel electron distributions are consistent, we leave # the mean value at 0, but to be conservative, we set the error to 10%. # - The energy resolution for muons was pretty much spot on, and so we expect # the same from electrons. In addition, the Michel spectrum is consistent so # at that energy level we don't see anything which leads us to expect a major # difference. To be conservative, and because I don't think it really affects # the analysis at all, I'll leave the uncertainty here at 10% anyways. # - For the electron + muon energy resolution, I don't have any real good way # to estimate this, so we are conservative and set the error to 10%. PRIORS = [ 1.0, # Atmospheric Neutrino Scale 0.0, # Electron energy scale 0.0, # Electron energy resolution 0.066, # Muon energy scale 0.0, # Muon energy resolution 0.0, # Muon scale ] PRIOR_UNCERTAINTIES = [ 0.2 , # Atmospheric Neutrino Scale 0.1, # Electron energy scale 0.1, # Electron energy resolution 0.011, # Muon energy scale 0.014, # Muon energy resolution 10.0, # Muon scale ] # Lower bounds for the fit parameters PRIORS_LOW = [ EPSILON, -10, EPSILON, -10, EPSILON, EPSILON ] # Upper bounds for the fit parameters PRIORS_HIGH = [ 10, 10, 10, 10, 10, 10 ] particle_id = {20: 'e', 22: r'\mu'} def plot_hist2(hists, bins, color=None): for id in (20,22,2020,2022,2222): if id == 20: plt.subplot(2,3,1) elif id == 22: plt.subplot(2,3,2) elif id == 2020: plt.subplot(2,3,4) elif id == 2022: plt.subplot(2,3,5) elif id == 2222: plt.subplot(2,3,6) bincenters = (bins[1:] + bins[:-1])/2 plt.hist(bincenters, bins=bins, histtype='step', weights=hists[id],color=color) plt.gca().set_xscale("log") plt.xlabel("Energy (MeV)") plt.title('$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(id),2)]) + '$') if len(hists): plt.tight_layout() def get_mc_hists(data,x,bins,scale=1.0,reweight=False): """ Returns the expected Monte Carlo histograms for the atmospheric neutrino background. Args: - data: pandas dataframe of the Monte Carlo events - x: fit parameters - bins: histogram bins - scale: multiply histograms by an overall scale factor This function does two basic things: 1. apply the energy bias and resolution corrections 2. histogram the results Returns a dictionary mapping particle id combo -> histogram. """ df_dict = {} for id in (20,22,2020,2022,2222): df_dict[id] = data[data.id == id] return get_mc_hists_fast(df_dict,x,bins,scale,reweight) def get_mc_hists_fast(df_dict,x,bins,scale=1.0,reweight=False): """ Same as get_mc_hists() but the first argument is a dictionary mapping particle id -> dataframe. This is much faster than selecting the events from the dataframe every time. """ mc_hists = {} for id in (20,22,2020,2022,2222): df = df_dict[id] if id == 20: ke = df.energy1.values*(1+x[1]) resolution = df.energy1.values*max(EPSILON,x[2]) elif id == 2020: ke = df.energy1.values*(1+x[1]) + df.energy2.values*(1+x[1]) resolution = np.sqrt((df.energy1.values*max(EPSILON,x[2]))**2 + (df.energy2.values*max(EPSILON,x[2]))**2) elif id == 22: ke = df.energy1.values*(1+x[3]) resolution = df.energy1.values*max(EPSILON,x[4]) elif id == 2222: ke = df.energy1.values*(1+x[3]) + df.energy2.values*(1+x[3]) resolution = np.sqrt((df.energy1.values*max(EPSILON,x[4]))**2 + (df.energy2.values*max(EPSILON,x[4]))**2) elif id == 2022: ke = df.energy1.values*(1+x[1]) + df.energy2.values*(1+x[3]) resolution = np.sqrt((df.energy1.values*max(EPSILON,x[2]))**2 + (df.energy2.values*max(EPSILON,x[4]))**2) if reweight: cdf = fast_cdf(bins[:,np.newaxis],ke,resolution)*df.weight.values else: cdf = fast_cdf(bins[:,np.newaxis],ke,resolution) mc_hists[id] = np.sum(cdf[1:] - cdf[:-1],axis=-1) mc_hists[id] *= scale return mc_hists def get_data_hists(data,bins,scale=1.0): """ Returns the data histogrammed into `bins`. """ data_hists = {} for id in (20,22,2020,2022,2222): data_hists[id] = np.histogram(data[data.id == id].ke.values,bins=bins)[0]*scale return data_hists def make_nll(data, muons, mc, atmo_scale_factor, muon_scale_factor, bins, print_nll=False): df_dict = {} for id in (20,22,2020,2022,2222): df_dict[id] = mc[mc.id == id] df_dict_muon = {} for id in (20,22,2020,2022,2222): df_dict_muon[id] = muons[muons.id == id] data_hists = get_data_hists(data,bins) def nll(x, grad=None): if any(x[i] < 0 for i in (0,2,4,5)): return np.inf # Get the Monte Carlo histograms. We need to do this within the # likelihood function since we apply the energy resolution parameters # to the Monte Carlo. mc_hists = get_mc_hists_fast(df_dict,x,bins,scale=1/atmo_scale_factor) muon_hists = get_mc_hists_fast(df_dict_muon,x,bins,scale=1/muon_scale_factor) # Calculate the negative log of the likelihood of observing the data # given the fit parameters nll = 0 for id in data_hists: oi = data_hists[id] ei = mc_hists[id]*x[0] + muon_hists[id]*x[5] + EPSILON N = ei.sum() nll -= -N - np.sum(gammaln(oi+1)) + np.sum(oi*np.log(ei)) # Add the priors nll -= norm.logpdf(x,PRIORS,PRIOR_UNCERTAINTIES).sum() if print_nll: # Print the result print("nll = %.2f" % nll) return nll return nll def get_mc_hists_posterior(data_mc,muon_hists,data_hists,atmo_scale_factor,muon_scale_factor,x,bins): """ Returns the posterior on the Monte Carlo histograms. Basically this function just histograms the Monte Carlo data. However, there is one extra thing it does. In general when doing a fit, the Monte Carlo histograms have some uncertainty since you can never simulate an infinite number of statistics. I don't think I've ever really seen anyone properly treat this. Since the uncertainty on the central value in each bin is just given by the Dirichlet distribution, we treat the problem of finding the best value of the posterior as a problem in which you're prior is equal to the expected number of events from the Monte Carlo, and then you actually see the data. Since the likelihood on the true mean in each bin is a multinomial, the posterior is also a dirichlet where the alpha parameters are given by a sum of the prior and observed counts. All that is a long way of saying we calculate the posterior as the sum of the Monte Carlo events (unscaled) and the observed events. In the limit of infinite statistics, this is just equal to the Monte Carlo predicted histogram, but deals with the fact that we don't have infinite statistics, and so a single outlier event isn't necessarily a problem with the model. Returns a dictionary mapping particle id combo -> histogram. """ mc_hists = get_mc_hists(data_mc,x,bins,reweight=True) for id in (20,22,2020,2022,2222): mc_hists[id] = get_mc_hist_posterior(mc_hists[id],data_hists[id],norm=x[0]/atmo_scale_factor) # FIXME: does the orering of when we add the muons matter here? mc_hists[id] += muon_hists[id]*x[5]/muon_scale_factor return mc_hists def get_multinomial_prob(data, data_muon, data_mc, weights, atmo_scale_factor, muon_scale_factor, id, x_samples, bins, percentile=50.0, size=10000): """ Returns the p-value that the histogram of the data is drawn from the MC histogram. The p-value is calculated by first sampling the posterior of the fit parameters `size` times. For each iteration we calculate a p-value. We then return the `percentile` percentile of all the p-values. This approach is similar to both the supremum and posterior predictive methods of calculating a p-value. For more information on different methods of calculating p-values see https://cds.cern.ch/record/1099967/files/p23.pdf. Arguments: data: 1D array of KE values data_mc: 1D array of MC KE values x_samples: MCMC samples of the floated parameters in the fit bins: bins used to bin the mc histogram size: number of values to compute """ df_dict_muon = {} for _id in (20,22,2020,2022,2222): df_dict_muon[_id] = data_muon[data_muon.id == _id] data_hists = get_data_hists(data,bins) # Get the total number of "universes" simulated in the GENIE reweight tool if len(data_mc): nuniverses = weights['universe'].max()+1 else: nuniverses = 0 ps = [] for i in range(size): x = x_samples[np.random.randint(x_samples.shape[0])] muon_hists = get_mc_hists_fast(df_dict_muon,x,bins) if nuniverses > 0: universe = np.random.randint(nuniverses) data_mc_with_weights = pd.merge(data_mc,weights[weights.universe == universe],how='left',on=['run','evn']) data_mc_with_weights.weight = data_mc_with_weights.weight.fillna(1.0) else: data_mc_with_weights = data_mc.copy() data_mc['weight'] = 1.0 mc = get_mc_hists_posterior(data_mc_with_weights,muon_hists,data_hists,atmo_scale_factor,muon_scale_factor,x,bins)[id] N = mc.sum() # Fix a bug in scipy(). See https://github.com/scipy/scipy/issues/8235 (I think). mc = mc + 1e-10 p = mc/mc.sum() chi2_data = nllr(data_hists[id],mc) # To draw the multinomial samples we first draw the expected number of # events from a Poisson distribution and then loop over the counts and # unique values. The reason we do this is that you can't call # multinomial.rvs with a multidimensional `n` array, and looping over every # single entry takes forever ns = np.random.poisson(N,size=1000) samples = [] for n, count in zip(*np.unique(ns,return_counts=True)): samples.append(multinomial.rvs(n,p,size=count)) samples = np.concatenate(samples) # Calculate the negative log likelihood ratio for the data simulated under # the null hypothesis
#!/usr/bin/env python
# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <https://www.gnu.org/licenses/>.

from __future__ import print_function, division
import yaml
try:
    from yaml import CLoader as Loader
except ImportError:
    from yaml.loader import SafeLoader as Loader
import string
from os.path import split, join, abspath
import uuid
from subprocess import check_call
import os
import sys

# Next two functions are a backport of the shutil.which() function from Python
# 3.3 from Lib/shutil.py in the CPython code See
# https://github.com/python/cpython/blob/master/Lib/shutil.py.

# Check that a given file can be accessed with the correct mode.
# Additionally check that `file` is not a directory, as on Windows
# directories pass the os.access check.
def _access_check(fn, mode):
    return (os.path.exists(fn) and os.access(fn, mode) and not os.path.isdir(fn))

def which(cmd, mode=os.F_OK | os.X_OK, path=None):
    """Given a command, mode, and a PATH string, return the path which
    conforms to the given mode on the PATH, or None if there is no such
    file.
    `mode` defaults to os.F_OK | os.X_OK. `path` defaults to the result
    of os.environ.get("PATH"), or can be overridden with a custom search
    path.
    """
    # If we're given a path with a directory part, look it up directly rather
    # than referring to PATH directories. This includes checking relative to the
    # current directory, e.g. ./script
    if os.path.dirname(cmd):
        if _access_check(cmd, mode):
            return cmd
        return None

    if path is None:
        path = os.environ.get("PATH", None)
        if path is None:
            try:
                path = os.confstr("CS_PATH")
            except (AttributeError, ValueError):
                # os.confstr() or CS_PATH is not available
                path = os.defpath
        # bpo-35755: Don't use os.defpath if the PATH environment variable is
        # set to an empty string

    # PATH='' doesn't match, whereas PATH=':' looks in the current directory
    if not path:
        return None

    path = path.split(os.pathsep)

    if sys.platform == "win32":
        # The current directory takes precedence on Windows.
        curdir = os.curdir
        if curdir not in path:
            path.insert(0, curdir)

        # PATHEXT is necessary to check on Windows.
        pathext = os.environ.get("PATHEXT", "").split(os.pathsep)
        # See if the given file matches any of the expected path extensions.
        # This will allow us to short circuit when given "python.exe".
        # If it does match, only test that one, otherwise we have to try
        # others.
        if any(cmd.lower().endswith(ext.lower()) for ext in pathext):
            files = [cmd]
        else:
            files = [cmd + ext for ext in pathext]
    else:
        # On other platforms you don't have things like PATHEXT to tell you
        # what file suffixes are executable, so just pass on cmd as-is.
        files = [cmd]

    seen = set()
    for dir in path:
        normdir = os.path.normcase(dir)
        if not normdir in seen:
            seen.add(normdir)
            for thefile in files:
                name = os.path.join(dir, thefile)
                if _access_check(name, mode):
                    return name
    return None

CONDOR_TEMPLATE = \
"""
# We need the job to run our executable script, with the
#  input.txt filename as an argument, and to transfer the
#  relevant input and output files:
executable = @executable
arguments = @args
transfer_input_files = @transfer_input_files

transfer_output_files = @transfer_output_files

error = @error
output = @output
log = @log

# The below are good base requirements for first testing jobs on OSG, 
#  if you don't have a good idea of memory and disk usage.
requirements = (HAS_MODULES == True) && (OSGVO_OS_STRING == "RHEL 7") && (OpSys == "LINUX")
request_cpus = 1
request_memory = 1 GB
request_disk = 1 GB

# Queue one job with the above specifications.
queue 1

+ProjectName = "SNOplus"
""".strip()

# all files required to run the fitter (except the DQXX files)
INPUT_FILES = ["muE_water_liquid.txt","pmt_response_qoca_d2o_20060216.dat","rsp_rayleigh.dat","e_water_liquid.txt","pmt_pcath_response.dat","pmt.txt","muE_deuterium_oxide_liquid.txt","pmt_response.dat","proton_water_liquid.txt"]

# generate a UUID to append to all the filenames so that if we run the same job
# twice we don't overwrite the first job
ID = uuid.uuid1()

class MyTemplate(string.Template):
    delimiter = '@'

def splitext(path):
    """
    Like os.path.splitext() except it returns the full extension if the
    filename has multiple extensions, for example:

        splitext('foo.tar.gz') -> 'foo', '.tar.gz'
    """
    full_root, full_ext = os.path.splitext(path)
    while True:
        root, ext = os.path.splitext(full_root)
        if ext:
            full_ext = ext + full_ext
            full_root = root
        else:
            break

    return full_root, full_ext

def submit_job(filename, run, gtid, dir, dqxx_dir, min_nhit, max_particles, particle_combo=None):
    print("submitting job for %s gtid %i" % (filename, gtid))
    head, tail = split(filename)
    root, ext = splitext(tail)

    # all output files are prefixed with FILENAME_GTID_UUID
    if particle_combo:
        prefix = "%s_%08i_%i_%s" % (root,gtid,particle_combo,ID.hex)
    else:
        prefix = "%s_%08i_%s" % (root,gtid,ID.hex)

    # fit output filename
    output = "%s.hdf5" % prefix
    # condor submit filename
    condor_submit = "%s.submit" % prefix

    # set up the arguments for the template
    executable = which("fit")
    wrapper = which("fit-wrapper")

    if executable is None:
        print("couldn't find fit in path!",file=sys.stderr)
        sys.exit(1)

    if wrapper is None:
        print("couldn't find fit-wrapper in path!",file=sys.stderr)
        sys.exit(1)

    args = [tail,"-o",output,"--gtid",gtid,"--min-nhit",min_nhit,"--max-particles",max_particles]
    if particle_combo:
        args += ["-p",particle_combo]
    transfer_input_files = ",".join([executable,filename,join(dqxx_dir,"DQXX_%010i.dat" % run)] + [join(dir,filename) for filename in INPUT_FILES])
    transfer_output_files = ",".join([output])

    condor_error = "%s.error" % prefix
    condor_output = "%s.output" % prefix
    condor_log = "%s.log" % prefix

    template = MyTemplate(CONDOR_TEMPLATE)

    submit_string = template.safe_substitute(
            executable=wrapper,
            args=" ".join(map(str,args)),
            transfer_input_files=transfer_input_files,
            transfer_output_files=transfer_output_files,
            error=condor_error,
            output=condor_output,
            log=condor_log)
    
    # write out the formatted template
    with open(condor_submit, "w") as f:
        f.write(submit_string)

    # submit the job
    check_call(["condor_submit",condor_submit])

def array_to_particle_combo(combo):
    particle_combo = 0
    for i, id in enumerate(combo[::-1]):
        particle_combo += id*100**i
    return particle_combo

if __name__ == '__main__':
    import argparse
    from subprocess import check_call
    import os
    import tempfile
    import h5py
    from itertools import combinations_with_replacement

    parser = argparse.ArgumentParser("submit grid jobs")
    parser.add_argument("filenames", nargs='+', help="input files")
    parser.add_argument("--min-nhit", type=int, help="minimum nhit to fit an event", default=100)
    parser.add_argument("--max-particles", type=int, help="maximum number of particles to fit for", default=3)
    parser.add_argument("--skip-second-event", action='store_true', help="only fit the first event after a MAST bank", default=False)
    parser.add_argument("--state", type=str, help="state file", default=None)
    args = parser.parse_args()

    if args.state is None:
        home = os.path.expanduser("~")
        args.state = join(home,'state.txt')

    if 'SDDM_DATA' not in os.environ:
        print("Please set the SDDM_DATA environment variable to point to the fitter source code location", file=sys.stderr)
        sys.exit(1)

    dir = os.environ['SDDM_DATA']

    if 'DQXX_DIR' not in os.environ:
        print("Please set the DQXX_DIR environment variable to point to the directory with the DQXX files", file=sys.stderr)
        sys.exit(1)

    dqxx_dir = os.environ['DQXX_DIR']

    args.state = abspath(args.state)

    # get the current working directory
    home_dir = os.getcwd()

    # get absolute paths since we are going to create a new directory
    dir = abspath(dir)
    dqxx_dir = abspath(dqxx_dir)

    zdab_cat = which("zdab-cat")

    if zdab_cat is None:
        print("couldn't find zdab-cat in path!",file=sys.stderr)
        sys.exit(1)

    for filename in args.filenames:
        filename = abspath(filename)

        head, tail = split(filename)
        root, ext = splitext(tail)

        try:
            with open(args.state) as f:
                state = yaml.load(f.read(),Loader=Loader)
        except IOError:
            state = None

        if state is None:
            state = []

        if tail in state:
            print("skipping %s" % filename)
            continue

        with open(os.devnull, 'w') as f:
            # Create a temporary file to store the events. Docs recommended
            # this method instead of mkstemp(), but I think they're the same.
            output = tempfile.NamedTemporaryFile(suffix='.hdf5',delete=False)
            output.close()

            if args.skip_second_event:
                check_call([zdab_cat,"--skip-second-event",filename,"-o",output.name],stderr=f)
            else:
                check_call([zdab_cat,filename,"-o",output.name],stderr=f)

        new_dir = "%s_%s" % (root,ID.hex)

        os.mkdir(new_dir)
        os.chdir(new_dir)

        with h5py.File(output.name) as f:
            for ev in f['ev']:
                if ev['nhit'] >= args.min_nhit:
                    for i in range(1,args.max_particles+1):
                        for particle_combo in map(array_to_particle_combo,combinations_with_replacement([20,22],i)):
                            submit_job(filename, ev['run'], ev['gtid'], dir, dqxx_dir, args.min_nhit, args.max_particles, particle_combo)

        # Delete temporary HDF5 file
        os.unlink(output.name)

        state.append(tail)

        with open(args.state,"w") as f:
            f.write(yaml.dump(state,default_flow_style=False))

        os.chdir(home_dir)