#!/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 . 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 numpy as np from scipy.stats import iqr from matplotlib.lines import Line2D # from https://stackoverflow.com/questions/287871/how-to-print-colored-text-in-terminal-in-python class bcolors: HEADER = '\033[95m' OKBLUE = '\033[94m' OKGREEN = '\033[92m' WARNING = '\033[93m' FAIL = '\033[91m' ENDC = '\033[0m' BOLD = '\033[1m' UNDERLINE = '\033[4m' # on retina screens, the default plots are way too small # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 # Qt5 will scale everything using the dpi in ~/.Xresources import matplotlib matplotlib.use("Qt5Agg") SNOMAN_MASS = { 20: 0.511, 21: 0.511, 22: 105.658, 23: 105.658 } AV_RADIUS = 600.0 # Data cleaning bitmasks. DC_MUON = 0x1 DC_JUNK = 0x2 DC_CRATE_ISOTROPY = 0x4 DC_QVNHIT = 0x8 DC_NECK = 0x10 DC_FLASHER = 0x20 DC_ESUM = 0x40 DC_OWL = 0x80 DC_OWL_TRIGGER = 0x100 DC_FTS = 0x200 def plot_hist(df, title=None): for id, df_id in sorted(df.groupby('id')): if id == 20: plt.subplot(3,4,1) elif id == 22: plt.subplot(3,4,2) elif id == 2020: plt.subplot(3,4,5) elif id == 2022: plt.subplot(3,4,6) elif id == 2222: plt.subplot(3,4,7) elif id == 202020: plt.subplot(3,4,9) elif id == 202022: plt.subplot(3,4,10) elif id == 202222: plt.subplot(3,4,11) elif id == 222222: plt.subplot(3,4,12) plt.hist(df_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step') plt.xlabel("Energy (MeV)") plt.title(str(id)) if title: plt.suptitle(title) if len(df): plt.tight_layout() def chunks(l, n): """Yield successive n-sized chunks from l.""" for i in range(0, len(l), n): yield l[i:i + n] def print_warning(msg): print(bcolors.FAIL + msg + bcolors.ENDC,file=sys.stderr) if __name__ == '__main__': import argparse import matplotlib.pyplot as plt import numpy as np import pandas as pd import sys import h5py parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") args = parser.parse_args() ev = pd.concat([pd.read_hdf(filename, "ev") for filename in args.filenames]) fits = pd.concat([pd.read_hdf(filename, "fits") for filename in args.filenames]) # This is a bit of a hack. It appears that many times the fit will # actually do much better by including a very low energy electron or # muon. I believe the reason for this is that of course my likelihood # function is not perfect (for example, I don't include the correct # angular distribution for Rayleigh scattered light), and so the fitter # often wants to add a very low energy electron or muon to fix things. # # Ideally I would fix the likelihood function, but for now we just # discard any fit results which have a very low energy electron or # muon. # # FIXME: Test this since query() is new to pandas fits = fits.query('not (n > 1 and ((id1 == 20 and energy1 < 20) or (id2 == 20 and energy2 < 20) or (id3 == 20 and energy3 < 20)))') fits

#!/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 numpy as np

# on retina screens, the default plots are way too small
# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
# Qt5 will scale everything using the dpi in ~/.Xresources
import matplotlib
matplotlib.use("Qt5Agg")

if __name__ == '__main__':
    import argparse
    from mpl_toolkits.mplot3d import axes3d
    import matplotlib.pyplot as plt

    parser = argparse.ArgumentParser("plot likelihood function")
    parser.add_argument("filenames", nargs='+', help="input files")
    args = parser.parse_args()

    for filename in args.filenames:
        print(filename)
        data = np.genfromtxt(filename)

        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')

        X = data[:,0].reshape((50,50))
        Y = data[:,1].reshape((50,50))
        Z = data[:,2].reshape((50,50))

        ax.plot_wireframe(X, Y, Z)

    plt.show()

plt.hist(muons.ke.values, bins=np.logspace(3,7,100), histtype='step') plt.xlabel("Energy (MeV)") plt.gca().set_xscale("log") plt.subplot(2,1,2) plt.hist(np.cos(muons.theta.values), bins=np.linspace(-1,1,100), histtype='step') plt.xlabel(r"$\cos(\theta)$") plt.suptitle("Muons") # For the Michel energy plot, we only look at the single particle electron # fit michel = michel[michel.id == 20] plt.figure() plt.hist(michel.ke.values, bins=np.linspace(20,100,100), histtype='step', label="My fitter") if michel.size: plt.hist(michel[~np.isnan(michel.rsp_energy.values)].rsp_energy.values, bins=np.linspace(20,100,100), histtype='step',label="RSP") plt.xlabel("Energy (MeV)") plt.title("Michel Electrons") plt.legend() plt.show()