diff --git a/README.md b/README.md index 1dce382..2dc83a9 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,7 @@ In order to run this analysis, you should find the alert event time (sent in the Once you have these details, you can run the analysis with ```console -python run_track_followup.py --skymap=/home/followup/output_plots/run{RUNID}.evt{EVENTID}.HESE.skymap_nside_512.fits.gz --time={ALERT_MJD} --gcn_notice_num={GCN_CIRCULAR_NUMBER} --alert_id={RUNID}:{EVENTID} +python run_track_followup.py --skymap=/home/followup/output_plots/run{RUNID}.evt{EVENTID}.neutrino_8_64_512_1024.skymap_nside_1024.fits.gz --time={ALERT_MJD} --alert_circ={GCN_CIRCULAR_NUMBER} --alert_id={RUNID}:{EVENTID} ``` This will run two analyses, one with a time window of 1000 seconds and one with a time window of 2 days, both centered on the alert time. It will also remove the alert event from the sample, and it will assume a spectral index of -2.5, which is different than the -2 used for the `run_external_followup.py` script. @@ -141,7 +141,7 @@ To run the analysis to follow up a graviational wave event, you will need to nav As an example, there is a sample map included here: `fast_response/sample_skymaps/S191216ap_update.xml`. To run this event, the input would look like: ```console -python run_gw_followup.py --name="S191216ap Update" --time=58833.893 --skymap="https://gracedb.ligo.org/api/superevents/S191216ap/files/LALInference.fits.gz,0" +python run_gw_followup.py --name="S191216ap Update" --time=58833.89836 --skymap="https://gracedb.ligo.org/api/superevents/S191216ap/files/LALInference.fits.gz,0" ``` Additionally, you can choose to pass the argument `--allow_neg_ts` (bool) if you want to use the convention where a negative TS is allowed. The default is to use the convention TS>=0. diff --git a/doc/source/conf.py b/doc/source/conf.py index 1b18cfd..ea1a696 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -16,7 +16,7 @@ project = 'FastResponseAnalysis' copyright = '2023, Alex Pizzuto, Jessie Thwaites' author = 'Alex Pizzuto, Jessie Thwaites' -release = '1.2.2' +release = '1.3.0' # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration diff --git a/fast_response/AlertFollowup.py b/fast_response/AlertFollowup.py index f574778..7921e89 100644 --- a/fast_response/AlertFollowup.py +++ b/fast_response/AlertFollowup.py @@ -25,6 +25,7 @@ class AlertFollowup(PriorFollowup): _fix_index = True _float_index = not _fix_index _index = 2.5 + _llh_map = False def run_background_trials(self, ntrials = 1000): r"""For alert events with specific time windows, @@ -257,7 +258,22 @@ class TrackFollowup(AlertFollowup): """ Class for followup of track alert events By default, uses a fixed index of 2.5 in LLH. - By default, converts milipede LLH map to a PDF. + Assumes a probability skymap for the follow-up. + Built on AlertFollowup and PriorFollowup base classes + + """ + _smear = False + _dataset = "GFUOnline_v001p02" + _fix_index = True + _float_index = not _fix_index + _index = 2.5 + _llh_map = False + +class TrackFollowupLLH(AlertFollowup): + """ + [Old] Class for followup of track alert events + By default, uses a fixed index of 2.5 in LLH + and converts milipede LLH map to a PDF. Built on AlertFollowup and PriorFollowup base classes """ @@ -266,6 +282,7 @@ class TrackFollowup(AlertFollowup): _fix_index = True _float_index = not _fix_index _index = 2.5 + _llh_map = True def format_skymap(self, skymap): """ @@ -370,4 +387,5 @@ class CascadeFollowup(AlertFollowup): _dataset = "GFUOnline_v001p02" _fix_index = True _float_index = not _fix_index - _index = 2.5 \ No newline at end of file + _index = 2.5 + _llh_map = False \ No newline at end of file diff --git a/fast_response/FastResponseAnalysis.py b/fast_response/FastResponseAnalysis.py index 0cb7b3d..f535b75 100644 --- a/fast_response/FastResponseAnalysis.py +++ b/fast_response/FastResponseAnalysis.py @@ -62,8 +62,8 @@ def __init__(self, name, tstart, tstop, skipped=None, seed=None, outdir=None, save=True, extension=None): self.name = name - if seed is not None: - self.llh_seed(seed) + # if seed is not None: + # self.llh_seed(seed) if outdir is None: outdir = os.environ.get('FAST_RESPONSE_OUTPUT') if outdir is None: @@ -581,8 +581,9 @@ def plot_skymap_zoom(self, with_contour=False, contour_files=None): cont = np.loadtxt(c_file, skiprows=1) cont_ra = cont.T[0] cont_dec = cont.T[1] - label = 'Millipede 50\%, 90\% (160427A syst.)' \ - if contour_counter == 0 else '' + cont_label = 'Millipede 50\%, 90\% (160427A syst.)' if self._llh_map else 'Skymap 50\%, 90\%' + label = cont_label if contour_counter == 0 else '' + hp.projplot(np.pi/2. - cont_dec, cont_ra, linewidth=3., color='k', linestyle=cont_ls[contour_counter], coord='C', label=label) @@ -872,7 +873,7 @@ def run_background_trials(self, ntrials=1000): self.tsd = tsd self.save_items['tsd'] = tsd - def find_coincident_events(self): + def find_coincident_events(self, print_events=False): r"""Find coincident events for a skymap based analysis. These are ontime events that are also in the 90% contour of the skymap @@ -885,13 +886,16 @@ def find_coincident_events(self): overlap = np.isin(exp_pix, self.ipix_90) events = events[overlap] - # print nearby events, as a check (if needed) - # msk1 = (self.llh.exp[t_mask]['ra'] < (self.skymap_fit_ra+np.radians(5)))*(self.llh.exp[t_mask]['ra'] > (self.skymap_fit_ra-np.radians(5))) - # msk2 = (self.llh.exp[t_mask]['dec'] < (self.skymap_fit_dec+np.radians(5)))*((self.llh.exp[t_mask]['dec'] > self.skymap_fit_dec-np.radians(5))) - # msk3 = msk1*msk2 - # print('Nearby events:') - # print("[run, event, ra, dec, sigma, logE, time]") - # for e in self.llh.exp[t_mask][msk3]: print([e[k] for k in ['run', 'event', 'ra', 'dec', 'sigma', 'logE', 'time']]) + if print_events: + # print nearby events, as a check (if needed) + msk1 = (self.llh.exp[t_mask]['ra'] < (self.skymap_fit_ra+np.radians(5)))*(self.llh.exp[t_mask]['ra'] > (self.skymap_fit_ra-np.radians(5))) + msk2 = (self.llh.exp[t_mask]['dec'] < (self.skymap_fit_dec+np.radians(5)))*((self.llh.exp[t_mask]['dec'] > self.skymap_fit_dec-np.radians(5))) + msk3 = msk1*msk2 + if np.count_nonzero(msk3) > 0: + print('Events within 5 deg of best-fit:') + print("[run, event, ra, dec, sigma, logE, time]") + for e in self.llh.exp[t_mask][msk3]: + print([e[k] for k in ['run', 'event', 'ra', 'dec', 'sigma', 'logE', 'time']]) if len(events) == 0: coincident_events = [] diff --git a/fast_response/listeners/gw_gcn_listener.py b/fast_response/listeners/gw_gcn_listener.py index 6a1ccd4..ff4fc67 100755 --- a/fast_response/listeners/gw_gcn_listener.py +++ b/fast_response/listeners/gw_gcn_listener.py @@ -177,14 +177,17 @@ def process_gcn(payload, root): #'--allow_neg_ts=True'] ) - output = os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'),eventtime[0:10].replace('-','_')+'_'+name) + analysis_start = Time(event_mjd - 500./86400., format='mjd').iso + output = os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'), + analysis_start[0:10].replace('-','_')+'_'+name) #update webpages webpage_update = os.path.join(analysis_path,'document.py') if not mock and root.attrib['role'] == 'observation': try: subprocess.call([webpage_update, '--gw', f'--path={output}']) - wp_link = 'https://user-web.icecube.wisc.edu/~jthwaites/FastResponse/gw-webpage/output/{}.html'.format(eventtime[0:10].replace('-','_')+'_'+name) + wp_link = 'https://user-web.icecube.wisc.edu/~jthwaites/FastResponse/gw-webpage/output/'+\ + '{}.html'.format(analysis_start[0:10].replace('-','_')+'_'+name) slack_message = "UML GW analysis finished running for event {}: <{}|link>.".format(name, wp_link) with open('../slack_posters/internal_alert_slackbot.txt') as f: chan = f.readline().rstrip('\n') diff --git a/fast_response/make_ontime_plots.py b/fast_response/make_ontime_plots.py index 119a1b6..a07303b 100644 --- a/fast_response/make_ontime_plots.py +++ b/fast_response/make_ontime_plots.py @@ -174,7 +174,6 @@ def make_rate_plots(time_window, run_table, query_events, dirname, season='neutr format='mjd').iso #first 35 days rates1 = icecube.realtime_tools.live.get_rates(run_table[0]['start'], midpt) rates2 = icecube.realtime_tools.live.get_rates(midpt, run_table[-1]['stop']) - #avoid double counting last entry rates = np.concatenate([rates1, rates2]) else: rates = icecube.realtime_tools.live.get_rates(run_table[0]['start'], run_table[-1]['stop']) diff --git a/fast_response/reports/ReportGenerator.py b/fast_response/reports/ReportGenerator.py index 9b34397..3174eea 100644 --- a/fast_response/reports/ReportGenerator.py +++ b/fast_response/reports/ReportGenerator.py @@ -203,6 +203,8 @@ def ontime_table(query_dict): newevent = event['value']['data'] for key,val in list(newevent['reco']['splinempe'].items()): newevent['splinempe_'+key]=val + if '+' in newevent['eventtime']: + newevent['eventtime'] = newevent['eventtime'].split('+')[0] if Time(newevent['eventtime'],scale='utc',format='iso') >= Time("2018-07-10 17:52:03.34", format='iso',scale='utc'): newevent['muex'] = newevent['reco']['energy']['mpe_muex'] del newevent['reco'] diff --git a/fast_response/scripts/combine_results_kafka.py b/fast_response/scripts/combine_results_kafka.py index 0e72508..ce8f396 100644 --- a/fast_response/scripts/combine_results_kafka.py +++ b/fast_response/scripts/combine_results_kafka.py @@ -196,7 +196,7 @@ def parse_notice(record, wait_for_llama=False, heartbeat=False): return collected_results = {} - collected_results["$schema"]= "https://gcn.nasa.gov/schema/v4.0.0/gcn/notices/icecube/lvk_nu_track_search.schema.json" + collected_results["$schema"]= "https://gcn.nasa.gov/schema/v4.1.0/gcn/notices/icecube/lvk_nu_track_search.schema.json" collected_results["type"]= "IceCube LVK Alert Nu Track Search" eventtime = record.find('.//ISOTime').text @@ -235,7 +235,8 @@ def parse_notice(record, wait_for_llama=False, heartbeat=False): ) while results_done == False: - start_date = Time(dateutil.parser.parse(eventtime)).datetime + start_time = Time(Time(eventtime, format='isot').mjd - 500./86400., format='mjd').isot + start_date = Time(dateutil.parser.parse(start_time)).datetime start_str = f'{start_date.year:02d}_{start_date.month:02d}_{start_date.day:02d}' uml_results_path = os.path.join(fra_results_location, start_str + '_' + name.replace(' ', '_') \ diff --git a/fast_response/scripts/run_external_followup.py b/fast_response/scripts/run_external_followup.py index b0135e6..b43bde1 100755 --- a/fast_response/scripts/run_external_followup.py +++ b/fast_response/scripts/run_external_followup.py @@ -80,9 +80,9 @@ def run_analysis(args): type= lambda z:[ tuple(int(y) for y in x.split(':')) for x in z.split(',')], help="Event to exclude from the analyses, eg." "Example --skip-events=127853:67093193") - parser.add_argument('--ntrials', default=100, type=int, + parser.add_argument('--ntrials', default=200, type=int, help="Number of background trials to perform") - parser.add_argument('--n_per_sig', default=100, type=int, + parser.add_argument('--n_per_sig', default=200, type=int, help="Number of signal trials per injected ns to perform") parser.add_argument('--document', default=False, action='store_true') args = parser.parse_args() diff --git a/fast_response/scripts/run_track_followup.py b/fast_response/scripts/run_track_followup.py index 7960189..48cdf84 100755 --- a/fast_response/scripts/run_track_followup.py +++ b/fast_response/scripts/run_track_followup.py @@ -9,7 +9,7 @@ import os, argparse, glob from astropy.time import Time -from fast_response.AlertFollowup import TrackFollowup +from fast_response.AlertFollowup import TrackFollowup, TrackFollowupLLH import fast_response.web_utils as web_utils parser = argparse.ArgumentParser(description='Fast Response Analysis') @@ -17,7 +17,7 @@ help='path to skymap') parser.add_argument('--time', type=float, default=None, help='Time of the alert event (mjd)') -parser.add_argument('--gcn_notice_num', default=0, type=int, +parser.add_argument('--alert_circ', default=0, type=int, help="GCN Circular NUMBER for updated circular (not Bacodine/Notice)") parser.add_argument('--alert_id', default=None, type= lambda z:[ tuple(int(y) for y in x.split(':')) for x in z.split(',')], @@ -27,6 +27,10 @@ parser.add_argument('--suffix', type=str, default='A', help="letter to differentiate multiple alerts on the same day (default = A)." "Event name given by IceCube-yymmdd + suffix.") +parser.add_argument('--print_nearby', action='store_true', default=False, + help='Option to print events within 5 degrees of best-fit to the screen (default False)') +parser.add_argument('--deltallh_skymap', action='store_true', default=False, + help='Option MUST be raised if using deltaLLH map (default: False, assumes probability map)') args = parser.parse_args() track_time = Time(args.time, format='mjd') @@ -53,14 +57,15 @@ for containment in ['50', '90']: contour_fs = contour_fs + glob.glob(base_skymap_path + f'run{run_id:08d}.evt{ev_id:012d}.*.contour_{containment}.txt') - # contour_fs = [base_skymap_path \ - # + f'run{run_id:08d}.evt{ev_id:012d}.HESE.contour_{containment}.txt' \ - # for containment in ['50', '90']] - # contour_fs = [f for f in contour_fs if os.path.exists(f)] if len(contour_fs) == 0: contour_fs = None - f = TrackFollowup(name, args.skymap, start, stop, skipped=args.alert_id) + if args.deltallh_skymap: + print('Initializing Track follow up with deltaLLH map') + f = TrackFollowupLLH(name, args.skymap, start, stop, skipped=args.alert_id) + else: + print('Initializing Track follow up with probability skymap') + f = TrackFollowup(name, args.skymap, start, stop, skipped=args.alert_id) f.unblind_TS() f.plot_ontime(contour_files=contour_fs) @@ -68,12 +73,12 @@ f.make_dNdE() f.plot_tsd() f.upper_limit() - f.find_coincident_events() + f.find_coincident_events(print_events=args.print_nearby) results = f.save_results() f.generate_report() all_results[delta_t] = results -all_results[1000.]['gcn_num'] = args.gcn_notice_num +all_results[1000.]['gcn_num'] = args.alert_circ # Write circular to the output directory of the 2 day analysis f.write_circular(all_results) diff --git a/fast_response/web_utils.py b/fast_response/web_utils.py index a8a0a65..59a74fa 100644 --- a/fast_response/web_utils.py +++ b/fast_response/web_utils.py @@ -112,7 +112,7 @@ def updateDataFrame(analysis, gw=False, make_df=False): num = np.count_nonzero(df.index == analysis['name']) analysis['name'] += '_{}'.format(num) df.loc[analysis['name']] = new_list - + if gw: df.to_pickle(f'{base_path}/results_dataframe_gw.pkl') else: diff --git a/setup.py b/setup.py index f4e47f0..cb432df 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ import setuptools long_message = 'Fast Response Analysis' -version = "1.2.2" +version = "1.3.0" setuptools.setup( name="fast_response",