diff --git a/utils/__pycache__/plot.cpython-37.pyc b/utils/__pycache__/plot.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..c3c80f2961467f6844fb0ad76b0358ea4bab6608
Binary files /dev/null and b/utils/__pycache__/plot.cpython-37.pyc differ
diff --git a/utils/plot.py b/utils/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..ecd1cd6a6e60ebb471b461ddb29a172c43046872
--- /dev/null
+++ b/utils/plot.py
@@ -0,0 +1,75 @@
+import numpy as np
+import pandas as pd
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+# PLOTTING OPTIONS
+
+def remove_space(df):
+    df.columns = [c.replace(' ','') for c in df.columns]
+    return df
+
+def readlnb(path):
+    d221 = pd.read_csv(path + "/221/GW150914_PROD1_Kerr_221_0M/Nested_sampler/Evidence.txt",sep=' ')
+    d220 = pd.read_csv(path + "/220/GW150914_PROD1_Kerr_220_0M/Nested_sampler/Evidence.txt",sep=' ')
+    d221 = remove_space(d221)
+    d220 = remove_space(d220)
+    return d221['lnB'].values[0] - d220['lnB'].values[0]
+
+def read_posA1(path,lmn='221'):
+    d = pd.read_csv(path+'/'+lmn+'/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat',sep=' ')
+    d.columns = [d.columns[(i+1)%len(d.columns)] for i in range(len(d.columns))]
+    d.drop(columns='#',inplace=True)
+    d = remove_space(d)
+    return d['A2221'].values
+
+def plot_A221violin(time,rootpath):
+    '''
+    Plot A221 violin figures
+    ========
+    time: a time array
+    rootpath: the root path of PyRing outputs
+
+    '''
+    A221 = {}
+    for i,t in enumerate(time):
+        path = rootpath + '/t'+str(i)
+        A221[i] = read_posA1(path)
+    x,y = zip(*A221.items())
+
+    violinparts = ax.violinplot(y, 
+          showmeans=False, 
+          showmedians=True, 
+          showextrema=False, 
+          vert=True, 
+          widths=0.7)
+    for pc in violinparts['bodies']:
+        pc.set_facecolor('tab:blue')
+        pc.set_edgecolor('tab:blue')
+        pc.set_alpha(0.4)
+    
+    #plot error bar:
+    for i,d in enumerate(y):
+        x = np.median(d)
+        #print i,x                 
+        dx_l,dx_u = np.percentile(d,[5]),np.percentile(d,[95])
+        err_low = np.abs(x-dx_l)
+        err_high = np.abs(dx_u-x)
+        ax.errorbar(i+1,x,yerr=[err_low,err_high],ecolor='tab:blue',alpha=0.4,
+                            marker=None, capthick=2, capsize=18, linewidth=0.01)
+
+    # set x axis
+    ax.get_xaxis().set_tick_params(direction='out')
+    ax.xaxis.set_ticks_position('bottom')
+    labels = ['{:.2f}'.format(i) for i in time]
+    ax.set_xticks(np.arange(1, len(labels) + 1))
+    ax.set_xticklabels(labels)
+    ax.set_xlim(0.25, len(labels) + 0.75)
+    
+
+    for tick in ax.get_xticklabels():
+        tick.set_rotation(60)
+
+    ax.set_title('GW150914')
+    ax.set_ylabel('A_{221}')
+    ax.set_ylim(0,30)
+    ax.set_xlabel('$t-t_\mathrm{ref} (tH1=1126259462.42323)$ ms')
\ No newline at end of file