累积引力波事件率图的 python 实现

An optional description of the image for screen readers. Image credit: LIGO-G2001862. Reproduced by Herb.

数据来源 (Data Sources):

LIGO-Virgo Cumulative Event Rate Plot O1-O3a (added a separate plot to include O3b public alerts)

Note that G1901322 contains confirmed events and public alerts. And original cumulative event rate plot was coded in MATLAB.


Loading packages and data

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)

Total number of GW events: 73

Pre-processing the data for drawing

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))

Total of days: 760

Number of days in O1/O2/O3a/O3b: 129/268/182/181

Number of events in O1/O2/O3a/O3b: 3/8/39/23

Figure: Cumulative Count of GW Events by dates.

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()

Figure: Cumulative Count of GW Events by days of running for each event.

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()

Next
Previous

Related