#!/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 # 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") IDP_E_MINUS = 20 IDP_MU_MINUS = 22 SNOMAN_MASS = { 20: 0.511, 21: 0.511, 22: 105.658, 23: 105.658 } def plot_hist(x, label=None): # determine the bin width using the Freedman Diaconis rule # see https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule h = 2*iqr(x)/len(x)**(1/3) n = max(int((np.max(x)-np.min(x))/h),10) bins = np.linspace(np.min(x),np.max(x),n) plt.hist(x, bins=bins, histtype='step', label=label) def plot_legend(n): plt.figure(n) ax = plt.gca() handles, labels = ax.get_legend_handles_labels() new_handles = [Line2D([],[],c=h.get_edgecolor()) for h in handles] plt.legend(handles=new_handles,labels=labels) def get_stats(x): """ Returns a tuple (mean, error mean, std, error std) for the values in x. The formula for the standard error on the standard deviation comes from https://stats.stackexchange.com/questions/156518. """ mean = np.mean(x) std = np.std(x) n = len(x) u4 = np.mean((x-mean)**4) error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std) return mean, std/np.sqrt(n), std, error if __name__ == '__main__': import argparse import matplotlib.pyplot as plt import numpy as np import h5py import pandas as pd parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") args = parser.parse_args() for filename in args.filenames: print(filename) with h5py.File(filename) as f: ev = pd.read_hdf(filename, "ev") mcgn = pd.read_hdf(filename, "mcgn") fits = pd.read_hdf(filename, "fits") # get rid of 2nd events like Michel electrons ev = ev.sort_values(['run','gtid']).groupby(['evn'],as_index=False).nth(0) # Now, we merge all three datasets together to produce a single # dataframe. To do so, we join the ev dataframe with the mcgn frame # on the evn column, and then join with the fits on the run and # gtid columns. # # At the end we will have a single dataframe with one row for each # fit, i.e. it will look like: #

/* 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/>.
 */

#ifndef QUANTUM_EFFICIENCY_H
#define QUANTUM_EFFICIENCY_H

double get_quantum_efficiency(double wavelength);

#endif