m_F19 = 18.9984031621
dm_F19 = 0.0000000009
m_Ca42 = 41.95861778
dm_Ca42 = 0.00000016
m_e = 0.0005485799

m_CaF = m_Ca42 + m_F19
m_CaF_anion = m_Ca42 + m_F19 - m_e

m_Rb85 =       84.9117897360
m_Rb85_anion = 84.9117897360 - m_e
dm_Rb85 =       0.0000000054 

amu_keV = 931494.102
ME_NSCL2022 = -56355.5 # keV
dME_NSCL2022 = 5 # keV
m_NSCL2022 = 61 + ME_NSCL2022 / amu_keV
dm_NSCL2022 = dME_NSCL2022 / amu_keV

ToF_ICR_stat_err = 1.1691666854510232e-06
ToF_ICR_syst_err = 60.9395042585674 * 8e-9
ToF_ICR_err = ( ToF_ICR_stat_err**2 + ToF_ICR_syst_err**2 )**0.5

PI_ICR_stat_err = 6.71562072700104e-07
PI_ICR_err = ( PI_ICR_stat_err**2 + ToF_ICR_syst_err**2 )**0.5


def m_Zn61(m_ratio):
    return (m_CaF-m_e)*m_ratio + m_e
def dm_Zn61(m_ratio, dm_ratio):
    err_stat = (m_CaF-m_e)*dm_ratio
    err_res_syst = (m_CaF-m_e) * m_ratio * 8e-9
    err_mass_shift = (m_CaF-m_CaF*m_ratio) * 1.6e-10
    err_ref = m_ratio * ( dm_Ca42**2 + dm_F19**2 )**0.5
    return ( err_stat**2 + err_ref**2 + err_res_syst**2 + err_mass_shift**2 )**0.5
def ME_Zn61(m):
    return (m-61) * amu_keV
def dME_Zn61(dm):
    return dm * amu_keV

def m_Zn61_Rb(m_ratio):
    return (m_Rb85-m_e)*m_ratio + m_e
def dm_Zn61_Rb(m_ratio, dm_ratio):
    r_err_stat =       dm_ratio
    r_err_res_syst =   m_ratio * 8e-9
    r_err_mass_shift = m_ratio * (m_Rb85-m_Rb85*m_ratio) * 1.6e-10
    # r_err_syst = ( r_err_res_syst**2 + r_err_mass_shift**2 )**0.5 # quadratic addition
    r_err_syst = abs(r_err_res_syst) + abs(r_err_mass_shift) # linear addition
    m_err_ref = m_ratio * dm_Rb85
    return ( (m_Rb85-m_e)**2 * (r_err_stat**2 + r_err_syst**2) 
            + m_err_ref**2 )**0.5

###############################################################################

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 8})
plt.rcParams.update({'pdf.fonttype': 42})
plt.close('all')

cm = 1 / 2.54
figsize = (8.5*cm, 5*cm)
fig, ax = plt.subplots(figsize=figsize)

plt.ylim(-56349-31, -56349+21)

color_salm013 = (0.99745098, 0.9352549 , 0.92811765)
plt.fill_between((0.5, 4.5), -56349-16, -56349+16, 
                 color=color_salm013, lw=0.3, ec='salmon')
plt.text(2.5, -56349+17, 
         'AME 2020 value', color='salmon', ha='center')#, va='top')

color_olivedrab02 = (0.88392157, 0.91137255, 0.82745098)
plt.fill_between((0.5, 4.5), ME_NSCL2022-dME_NSCL2022, ME_NSCL2022+dME_NSCL2022, 
                 color=color_olivedrab02, lw=0.3, ec='olivedrab')
plt.text(2.5, ME_NSCL2022+dME_NSCL2022+1,
         'NSCL 2022 (ToF-ICR)', color='olivedrab', ha='center')

plt.errorbar(1, -56355.7, 22.3, color='firebrick', capsize=10, marker='_') # MR-ToF
plt.errorbar(2, (m_Zn61(0.9997125272)-61)*amu_keV, yerr=dm_Zn61(0.9997125272, 0.0000000580)*amu_keV, color='steelblue', capsize=10, marker='_') # ChiSq scaled errors, all-in-one
plt.errorbar(3, (m_Zn61_Rb(0.7176782350)-61)*amu_keV, yerr=dm_Zn61_Rb(0.7176782350, 0.0000000354)*amu_keV, color='orange', capsize=10, marker='_') # 0 order, ChiSq scaled errors, all-in-one
plt.errorbar(4, (m_Zn61(0.9997125737)-61)*amu_keV, yerr=dm_Zn61(0.9997125737, 0.0000000230)*amu_keV, color='purple', capsize=10, marker='_') # errors scaled with individual redChi^2, all-in-one

plt.xticks(ticks=[1,2,3,4],
           labels=['MR-ToF\n$^{61}$Cu$^{+}, ^{42}$CaF$^{+}$',
                   'ToF-ICR\n$^{42}$CaF$^{+}$', 
                   'ToF-ICR\n$^{85}$Rb$^{+}$', 
                   'PI-ICR\n$^{42}$CaF$^{+}$',])

plt.xlim(0.5,4.5)
plt.ylabel('mass excess / keV')
plt.minorticks_on()
ax.tick_params(axis='x', which='minor', bottom=False)

plt.savefig('comparison.pdf', bbox_inches='tight', pad_inches=0.01)