Skip to content

Commit

Permalink
Make pycbc_bin_trigger_rates_dq use a uniform f_low instead of gettin…
Browse files Browse the repository at this point in the history
…g f_low from bank. Also improved code style.
  • Loading branch information
maxtrevor committed Aug 31, 2023
1 parent 86a83a1 commit 3630218
Showing 1 changed file with 29 additions and 28 deletions.
57 changes: 29 additions & 28 deletions bin/all_sky_search/pycbc_bin_trigger_rates_dq
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,13 @@ from pycbc.version import git_verbose_msg as version
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--version', action='version', version=version)
parser.add_argument('--verbose', action="store_true")
parser.add_argument("--ifo", type=str,required=True)
parser.add_argument("--ifo", type=str, required=True)
parser.add_argument("--trig-file", required=True)
parser.add_argument("--stat-threshold", type=float,
help="Only consider triggers with statistic value "
"above this threshold")
parser.add_argument("--dq-file", required=True,nargs='+')
parser.add_argument("--dq-channel", required=False,type=str,
parser.add_argument("--dq-file", required=True, nargs='+')
parser.add_argument("--dq-channel", required=False, type=str,
help='name of channel to read in from '
'provided dq file. Required if '
'provided file is .hdf5 format')
Expand All @@ -37,16 +37,17 @@ parser.add_argument("--prune-window", type=float, default=0.1,
"which is loudest in each split, default 0.1s")

pystat.insert_statistic_option_group(parser,
default_ranking_statistic='single_ranking_only')
default_ranking_statistic='single_ranking_only')
args = parser.parse_args()
pycbc.init_logging(args.verbose)

with h5.File(args.bank_file, 'r') as bank:
if args.background_bins:
logging.info('Sorting bank into bins...')
data = {'mass1': bank['mass1'][:], 'mass2': bank['mass2'][:],
'spin1z': bank['spin1z'][:], 'spin2z': bank['spin2z'][:],
'f_lower': bank['f_lower']}
'spin1z': bank['spin1z'][:], 'spin2z': bank['spin2z'][:]}
f_low = np.ones(len(data['mass1'])) * 15.
data['f_lower'] = f_low
locs_dict = pycbc.events.background_bin_from_string(
args.background_bins, data)
del data
Expand All @@ -57,11 +58,11 @@ with h5.File(args.bank_file, 'r') as bank:

logging.info('Reading trigger file...')
ifo = args.ifo
with h5.File(args.trig_file,'r') as trig_file:
with h5.File(args.trig_file, 'r') as trig_file:
trig_times = trig_file[ifo+'/end_time'][:]
trig_ids = trig_file[ifo+'/template_id'][:]

if args.stat_threshold or args.prune_number>0:
if args.stat_threshold or args.prune_number > 0:
logging.info('Calculating stat and filtering...')
rank_method = pystat.get_statistic_from_opts(args, [ifo])
stat = rank_method.get_sngl_ranking(trig_file[ifo])
Expand All @@ -71,7 +72,7 @@ with h5.File(args.trig_file,'r') as trig_file:
trig_times = trig_times[abovethresh]
stat = stat[abovethresh]
del abovethresh
if args.prune_number<1:
if args.prune_number < 1:
del stat

trig_times_int = trig_times.astype('int')
Expand All @@ -82,26 +83,26 @@ dq_logl = np.array([])
for filename in args.dq_file:
logging.info('Reading DQ file %s...', filename)
dq_data = load_timeseries(filename, group=args.dq_channel)
dq_logl = np.concatenate((dq_logl,dq_data[:]))
dq_times = np.concatenate((dq_times,dq_data.sample_times))
dq_logl = np.concatenate((dq_logl, dq_data[:]))
dq_times = np.concatenate((dq_times, dq_data.sample_times))
del dq_data

# todo: make this configurable
percent_bin = 0.5
n_bins = int(100./percent_bin)
percentiles = np.linspace(0,100,n_bins+1)
percentiles = np.linspace(0, 100, n_bins+1)
bin_times = np.zeros(n_bins)
dq_percentiles = np.percentile(dq_logl,percentiles)[1:]
dq_percentiles = np.percentile(dq_logl, percentiles)[1:]


# seconds bin tells what bin each second ends up
seconds_bin = np.array([n_bins - \
len(dq_percentiles[dq_percentiles >= dq_ll]) \
seconds_bin = np.array([n_bins -
len(dq_percentiles[dq_percentiles >= dq_ll])
for dq_ll in dq_logl]).astype('int')
del dq_percentiles

# bin times tells how much time ends up in each bin
bin_times = np.array([len(seconds_bin[seconds_bin==(i)]) \
bin_times = np.array([len(seconds_bin[seconds_bin == i])
for i in range(n_bins)]).astype('float')
full_time = float(len(seconds_bin))
times_nz = (bin_times > 0)
Expand All @@ -111,7 +112,7 @@ del dq_logl
dq_percentiles_time = dict(zip(dq_times, seconds_bin*percent_bin/100))
del dq_times

if args.prune_number>0:
if args.prune_number > 0:
for bin_name in locs_names:
logging.info('Processing bin %s...', bin_name)
bin_locs = locs_dict[bin_name]
Expand All @@ -121,9 +122,9 @@ if args.prune_number>0:
for j in range(args.prune_number):
max_stat_arg = np.argmax(trig_stats_bin)
remove = np.nonzero(abs(trig_times_bin[max_stat_arg] - trig_times)
< args.prune_window)[0]
remove_inbin = np.nonzero(abs(trig_times_bin[max_stat_arg] \
- trig_times_bin) < args.prune_window)[0]
< args.prune_window)[0]
remove_inbin = np.nonzero(abs(trig_times_bin[max_stat_arg]
- trig_times_bin) < args.prune_window)[0]
stat[remove] = 0
trig_stats_bin[remove_inbin] = 0
keep = np.nonzero(stat)[0]
Expand All @@ -134,27 +135,27 @@ if args.prune_number>0:

del trig_times

with h5.File(args.output_file,'w') as f:
with h5.File(args.output_file, 'w') as f:
for bin_name in locs_names:
bin_locs = locs_dict[bin_name]
trig_times_bin = trig_times_int[np.isin(trig_ids, bin_locs)]
trig_percentile = np.array([dq_percentiles_time[t] \
trig_percentile = np.array([dq_percentiles_time[t]
for t in trig_times_bin])
logging.info('Processing %d triggers...', len(trig_percentile))
(counts, bins) = np.histogram(trig_percentile, bins = (percentiles)/100)

(counts, bins) = np.histogram(trig_percentile, bins=(percentiles)/100)
counts = counts.astype('float')
rates = np.zeros(len(bin_times))
rates[times_nz] = counts[times_nz]/bin_times[times_nz]
mean_rate = len(trig_percentile) / full_time
if mean_rate > 0.:
rates = rates / mean_rate

logging.info('Writing rates to output file %s...', args.output_file)
grp = f.create_group(bin_name)
grp['rates']=rates
grp['locs']=locs_dict[bin_name]
grp['rates'] = rates
grp['locs'] = locs_dict[bin_name]

f.attrs['names'] = locs_names

logging.info('Done!')

0 comments on commit 3630218

Please sign in to comment.