diff options
Diffstat (limited to 'utils/plot-energy')
| -rwxr-xr-x | utils/plot-energy | 420 |
1 files changed, 363 insertions, 57 deletions
diff --git a/utils/plot-energy b/utils/plot-energy index 8539040..5079c15 100755 --- a/utils/plot-energy +++ b/utils/plot-energy @@ -13,12 +13,29 @@ # # You should have received a copy of the GNU General Public License along with # this program. If not, see <https://www.gnu.org/licenses/>. - +""" +Script to plot final fit results along with sidebands for the dark matter +analysis. To run it just run: + + $ ./plot-energy [list of fit results] + +Currently it will plot energy distributions for external muons, michel +electrons, atmospheric events with neutron followers, and prompt signal like +events. Each of these plots will have a different subplot for the particle ID +of the best fit, i.e. single electron, single muon, double electron, electron + +muon, or double muon. + +When run with the --dc command line argument it instead produces corner plots +showing the distribution of the high level variables used in the contamination +analysis for all the different instrumental backgrounds and external muons. +""" 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 +from scipy.special import spence +from itertools import izip_longest PSUP_RADIUS = 840.0 @@ -67,7 +84,39 @@ DC_FTS = 0x200 DC_ITC = 0x400 DC_BREAKDOWN = 0x800 -def plot_hist(df, title=None, muons=False): +particle_id = {20: 'e', 22: r'\mu'} + +def grouper(iterable, n, fillvalue=None): + "Collect data into fixed-length chunks or blocks" + # grouper('ABCDEFG', 3, 'x') --> ABC DEF Gxx + args = [iter(iterable)] * n + return izip_longest(fillvalue=fillvalue, *args) + +def plot_hist2(df, muons=False): + for id, df_id in sorted(df.groupby('id')): + 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) + + if muons: + plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') + plt.xlabel("log10(Energy (GeV))") + else: + plt.hist(df_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step') + plt.xlabel("Energy (MeV)") + plt.title('$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(id),2)]) + '$') + + if len(df): + plt.tight_layout() + +def plot_hist(df, muons=False): for id, df_id in sorted(df.groupby('id')): if id == 20: plt.subplot(3,4,1) @@ -92,13 +141,10 @@ def plot_hist(df, title=None, muons=False): plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') plt.xlabel("log10(Energy (GeV))") else: - plt.hist(df_id.ke.values, bins=np.linspace(20,1000e3,100), histtype='step') + 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() @@ -234,8 +280,8 @@ def atmospheric_events(ev): # neutron followers have to obey stricter set of data cleaning cuts neutron = follower[follower.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ESUM | DC_OWL | DC_OWL_TRIGGER | DC_FTS) == 0] neutron = neutron[~np.isnan(neutron.ftp_x) & ~np.isnan(neutron.rsp_energy)] - r = np.sqrt(neutron.ftp_x**2 + neutron.ftp_y**2 + neutron.ftp_z**2) - neutron = neutron[r < AV_RADIUS] + # FIXME: What should the radius cut be here? AV? (r/r_psup)^3 < 0.9? + neutron = neutron[neutron.ftp_r < AV_RADIUS] neutron = neutron[neutron.rsp_energy > 4.0] # neutron events accepted after 20 microseconds and before 250 ms (50 ms during salt) @@ -279,6 +325,108 @@ def prompt_event(ev): ev.loc[ev.prompt,'prompt'] &= np.concatenate(([True],np.diff(ev[ev.prompt].gtr.values) > 250e6)) return ev +# Taken from https://raw.githubusercontent.com/mwaskom/seaborn/c73055b2a9d9830c6fbbace07127c370389d04dd/seaborn/utils.py +def despine(fig=None, ax=None, top=True, right=True, left=False, + bottom=False, offset=None, trim=False): + """Remove the top and right spines from plot(s). + + fig : matplotlib figure, optional + Figure to despine all axes of, default uses current figure. + ax : matplotlib axes, optional + Specific axes object to despine. + top, right, left, bottom : boolean, optional + If True, remove that spine. + offset : int or dict, optional + Absolute distance, in points, spines should be moved away + from the axes (negative values move spines inward). A single value + applies to all spines; a dict can be used to set offset values per + side. + trim : bool, optional + If True, limit spines to the smallest and largest major tick + on each non-despined axis. + + Returns + ------- + None + + """ + # Get references to the axes we want + if fig is None and ax is None: + axes = plt.gcf().axes + elif fig is not None: + axes = fig.axes + elif ax is not None: + axes = [ax] + + for ax_i in axes: + for side in ["top", "right", "left", "bottom"]: + # Toggle the spine objects + is_visible = not locals()[side] + ax_i.spines[side].set_visible(is_visible) + if offset is not None and is_visible: + try: + val = offset.get(side, 0) + except AttributeError: + val = offset + _set_spine_position(ax_i.spines[side], ('outward', val)) + + # Potentially move the ticks + if left and not right: + maj_on = any( + t.tick1line.get_visible() + for t in ax_i.yaxis.majorTicks + ) + min_on = any( + t.tick1line.get_visible() + for t in ax_i.yaxis.minorTicks + ) + ax_i.yaxis.set_ticks_position("right") + for t in ax_i.yaxis.majorTicks: + t.tick2line.set_visible(maj_on) + for t in ax_i.yaxis.minorTicks: + t.tick2line.set_visible(min_on) + + if bottom and not top: + maj_on = any( + t.tick1line.get_visible() + for t in ax_i.xaxis.majorTicks + ) + min_on = any( + t.tick1line.get_visible() + for t in ax_i.xaxis.minorTicks + ) + ax_i.xaxis.set_ticks_position("top") + for t in ax_i.xaxis.majorTicks: + t.tick2line.set_visible(maj_on) + for t in ax_i.xaxis.minorTicks: + t.tick2line.set_visible(min_on) + + if trim: + # clip off the parts of the spines that extend past major ticks + xticks = ax_i.get_xticks() + if xticks.size: + firsttick = np.compress(xticks >= min(ax_i.get_xlim()), + xticks)[0] + lasttick = np.compress(xticks <= max(ax_i.get_xlim()), + xticks)[-1] + ax_i.spines['bottom'].set_bounds(firsttick, lasttick) + ax_i.spines['top'].set_bounds(firsttick, lasttick) + newticks = xticks.compress(xticks <= lasttick) + newticks = newticks.compress(newticks >= firsttick) + ax_i.set_xticks(newticks) + + yticks = ax_i.get_yticks() + if yticks.size: + firsttick = np.compress(yticks >= min(ax_i.get_ylim()), + yticks)[0] + lasttick = np.compress(yticks <= max(ax_i.get_ylim()), + yticks)[-1] + ax_i.spines['left'].set_bounds(firsttick, lasttick) + ax_i.spines['right'].set_bounds(firsttick, lasttick) + newticks = yticks.compress(yticks <= lasttick) + newticks = newticks.compress(newticks >= firsttick) + ax_i.set_yticks(newticks) + def plot_corner_plot(ev, title, save=None): variables = ['r_psup','psi','z','udotr'] labels = [r'$(r/r_\mathrm{PSUP})^3$',r'$\psi$','z',r'$\vec{u}\cdot\vec{r}$'] @@ -287,7 +435,8 @@ def plot_corner_plot(ev, title, save=None): ev = ev.dropna(subset=variables) - f = plt.figure(figsize=(8,8)) + f = plt.figure(figsize=(FIGSIZE[0],FIGSIZE[0])) + despine(fig,trim=True) for i in range(len(variables)): for j in range(len(variables)): if j > i: @@ -328,8 +477,11 @@ def plot_corner_plot(ev, title, save=None): if i != j: plt.axhline(cuts[i],color='k',ls='--',alpha=0.5) + plt.tight_layout() + if save: - plt.savefig(save) + plt.savefig(save + ".pdf") + plt.savefig(save + ".eps") plt.suptitle(title) @@ -501,6 +653,37 @@ def std_err(x): error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std) return error +# Fermi constant +GF = 1.16637887e-5 # 1/MeV^2 +ELECTRON_MASS = 0.5109989461 # MeV +MUON_MASS = 105.6583745 # MeV +PROTON_MASS = 938.272081 # MeV +FINE_STRUCTURE_CONSTANT = 7.297352566417e-3 + +def f(x): + y = (5/(3*x**2) + 16*x/3 + 4/x + (12-8*x)*np.log(1/x-1) - 8)*np.log(MUON_MASS/ELECTRON_MASS) + y += (6-4*x)*(2*spence(x) - 2*np.log(x)**2 + np.log(x) + np.log(1-x)*(3*np.log(x)-1/x-1) - np.pi**2/3-2) + y += (1-x)*(34*x**2+(5-34*x**2+17*x)*np.log(x) - 22*x)/(3*x**2) + y += 6*(1-x)*np.log(x) + return y + +def michel_spectrum(T): + """ + Michel electron energy spectrum for a free muon. `T` should be the kinetic + energy of the electron or positron in MeV. + + Note: The result is not normalized. + + From https://arxiv.org/abs/1406.3575. + """ + E = T + ELECTRON_MASS + x = 2*E/MUON_MASS + mask = (x > 0) & (x < 1) + y = np.zeros_like(x,dtype=np.double) + y[mask] = GF**2*MUON_MASS**5*x[mask]**2*(6-4*x[mask]+FINE_STRUCTURE_CONSTANT*f(x[mask])/np.pi)/(192*np.pi**3) + y *= 2*MUON_MASS + return y + if __name__ == '__main__': import argparse import matplotlib.pyplot as plt @@ -512,6 +695,7 @@ if __name__ == '__main__': parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") parser.add_argument("--dc", action='store_true', default=False, help="plot corner plots for backgrounds") + parser.add_argument("--save", action='store_true', default=False, help="save corner plots for backgrounds") args = parser.parse_args() ev = pd.concat([pd.read_hdf(filename, "ev") for filename in args.filenames],ignore_index=True) @@ -618,6 +802,9 @@ if __name__ == '__main__': fits['r'] = np.sqrt(fits.x**2 + fits.y**2 + fits.z**2) fits['r_psup'] = (fits['r']/PSUP_RADIUS)**3 + ev['ftp_r'] = np.sqrt(ev.ftp_x**2 + ev.ftp_y**2 + ev.ftp_z**2) + ev['ftp_r_psup'] = (ev['ftp_r']/PSUP_RADIUS)**3 + print("number of events = %i" % len(ev)) # Now, select prompt events. @@ -646,6 +833,47 @@ if __name__ == '__main__': # retrigger cut ev = ev.groupby('run',group_keys=False).apply(retrigger_cut) + if args.save: + # default \textwidth for a fullpage article in Latex is 16.50764 cm. + # You can figure this out by compiling the following TeX document: + # + # \documentclass{article} + # \usepackage{fullpage} + # \usepackage{layouts} + # \begin{document} + # textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth} + # \end{document} + + width = 16.50764 + width /= 2.54 # cm -> inches + # According to this page: + # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf, + # Tufte suggests an aspect ratio of 1.5 - 1.6. + height = width/1.5 + FIGSIZE = (width,height) + + import matplotlib.pyplot as plt + + font = {'family':'serif', 'serif': ['computer modern roman']} + plt.rc('font',**font) + + plt.rc('text', usetex=True) + else: + # 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") + + import matplotlib.pyplot as plt + + # Default figure size. Currently set to my monitor width and height so that + # things are properly formatted + FIGSIZE = (13.78,7.48) + + # Make the defalt font bigger + plt.rc('font', size=22) + if args.dc: ev = ev[ev.prompt] @@ -684,17 +912,29 @@ if __name__ == '__main__': print(signal) # save as PDF b/c EPS doesn't support alpha values - plot_corner_plot(breakdown,"Breakdowns",save="breakdown_corner_plot.pdf") - plot_corner_plot(muon,"Muons",save='muon_corner_plot.pdf') - plot_corner_plot(flashers,"Flashers",save="flashers_corner_plot.pdf") - plot_corner_plot(neck,"Neck",save="neck_corner_plot.pdf") - plot_corner_plot(noise,"Noise",save="noise_corner_plot.pdf") - plot_corner_plot(signal,"Signal",save="signal_corner_plot.pdf") - - plt.figure() - plot_hist(flashers,"Flashers") - plt.figure() - plot_hist(muon,"Muons",muons=True) + if args.save: + plot_corner_plot(breakdown,"Breakdowns",save="breakdown_corner_plot") + plot_corner_plot(muon,"Muons",save="muon_corner_plot") + plot_corner_plot(flashers,"Flashers",save="flashers_corner_plot") + plot_corner_plot(neck,"Neck",save="neck_corner_plot") + plot_corner_plot(noise,"Noise",save="noise_corner_plot") + plot_corner_plot(signal,"Signal",save="signal_corner_plot") + else: + plot_corner_plot(breakdown,"Breakdowns") + plot_corner_plot(muon,"Muons") + plot_corner_plot(flashers,"Flashers") + plot_corner_plot(neck,"Neck") + plot_corner_plot(noise,"Noise") + plot_corner_plot(signal,"Signal") + + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(flashers) + despine(fig,trim=True) + plt.suptitle("Flashers") + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(muon,muons=True) + despine(fig,trim=True) + plt.suptitle("Muons") plt.show() sys.exit(0) @@ -816,25 +1056,57 @@ if __name__ == '__main__': # Note: Need to design and apply a psi based cut here - plt.figure() - plot_hist(prompt,"Without Neutron Follower") - plt.figure() - plot_hist(atm,"With Neutron Follower") - plt.figure() - plot_hist(michel_best_fit,"Michel Electrons") - plt.figure() - plot_hist(muon_best_fit,"External Muons") + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(prompt) + despine(fig,trim=True) + if args.save: + plt.savefig("prompt.pdf") + plt.savefig("prompt.eps") + else: + plt.suptitle("Without Neutron Follower") + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(atm) + despine(fig,trim=True) + if args.save: + plt.savefig("atm.pdf") + plt.savefig("atm.eps") + else: + plt.suptitle("With Neutron Follower") + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(michel_best_fit) + despine(fig,trim=True) + if args.save: + plt.savefig("michel_electrons.pdf") + plt.savefig("michel_electrons.eps") + else: + plt.suptitle("Michel Electrons") + fig = plt.figure(figsize=FIGSIZE) + plot_hist2(muon_best_fit,muons=True) + despine(fig,trim=True) + if len(muon_best_fit): + plt.tight_layout() + if args.save: + plt.savefig("external_muons.pdf") + plt.savefig("external_muons.eps") + else: + plt.suptitle("External Muons") # Plot the energy and angular distribution for external muons - plt.figure() + fig = plt.figure(figsize=FIGSIZE) plt.subplot(2,1,1) 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') + despine(fig,trim=True) plt.xlabel(r"$\cos(\theta)$") - plt.suptitle("Muons") + plt.tight_layout() + if args.save: + plt.savefig("muon_energy_cos_theta.pdf") + plt.savefig("muon_energy_cos_theta.eps") + else: + plt.suptitle("Muons") # For the Michel energy plot, we only look at the single particle electron # fit @@ -842,41 +1114,75 @@ if __name__ == '__main__': stopping_muons = pd.merge(muons,michel,left_on=['run','gtid'],right_on=['run','muon_gtid'],suffixes=('','_michel')) - # project muon to PSUP - stopping_muons['dx'] = stopping_muons.apply(get_dx,axis=1) - # energy based on distance travelled - stopping_muons['T_dx'] = dx_to_energy(stopping_muons.dx) - stopping_muons['dT'] = stopping_muons['energy1'] - stopping_muons['T_dx'] - - plt.figure() - plt.hist((stopping_muons['energy1']-stopping_muons['T_dx'])*100/stopping_muons['T_dx'], bins=np.linspace(-100,100,200), histtype='step') - plt.xlabel("Fractional energy difference (\%)") - plt.title("Fractional energy difference for Stopping Muons") + if len(stopping_muons): + # project muon to PSUP + stopping_muons['dx'] = stopping_muons.apply(get_dx,axis=1) + # energy based on distance travelled + stopping_muons['T_dx'] = dx_to_energy(stopping_muons.dx) + stopping_muons['dT'] = stopping_muons['energy1'] - stopping_muons['T_dx'] + + fig = plt.figure(figsize=FIGSIZE) + plt.hist((stopping_muons['energy1']-stopping_muons['T_dx'])*100/stopping_muons['T_dx'], bins=np.linspace(-100,100,200), histtype='step') + despine(fig,trim=True) + plt.xlabel("Fractional energy difference (\%)") + plt.title("Fractional energy difference for Stopping Muons") + plt.tight_layout() + if args.save: + plt.savefig("stopping_muon_fractional_energy_difference.pdf") + plt.savefig("stopping_muon_fractional_energy_difference.eps") + else: + plt.title("Stopping Muon Fractional Energy Difference") - # 100 bins between 50 MeV and 10 GeV - bins = np.arange(50,10000,1000) + # 100 bins between 50 MeV and 10 GeV + bins = np.arange(50,10000,1000) - pd_bins = pd.cut(stopping_muons['energy1'],bins) + pd_bins = pd.cut(stopping_muons['energy1'],bins) - T = (bins[1:] + bins[:-1])/2 + T = (bins[1:] + bins[:-1])/2 - dT = stopping_muons.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err]) + dT = stopping_muons.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err]) - plt.figure() - plt.errorbar(T,dT['median']*100/T,yerr=dT['median_err']*100/T) - plt.xlabel("Kinetic Energy (MeV)") - plt.ylabel(r"Energy bias (\%)") + fig = plt.figure(figsize=FIGSIZE) + plt.errorbar(T,dT['median']*100/T,yerr=dT['median_err']*100/T) + despine(fig,trim=True) + plt.xlabel("Kinetic Energy (MeV)") + plt.ylabel(r"Energy bias (\%)") + plt.tight_layout() + if args.save: + plt.savefig("stopping_muon_energy_bias.pdf") + plt.savefig("stopping_muon_energy_bias.eps") + else: + plt.title("Stopping Muon Energy Bias") - plt.figure() - plt.errorbar(T,dT['iqr_std']*100/T,yerr=dT['iqr_std_err']*100/T) - plt.xlabel("Kinetic Energy (MeV)") - plt.ylabel(r"Energy resolution (\%)") + fig = plt.figure(figsize=FIGSIZE) + plt.errorbar(T,dT['iqr_std']*100/T,yerr=dT['iqr_std_err']*100/T) + despine(fig,trim=True) + plt.xlabel("Kinetic Energy (MeV)") + plt.ylabel(r"Energy resolution (\%)") + plt.tight_layout() + if args.save: + plt.savefig("stopping_muon_energy_resolution.pdf") + plt.savefig("stopping_muon_energy_resolution.eps") + else: + plt.title("Stopping Muon Energy Resolution") - plt.figure() - plt.hist(michel.ke.values, bins=np.linspace(0,100,100), histtype='step', label="My fitter") + fig = plt.figure(figsize=FIGSIZE) + bins=np.linspace(0,100,100) + plt.hist(michel.ke.values, bins=bins, histtype='step', label="Dark Matter 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") + x = np.linspace(0,100,1000) + y = michel_spectrum(x) + y /= np.trapz(y,x=x) + N = len(michel) + plt.plot(x, N*y*(bins[1]-bins[0]), ls='--', color='k', label="Michel Spectrum") + despine(fig,trim=True) plt.xlabel("Energy (MeV)") - plt.title("Michel Electrons") + plt.tight_layout() plt.legend() + if args.save: + plt.savefig("michel_electrons_ke.pdf") + plt.savefig("michel_electrons_ke.eps") + else: + plt.title("Michel Electrons") plt.show() |