Skip to content

Instantly share code, notes, and snippets.

@iphysresearch
Last active November 5, 2020 09:42
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save iphysresearch/467d529f504e0c7a90eb5c1f3960ffa4 to your computer and use it in GitHub Desktop.
Save iphysresearch/467d529f504e0c7a90eb5c1f3960ffa4 to your computer and use it in GitHub Desktop.
Python script to plot cumulative number of events vs days in O1, O2, O3 (not including breaks). See more: https://iphysresearch.github.io/blog/post/gw/cumulative_numbe_of_gwevents/
import numpy as np
import datetime
import matplotlib.pyplot as plt
gw_event = [20150914,20151012,20151226, # O1 events
20170104,20170608,20170729,20170809,20170814,20170817,20170818,20170823,
20190408,20190412,20190413,20190413,20190421,20190424,20190425,20190426,
20190503,20190512,20190513,20190514,20190517,20190519,20190521,20190521,
20190527,20190602,20190620,20190630,20190701,20190706,20190707,20190708,
20190719,20190720,20190727,20190728,20190731,20190803,20190814,20190828,
20190828,20190909,20190910,20190915,20190924,20190929,20190930, # gap between O3a and O3b
20191105,20191109,20191129,20191204,20191205,20191213,20191215,
20191216,20191222,20200105,20200112,20200114,20200115,20200128,20200129,
20200208,20200213,20200219,20200224,20200225,20200302,20200311,20200316];
assert sorted(gw_event) == gw_event
datetime_event = [datetime.datetime.strptime(str(t), '%Y%m%d') for t in gw_event]
# Know more about Format Codes? see: https://docs.python.org/zh-cn/3/library/datetime.html#strftime-strptime-behavior
num_event = len(datetime_event)
print('Total number of GW events:', num_event)
O1_start = datetime.datetime(2015,9,12)
O1_end = datetime.datetime(2016,1,19)
len_O1 = O1_end - O1_start
O2_start = datetime.datetime(2016,11,30)
O2_end = datetime.datetime(2017,8,25)
len_O2 = O2_end - O2_start
O3a_start = datetime.datetime(2019,4,1)
O3a_end = datetime.datetime(2019,9,30)
len_O3a = O3a_end - O3a_start
O3b_start = datetime.datetime(2019,11,1)
O3b_end = datetime.datetime(2020,4,30)
len_O3b = O3b_end - O3b_start
total_days = len_O1 + len_O2 + len_O3a + len_O3b
O1 = len_O1
O2 = len_O1 + len_O2
O3a = len_O1 + len_O2 + len_O3a
O3b = len_O1 + len_O2 + len_O3a + len_O3b
nev_O1 = sum((O1_start <= np.asarray(datetime_event)) & (np.asarray(datetime_event) <= O1_end))
nev_O2 = sum((O2_start <= np.asarray(datetime_event)) & (np.asarray(datetime_event) <= O2_end))
nev_O3a = sum((O3a_start <= np.asarray(datetime_event)) & (np.asarray(datetime_event) <= O3a_end))
nev_O3b = sum((O3b_start <= np.asarray(datetime_event)) & (np.asarray(datetime_event) <= O3b_end))
assert num_event == nev_O1 + nev_O2 + nev_O3a + nev_O3b
print('Total of days:', total_days.days)
print('Number of days in O1/O2/O3a/O3b:', '{}/{}/{}/{}'.format(len_O1.days,len_O2.days,len_O3a.days,len_O3b.days))
print('Number of events in O1/O2/O3a/O3b:', '{}/{}/{}/{}'.format(nev_O1, nev_O2, nev_O3a, nev_O3b))
plt.figure(figsize=(8,4))
plt.hist(datetime_event, bins=73, histtype='step', cumulative=True, color='k', linewidth=2)
plt.ylabel('Cumulative #Events/Candidates')
plt.fill_betweenx([-1,80], O1_start,O1_end , color=[230/256,179/255,179/255])
plt.fill_betweenx([-1,80], O2_start,O2_end , color=[179/256,230/255,181/255])
plt.fill_betweenx([-1,80], O3a_start,O3a_end , color=[179/256,179/255,228/255])
plt.fill_betweenx([-1,80], O3b_start,O3b_end , color=[255/256,179/255,84/255])
plt.ylim(-1,80)
plt.xlim(O1_start, O3b_end)
plt.title('''Cumulative Count of Events and (non-retracted Alerts)
O1={}, O2={}, O3a={}, O3b={}, Total={}'''.format(nev_O1,nev_O2,nev_O3a,nev_O3b,num_event))
plt.savefig('cumulative_events_by_date.png', dpi=300, bbox_inches='tight')
plt.show()
def getstart(t):
if (O1_start <= t) & (t <= O1_end):
return O1_start
elif (O2_start <= t) & (t <= O2_end):
return O2_start - O1
elif (O3a_start <= t) & (t <= O3a_end):
return O3a_start - O2
elif (O3b_start <= t) & (t <= O3b_end):
return O3b_start - O3a
else:
raise
x = [ (t-getstart(t)).days for t in datetime_event]
y = range(len(datetime_event))
plt.figure(figsize=(7,5))
plt.plot(x, y, drawstyle='steps-post', color='k', linewidth=2)
plt.fill_betweenx([-1,80], 0, O1.days , color=[230/256,179/255,179/255])
plt.fill_betweenx([-1,80], O1.days, O2.days , color=[179/256,230/255,181/255])
plt.fill_betweenx([-1,80], O2.days, O3a.days , color=[179/256,179/255,228/255])
plt.fill_betweenx([-1,80], O3a.days, O3b.days , color=[255/256,179/255,84/255], )
plt.ylim(-1,80)
plt.xlim(0, O3b.days)
plt.ylabel('Cumulative #Events/Candidates')
plt.xlabel('Time (Days)')
plt.title('''Cumulative Count of Events and (non-retracted Alerts)
O1={}, O2={}, O3a={}, O3b={}, Total={}'''.format(nev_O1,nev_O2,nev_O3a,nev_O3b,num_event))
plt.text(O1.days*0.3, num_event*0.6, 'O1', fontsize=15)
plt.text(O1.days+(O2.days-O1.days)*0.4, num_event*0.6, 'O2', fontsize=15)
plt.text(O2.days+(O3a.days-O2.days)*0.3, num_event*0.6, 'O3a', fontsize=15)
plt.text(O3a.days+(O3b.days-O3a.days)*0.3, num_event*0.6, 'O3b', fontsize=15)
plt.savefig('cumulative_events_by_days.png', dpi=300, bbox_inches='tight')
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment