Fix cases where ifar limit was not being applied (#4574)

* Fix cases where ifar limit was not being applied

* some more missed cases / long lines

* redoing a couple more long lines

* missed one more place

* unrelated bug where the IFOs were not in the right order for use with the significance_dict

bin/all_sky_search/pycbc_add_statmap

@@ -361,8 +361,14 @@ fg_fars_out = np.sum(isincombo_mask * fg_fars, axis=0)
 fg_fars_exc_out = np.sum(isincombo_mask * fg_fars_exc, axis=0)
 
 # Apply any limits as appropriate
-fg_fars_out = significance.apply_far_limit(fg_fars_out, significance_dict, combo=fg_coinc_type)
-fg_fars_exc_out = significance.apply_far_limit(fg_fars_exc_out, significance_dict, combo=fg_coinc_type)
+fg_fars_out = significance.apply_far_limit(
+    fg_fars_out,
+    significance_dict,
+    combo=fg_coinc_type)
+fg_fars_exc_out = significance.apply_far_limit(
+    fg_fars_exc_out,
+    significance_dict,
+    combo=fg_coinc_type)
 
 fg_ifar = conv.sec_to_year(1. / fg_fars_out)
 fg_ifar_exc = conv.sec_to_year(1. / fg_fars_exc_out)
@@ -562,6 +568,7 @@ while True:
     final_combined_fg = final_combined_fg + \
                             combined_fg_data.select(where_combined)
     combined_fg_data = combined_fg_data.remove(where_combined)
+    fg_coinc_type = np.delete(fg_coinc_type, where_combined)
     n_triggers -= 1
 
     logging.info('Removing background triggers at time {} within window '
@@ -605,6 +612,17 @@ while True:
             sep_bg_data[key].data['decimation_factor'],
             bg_t_y,
             **significance_dict[ifo_combo_key])
+        fg_far = significance.apply_far_limit(
+            fg_far,
+            significance_dict,
+            combo=key,
+        )
+        bg_far = significance.apply_far_limit(
+            bg_far,
+            significance_dict,
+            combo=key,
+        )
+
         sep_bg_data[key].data['ifar'] = 1. / bg_far
         sep_fg_data[key].data['ifar'] = 1. / fg_far
         sep_fg_data[key].data['fap'] = 1 - \
@@ -631,9 +649,15 @@ while True:
     isincombo_mask = np.array([list(is_in_combo_time[ct])
                                for ct in all_ifo_combos])
     fg_fars = np.array([list(far[ct]) for ct in all_ifo_combos])
+    fg_fars_out = np.sum(isincombo_mask * fg_fars, axis=0)
+    fg_fars_out = significance.apply_far_limit(
+        fg_fars_out,
+        significance_dict,
+        combo=fg_coinc_type,
+    )
     # Combine the FARs with the mask to obtain the new ifars
     combined_fg_data.data['ifar'] = conv.sec_to_year(
-        1. / np.sum(isincombo_mask * fg_fars, axis=0))
+        1. / fg_fars_out)
     fg_time -= args.cluster_window
     combined_fg_data.data['fap'] = 1 - \
         np.exp(-conv.sec_to_year(fg_time) / combined_fg_data.data['ifar'])

bin/all_sky_search/pycbc_coinc_statmap

@@ -257,10 +257,22 @@ bg_far_exc, fg_far_exc = significance.get_far(
     background_time_exc,
     **significance_dict[ifo_combo])
 
-fg_far = significance.apply_far_limit(fg_far, significance_dict, combo=ifo_combo)
-bg_far = significance.apply_far_limit(bg_far, significance_dict, combo=ifo_combo)
-fg_far_exc = significance.apply_far_limit(fg_far_exc, significance_dict, combo=ifo_combo)
-bg_far_exc = significance.apply_far_limit(bg_far_exc, significance_dict, combo=ifo_combo)
+fg_far = significance.apply_far_limit(
+    fg_far,
+    significance_dict,
+    combo=ifo_combo)
+bg_far = significance.apply_far_limit(
+    bg_far,
+    significance_dict,
+    combo=ifo_combo)
+fg_far_exc = significance.apply_far_limit(
+    fg_far_exc,
+    significance_dict,
+    combo=ifo_combo)
+bg_far_exc = significance.apply_far_limit(
+    bg_far_exc,
+    significance_dict,
+    combo=ifo_combo)
 
 f['background/ifar'] = conv.sec_to_year(1. / bg_far)
 f['background_exc/ifar'] = conv.sec_to_year(1. / bg_far_exc)
@@ -421,6 +433,17 @@ while numpy.any(ifar_foreground >= background_time):
         return_counts=True,
         **significance_dict[ifo_combo])
 
+    fg_far = significance.apply_far_limit(
+        fg_far,
+        significance_dict,
+        combo=ifo_combo
+    )
+    bg_far = significance.apply_far_limit(
+        bg_far,
+        significance_dict,
+        combo=ifo_combo,
+    )
+
     # Update the ifar_foreground criteria depending on whether foreground
     # triggers are being removed via inclusive or exclusive background.   
     if args.hierarchical_removal_against == 'inclusive':
