import numpy as np from matplotlib import pyplot as plt #import math #from random import * #import os, sys #from scipy import stats excess_xerr = 0.0666667/2 offset = np.linspace(-1.0+excess_xerr,1.0-excess_xerr, 30) angle_psrj1256 = 0.374 angle_psrj1334 = -0.264 angle_snr = 0.2602 # ALL ENERGIES excess_all = [521.1255, 548.7674, 580.0453, 597.8476, 787.0257, 923.8861, 878.4341, 1306.128, 1407.358, 1223.267, 1172.873, 740.7292, 439.8151, 328.8221, 366.598, 337.386, 364.2495, 363.6686, 325.6028, 414.9057, 275.7392, 175.2247, 174.4922, 142.5355, 145.9701, 57.54354, 1.308359, 42.69603, 95.1851, 76.44855] excess_all_error = [37.72767, 38.74391, 38.95618, 37.32576, 41.25492, 42.83484, 40.74205, 47.0988, 48.23293, 44.93249, 45.92638, 41.00099, 35.35579, 35.77574, 36.56149, 34.43562, 37.04926, 35.33955, 36.81458, 38.54799, 34.73767, 35.54197, 35.95124, 33.57102, 36.3788, 35.3527, 32.24815, 36.43485, 38.84224, 36.92802] # E > 5 TeV excess_5TeV = [22.42065, 27.71609, 30.87044, 41.61195, 53.70641, 62.37878, 68.86499, 108.5662, 135.9365, 175.4837, 153.3364, 99.02442, 58.34442, 39.57725, 40.21133, 47.88735, 48.11093, 62.53881, 81.11655, 92.22791, 70.2512, 37.73834, 13.89038, 20.28194, 14.43172, 7.151055, 1.500369, 5.044003, -0.4935792, 9.924324] excess_5TeV_error = [9.303891, 9.783146, 9.886861, 9.953833, 11.02839, 11.48784, 11.38929, 13.43432, 14.47171, 15.5168, 15.15442, 13.21378, 11.09352, 10.7116, 10.78714, 10.6974, 11.22607, 11.46398, 12.73333, 13.25088, 11.91597, 10.8714, 9.63815, 9.489243, 9.784595, 9.414605, 8.49435, 9.62996, 9.780595, 10.12466] # E > 10 TeV excess_10TeV = [9.739837, 19.75674, 13.54574, 4.442139, 17.29932, 19.41398, 25.92578, 43.8474, 46.4514, 78.29337, 70.23325, 38.30177, 26.52231, 25.2608, 14.53392, 8.518459, 25.94258, 30.08483, 32.38992, 42.80378, 36.67501, 25.40776, 7.471138, 5.622364, 9.07112, 6.886342, 0.07748762, 4.23087, 0.7578744, 7.978853] excess_10TeV_error = [6.820231, 7.651236, 7.173386, 6.055551, 7.407325, 7.553343, 7.672439, 9.109904, 9.308168, 10.71244, 10.64418, 9.006179, 7.981633, 8.324583, 7.66388, 6.822831, 8.498838, 8.423692, 8.918771, 9.590113, 8.904451, 8.578166, 7.332127, 6.734552, 7.473491, 7.305629, 6.280602, 7.32166, 7.485886, 7.832453] # E > 15 TeV excess_15TeV = [5.312231, 12.82182, 2.410167, 3.304955, 3.059156, 9.226726, 12.57422, 18.38869, 19.57457, 44.38956, 35.06401, 14.58524, 14.0258, 15.88909, 6.432868, 2.493559, 13.93614, 10.0936, 15.39778, 15.50967, 22.42336, 20.21723, 7.165771, 4.147635, 4.016008, 2.851274, 1.882082, 4.511755, -1.648037, 5.080554] excess_15TeV_error = [5.74442, 6.483953, 5.589754, 5.306178, 5.65918, 6.187753, 6.15414, 6.940955, 7.056964, 8.439388, 8.16096, 6.777316, 6.423609, 6.931186, 6.211331, 5.488662, 6.850272, 6.200501, 6.968262, 7.017451, 7.242217, 7.434093, 6.434485, 5.771478, 6.134125, 5.989775, 5.48604, 6.335043, 6.269552, 6.718818] # E > 20 TeV excess_20TeV = [1.837435, 5.516666, 0.3989009, 3.016002, 1.19327, 4.385159, 7.311587, 13.65967, 12.77408, 27.7379, 17.88296, 10.84613, 7.6423, 14.18009, 7.025218, 2.69179, 8.344194, 5.999498, 9.633989, 9.611696, 13.61534, 12.84163, 3.100312, 0.6406552, 4.077078, 1.874961, -0.5679542, 2.188177, -0.8166772, 3.876028] excess_20TeV_error = [5.139643, 5.615793, 5.09901, 4.998228, 5.176802, 5.472493, 5.415951, 6.277728, 6.220985, 7.099044, 6.671078, 6.138275, 5.558789, 6.468019, 5.891267, 5.140079, 6.035342, 5.491554, 6.145774, 6.167697, 6.208123, 6.465871, 5.60321, 4.959664, 5.644316, 5.394306, 4.754541, 5.608642, 5.796643, 6.096053] plt.figure(figsize=(10, 12), dpi=50, facecolor='w', edgecolor='k') plt.rcParams["font.weight"] = "bold" plt.subplot(5, 1, 1) plt.axvline(x=angle_psrj1256, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, label='PSR J1826$-$1256', color='red') plt.axvline(x=angle_psrj1334, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, label='PSR J1826$-$1334', color='blue') plt.axvline(x=angle_snr, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, label='SNR G18.45$-$0.42', color='green') plt.errorbar(offset, excess_all, xerr=excess_xerr, yerr=excess_all_error, label='Data (All energies)', fmt='o', ecolor='black', color='black') plt.legend(loc='upper left',fontsize=11) plt.ylabel(r'Excess', fontsize=12, fontweight='bold') plt.tick_params(axis='y', labelsize=12) plt.tick_params(axis='x', labelsize=12) plt.legend(loc='upper right', fontsize=10) plt.subplot(5, 1, 2) plt.axvline(x=angle_psrj1256, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='red') plt.axvline(x=angle_psrj1334, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='blue') plt.axvline(x=angle_snr, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='green') plt.errorbar(offset, excess_5TeV, xerr=excess_xerr, yerr=excess_5TeV_error, label='Data (E > 5 TeV)', fmt='o', ecolor='black', color='black') plt.legend(loc='upper left',fontsize=11) plt.ylabel(r'Excess', fontsize=12, fontweight='bold') plt.tick_params(axis='y', labelsize=12) plt.tick_params(axis='x', labelsize=12) plt.legend(loc='upper right', fontsize=10) plt.subplot(5, 1, 3) plt.axvline(x=angle_psrj1256, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='red') plt.axvline(x=angle_psrj1334, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='blue') plt.axvline(x=angle_snr, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='green') plt.errorbar(offset, excess_10TeV, xerr=excess_xerr, yerr=excess_10TeV_error, label='Data (E > 10 TeV)', fmt='o', ecolor='black', color='black') plt.legend(loc='upper left',fontsize=11) plt.ylabel(r'Excess', fontsize=12, fontweight='bold') plt.tick_params(axis='y', labelsize=12) plt.tick_params(axis='x', labelsize=12) plt.legend(loc='upper right', fontsize=10) plt.subplot(5, 1, 4) plt.axvline(x=angle_psrj1256, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='red') plt.axvline(x=angle_psrj1334, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='blue') plt.axvline(x=angle_snr, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='green') plt.errorbar(offset, excess_15TeV, xerr=excess_xerr, yerr=excess_15TeV_error, label='Data (E > 15 TeV)', fmt='o', ecolor='black', color='black') plt.legend(loc='upper left',fontsize=11) plt.ylabel(r'Excess', fontsize=12, fontweight='bold') plt.tick_params(axis='y', labelsize=12) plt.tick_params(axis='x', labelsize=12) plt.legend(loc='upper right', fontsize=10) plt.subplot(5, 1, 5) plt.axvline(x=angle_psrj1256, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='red') plt.axvline(x=angle_psrj1334, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='blue') plt.axvline(x=angle_snr, ymin=0, ymax=max(excess_all), ls='--', lw=2.0, color='green') plt.errorbar(offset, excess_20TeV, xerr=excess_xerr, yerr=excess_20TeV_error, label='Data (E > 20 TeV)', fmt='o', ecolor='black', color='black') plt.legend(loc='upper left',fontsize=11) plt.ylabel(r'Excess', fontsize=12, fontweight='bold') plt.xlabel(r'Projected angle ($^\circ$)', fontsize=12, fontweight='bold') plt.tick_params(axis='y', labelsize=12) plt.tick_params(axis='x', labelsize=12) plt.legend(loc='upper right', fontsize=10) #Save file plt.savefig('excess_slices.pdf', dpi=300)