@@ -435,6 +458,11 @@ while numpy.any(ifar_foreground >= background_time):
             exc_zero_trigs.decimation_factor,
             background_time_exc,
             **significance_dict[ifo_combo])
+        fg_far_exc = significance.apply_far_limit(
+            fg_far_exc,
+            significance_dict,
+            combo=ifo_combo
+        )
         ifar_foreground = 1. / fg_far_exc
     # ifar_foreground has been updated and the code can continue.
 

bin/all_sky_search/pycbc_coinc_statmap_inj

@@ -37,12 +37,21 @@ if args.verbose:
     log_level = logging.INFO
     logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level)
 
-ifo_key = ''.join(args.ifos)
-significance_dict = significance.digest_significance_options([ifo_key], args)
 
 window = args.cluster_window
 logging.info("Loading coinc zerolag triggers")
 zdata = pycbc.io.MultiifoStatmapData(files=args.zero_lag_coincs, ifos=args.ifos)
+
+if 'ifos' in zdata.attrs:
+    ifos = zdata.attrs['ifos'].split(' ')
+    logging.info('using ifos from file {}'.format(args.zero_lag_coincs[0]))
+else:
+    ifos = args.ifos
+    logging.info('using ifos from command line input')
+
+ifo_key = ''.join(ifos)
+significance_dict = significance.digest_significance_options([ifo_key], args)
+
 zdata = zdata.cluster(window)
 
 f = h5py.File(args.output_file, "w")
@@ -51,7 +60,7 @@ f.attrs['num_of_ifos'] = zdata.attrs['num_of_ifos']
 f.attrs['pivot'] = zdata.attrs['pivot']
 f.attrs['fixed'] = zdata.attrs['fixed']
 f.attrs['timeslide_interval'] = zdata.attrs['timeslide_interval']
-f.attrs['ifos'] = ' '.join(sorted(args.ifos))
+f.attrs['ifos'] = ' '.join(sorted(ifos))
 
 # Copy over the segment for coincs and singles
 for key in zdata.seg.keys():
@@ -90,7 +99,11 @@ if len(zdata) > 0:
         background_time,
         **significance_dict[ifo_key])
 
-    fg_far_exc = significance.apply_far_limit(fg_far_exc, significance_dict, combo=ifo_key)
+    fg_far_exc = significance.apply_far_limit(
+        fg_far_exc,
+        significance_dict,
+        combo=ifo_key,
+    )
 
     ifar_exc = 1. / fg_far_exc
     fap_exc = 1 - numpy.exp(- coinc_time / ifar_exc)

bin/all_sky_search/pycbc_sngls_statmap

@@ -146,8 +146,16 @@ bg_far, fg_far = significance.get_far(
     fg_time,
     **significance_dict[ifo])
 
-fg_far = significance.apply_far_limit(fg_far, significance_dict, combo=ifo)
-bg_far = significance.apply_far_limit(bg_far, significance_dict, combo=ifo)
+fg_far = significance.apply_far_limit(
+    fg_far,
+    significance_dict,
+    combo=ifo,
+)
+bg_far = significance.apply_far_limit(
+    bg_far,
+    significance_dict,
+    combo=ifo,
+)
 
 bg_ifar = 1. / bg_far
 fg_ifar = 1. / fg_far
@@ -341,6 +349,18 @@ while numpy.any(ifar_louder > hier_ifar_thresh_s):
         fg_time,
         **significance_dict[ifo])
 
+    fg_far = significance.apply_far_limit(
+        fg_far,
+        significance_dict,
+        combo=ifo,
+    )
+
+    bg_far = significance.apply_far_limit(
+        bg_far,
+        significance_dict,
+        combo=ifo,
+    )
+
     bg_ifar = 1. / bg_far
     fg_ifar = 1. / fg_far
 
@@ -359,6 +379,12 @@ while numpy.any(ifar_louder > hier_ifar_thresh_s):
             fg_time_exc,
             **significance_dict[ifo])
 
+        fg_far_exc = significance.apply_far_limit(
+            fg_far_exc,
+            significance_dict,
+            combo=ifo,
+        )
+
         fg_ifar_exc = 1. / fg_far_exc
 
         ifar_louder = fg_ifar_exc

bin/all_sky_search/pycbc_sngls_statmap_inj

@@ -117,8 +117,14 @@ bg_far_exc, fg_far_exc = significance.get_far(
     fg_time_exc,
     **significance_dict[ifo])
 
-fg_far_exc = significance.apply_far_limit(fg_far_exc, significance_dict, combo=ifo)
-bg_far_exc = significance.apply_far_limit(bg_far_exc, significance_dict, combo=ifo)
+fg_far_exc = significance.apply_far_limit(
+    fg_far_exc,
+    significance_dict,
+    combo=ifo)
+bg_far_exc = significance.apply_far_limit(
+    bg_far_exc,
+    significance_dict,
+    combo=ifo)
 
 fg_ifar_exc = 1. / fg_far_exc
 bg_ifar_exc = 1. / bg_far_exc