O4 realtime (#51)

* load up to 5000s to check for events in realtime

* add internal public webpage fxn

* event time def earlier

* get gcn in case p<0.1

* try for pickle load

* index at 1 instead of 0 for events on skymap

* add a try to convert from moc

* add email/sms notif

* internal alerts pub webpage updates

* formatting for alerts webpage

* increment versioning

* dont send most prob direction if no coinc events

* change header for public wp

* make plots w file modification time

* precomp gw bg trials

* fix document, slack poster

* print central energy range min,max

* fix documenting bugs

* send llama results for 1det event

* adding config for sender

* add a version of combine_results run on classic

* remove uneeded, switch to rel paths

* no config needed, switch back to kafka

* stop it timing out while waiting for data

* correct error processing

* add nside param

* generate high sig gcn templ, calc per event p when gt 0.1

* add note about plot in pub wp

* get type of event and pass to make_call

* should be a float for probs

* fix bug in 2d webpage link

* fix paths

* fix paths, handle llama res if no ontime neutrinos

* update sender to match schema

* update with new key names

* add calls to shifters if results not found; 35 min timeout

* ra_unc is a list

* move incoming alert printing so rev1 doesn't print

* seperate alert and gw trials scripts

* update submission script for both

* call shifters if results not found if not a mock

* remove try/except from debugging

* increment version

Author: jessiethw

fast_response/FastResponseAnalysis.py

@@ -995,6 +995,7 @@ def make_dNdE(self):
 
         self.energy_range = (np.min([low_5_min_dec, low_5_max_dec]),
                              np.max([high_5_min_dec, high_5_max_dec]))
+        print('90% Central Energy Range: {}, {} GeV'.format(round(self.energy_range[0]), round(self.energy_range[1])))
         self.save_items['energy_range'] = self.energy_range
 
     def ns_scan(self):
@@ -1246,6 +1247,7 @@ def make_dNdE(self):
         plt.legend(loc=4, fontsize=18)
         plt.savefig(self.analysispath + '/central_90_dNdE.png',bbox_inches='tight')
 
+        print('90% Central Energy Range: {}, {} GeV'.format(round(low_5), round(high_5)))
         self.save_items['energy_range'] = (self.low5, self.high5)
 
     def write_circular(self):

fast_response/GWFollowup.py

@@ -326,12 +326,12 @@ def write_circular(self):
         except:
             noticeID = 'NOTICEID'
 
-        if pvalue > 0.1:
+        if pvalue > 0.2:
             template_path = os.path.join(base, 'circular_templates/gw_gcn_template_low.txt')
         else:
             template_path = os.path.join(base, 'circular_templates/gw_gcn_template_high.txt')
 
-        if pvalue>0.1:
+        if pvalue>0.2:
             with open(template_path, 'r') as gcn_template:
 
                 gcn = gcn_template.read()
@@ -568,7 +568,7 @@ def per_event_pvalue(self):
             usemask=False
         )
 
-        if self.p < 0.05:
+        if self.p < 0.1:
             for i in range(self.events_rec_array.size):
                 ts, p = self.per_event_scan(self.events_rec_array[i])
                 self.events_rec_array['pvalue'][i] = p

fast_response/MonitoringAndMocks/Data_Display.py

@@ -22,7 +22,7 @@
 warnings.filterwarnings('ignore', module='astropy._erfa')
 
 def dial_up(who="jessie"):
-    cell_tower = "/cvmfs/icecube.opensciencegrid.org/users/jthwaites/"
+    cell_tower = "/home/jthwaites/private/"
     #halp = "https://icecube.wisc.edu/~jthwaites/FastResponse/error_call.xml"
     subprocess.call([cell_tower+"make_call.py", f"--{who}=True", '--troubleshoot=True'])#, "--call_file", halp])
     #print([cell_tower+"make_call.py", f"--{who}=True", '--troubleshoot=True'])
@@ -747,11 +747,8 @@ def make_bg_pval_dist(fontsize=15, lower_y_bound=-3.5):
     print('Loading %i mocks (may take a while)'%(len(saved_mock_pkl)))
     for mock in saved_mock_pkl:
         with open(mock,'rb') as f:
-            try:
-                result=pickle.load(f)
-            except:
-                print('skipped {}'.format(mock))
-                continue
+            result=pickle.load(f)
+            #print('skipped {}'.format(mock))
             all_mocks[result['name']]=result['p']
     print('Done loading mocks.')
 

fast_response/listeners/gcn_listener.py

@@ -13,8 +13,6 @@
         gcn.notice_types.ICECUBE_CASCADE)
 
 def process_gcn(payload, root):
-    print("INCOMING ALERT: ",datetime.utcnow())
-
     analysis_path = os.environ.get('FAST_RESPONSE_SCRIPTS')
     if analysis_path is None:
         try:
@@ -37,19 +35,22 @@ def process_gcn(payload, root):
     stream = params['Stream']
     eventtime = root.find('.//ISOTime').text
     if stream == '26':
+        print("INCOMING ALERT: ",datetime.utcnow())
         print("Detected cascade type alert, running cascade followup. . . ")
         alert_type='cascade'
         event_name='IceCube-Cascade_{}{}{}'.format(eventtime[2:4],eventtime[5:7],eventtime[8:10])
 
         skymap = params['skymap_fits']
     else:
-        print("Found track type alert, running track followup. . . ")
         alert_type='track'
         event_name='IceCube-{}{}{}'.format(eventtime[2:4],eventtime[5:7],eventtime[8:10]) 
 
         # IceCube sends 2: a notice and a revision, only want to run once
         if int(params['Rev']) !=0:
-            return        
+            return
+        
+        print("INCOMING ALERT: ",datetime.utcnow())
+        print("Found track type alert, running track followup. . . ")
 
     event_id = params['event_id']
     run_id = params['run_id']
@@ -111,14 +112,17 @@ def process_gcn(payload, root):
         '--suffix={}'.format(suffix)]
         )
 
+    event_name=event_name+suffix
+    doc = False
     if args.document:
         try:
-            dir_1000 = glob.glob(os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'),
+            dir_1000 = glob(os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'),
                                           '*{}_1.0e+03_s').format(event_name))
             subprocess.call([analysis_path+'document.py', '--path', dir_1000[0]])
-            dir_2d = glob.glob(os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'),
+            dir_2d = glob(os.path.join(os.environ.get('FAST_RESPONSE_OUTPUT'),
                                           '*{}_1.7e+05_s').format(event_name))
             subprocess.call([analysis_path+'document.py', '--path', dir_2d[0]])
+            doc=True
         except:
             print('Failed to document to private webpage')
 
@@ -130,13 +134,18 @@ def process_gcn(payload, root):
             if shifters['start'][i]<datetime.utcnow()<shifters['stop'][i]:
                 on_shift+='<@{}> '.format(shifters['slack_id'][i])
         link = 'https://user-web.icecube.wisc.edu/~jthwaites/FastResponse/webpage/output/'
-        wp_link_1000 = '{}{}_1.0e+03_s.html'.format(link, eventtime[0:10].replace('-','_')+'_'+event_name)
-        wp_link_2d   = '{}{}_1.7e+05_s.html'.format(link, eventtime[0:10].replace('-','_')+'_'+event_name)
-        bot.send_message(f'Done running FRA for {alert_type} alert, {event_name}.\n ' +
-                         "Results for 1000s: <{}|link>. \n ".format(wp_link_1000) +
-                         "Results for 2d: <{}|link>. \n".format(wp_link_2d) +
-                         + on_shift +'on shift',
+        wp_link_1000 = '{}{}_{}_1.0e+03_s.html'.format(link, eventtime[0:10].replace('-','_'),event_name)
+        
+        day_before = '{}'.format(int(eventtime[8:10])-1)
+        if len(day_before)==1: day_before='0'+day_before
+        str_2d = '{}_{}'.format(eventtime[0:7].replace('-','_'),day_before)
+        wp_link_2d   = '{}{}_{}_1.7e+05_s.html'.format(link, str_2d, event_name)
+        bot.send_message(f'Done running FRA for {alert_type} alert, {event_name}.\n '+ on_shift +'on shift',
                          'blanket_blob')
+        if doc:
+            bot.send_message("-Results for 1000s: <{}|link> \n-Results for 2d: <{}|link>".format(
+                              wp_link_1000, wp_link_2d),
+                             'blanket_blob')
         print(' - slack message sent \n')
     except Exception as e:
         print(e)

fast_response/listeners/gw_gcn_listener.py

@@ -61,14 +61,29 @@ def process_gcn(payload, root):
 
     print('\n' +'INCOMING ALERT FOUND: ',datetime.utcnow())
     log_file.flush()
+
+    #get type of event (burst, bbh, nsbh, bns)
+    try:
+        if params['Group'] == 'Burst': 
+            merger_type = 'Burst'
+        else:
+            k = ['BNS','NSBH','BBH']
+            probs = {j: float(params[j]) for j in k}
+            merger_type = max(zip(probs.values(), probs.keys()))[1]
+    except:
+        print('Could not determine type of event')
+        merger_type = None
 
     if root.attrib['role']=='observation' and not mock:
         ## Call everyone because it's a real event!
-        call_command=['/cvmfs/icecube.opensciencegrid.org/users/jthwaites/make_call.py']
+        call_command=['/home/jthwaites/private/make_call.py']
 
         call_args = ['--justin']
         for arg in call_args:
             call_command.append(arg+'=True')
+        if merger_type is not None:
+            call_command.append(f'--type={merger_type}')
+            
         try:
             subprocess.call(call_command)
             #print('Call here.')

fast_response/precomputed_background/glob_precomputed_trials.py

@@ -14,6 +14,8 @@
                     help='Time Window in seconds')
 parser.add_argument('--dir',type=str, default='./',
                     help='directory for where trials are, will save npz to dir+/glob_trials/')
+parser.add_argument('--nside',type=int, default=256,
+                    help='nside used when running trials (default 256)')
 args = parser.parse_args()
 
 def glob_allsky_scans(delta_t, rate, dir, low_stats=False):
@@ -24,7 +26,7 @@ def glob_allsky_scans(delta_t, rate, dir, low_stats=False):
     #jobs_per_window = {1000.: 20, 172800.: 100, 2678400.: 100}
 
     files = glob(dir+'/gw_{:.1f}_mHz_seed_*_delta_t_{:.1e}.npz'.format(rate, delta_t))
-    nside = 256
+    nside = 512
     npix = hp.nside2npix(nside)
     maps = sparse.csr_matrix((0, npix), dtype=float)
 
@@ -52,7 +54,7 @@ def glob_allsky_scans(delta_t, rate, dir, low_stats=False):
 
     return maps
 
-for rate in [6.0]:#, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
+for rate in [6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2]:
     for low_stats in [False]:#[True, False]:
         print("Rate: {} mHz, low stats: {}".format(rate, low_stats))
         maps = glob_allsky_scans(args.deltaT, rate, args.dir, low_stats=low_stats)

fast_response/precomputed_background/precompute_ts.py

@@ -9,10 +9,10 @@
 import healpy as hp
 import os, sys, argparse
 from astropy.time import Time
-#from numpy.lib.recfunctions import append_fields
+from numpy.lib.recfunctions import append_fields
 from scipy import sparse
+from glob import glob
 
-from fast_response.GWFollowup import GWFollowup
 from fast_response.AlertFollowup import AlertFollowup
 
 parser = argparse.ArgumentParser(description='Fast Response Analysis')
@@ -24,8 +24,6 @@
                         help="Expected background rate in mHz (defualt 6.4)")
 parser.add_argument('--seed', default=1, type=int, 
                         help='Unique seed for running on the cluster')
-parser.add_argument('--type', default='gw', type=str,
-                        help='type of follow up (gw or alert) to use')
 parser.add_argument('--outdir',type=str, default=None,
                         help='Output directory to save npz (default = FAST_RESPONSE_OUTPUT env variable or cwd)')
 args = parser.parse_args()
@@ -40,7 +38,8 @@
 if not os.path.exists(outdir+'/trials/'):
     os.mkdir(outdir+'/trials/')
 outdir=outdir+'/trials/'
-#skymap_files = glob('/data/ana/realtime/alert_catalog_v2/2yr_prelim/fits_files/Run13*.fits.gz')
+
+skymap_files = glob('/data/ana/realtime/alert_catalog_v2/2yr_prelim/fits_files/Run13*.fits.gz')
 
 start_mjd = 58484.0 #Jan 1, 2019
 stop_mjd = start_mjd + (args.deltaT / 86400.)
@@ -51,14 +50,9 @@
 trials_per_sig = args.ntrials
 seed_counter = args.seed
 
-#skymap required for either initialization, but not used here
-if args.type=='gw':
-    f = GWFollowup('Precompute_trials_test', '/data/user/jthwaites/o3-gw-skymaps/S191216ap.fits.gz',
-                start, stop, save=False)
-else: 
-    f = AlertFollowup('Precompute_trials_test', 
-                '/data/user/jthwaites/cascade_skymaps/hese_59546_run00135946.evt000006354173.fits',
-                start, stop, save=False)
+#skymap required for initialization, but not used here
+f = AlertFollowup('Precompute_trials_test', skymap_files[0],
+                start, stop, save=False)  
 f.llh.nbackground=args.bkg*args.deltaT/1000.
 #inj = f.initialize_injector(gamma=2.5) #just put this here to initialize f.spatial_prior
 #print f.llh.nbackground
@@ -80,9 +74,7 @@
         maps[jj, pixels] = val['TS']
 print("DONE")
 hp_sparse = maps.tocsr()
-if args.type=='gw': 
-    outfilename = outdir+'gw_{:.1f}_mHz_seed_{}_delta_t_{:.1e}.npz'.format(args.bkg, args.seed, args.deltaT)
-else:
-    outfilename = outdir+'{:.1f}_mHz_seed_{}_delta_t_{:.1e}.npz'.format(args.bkg, args.seed, args.deltaT)
+
+outfilename = outdir+'{:.1f}_mHz_seed_{}_delta_t_{:.1e}.npz'.format(args.bkg, args.seed, args.deltaT)
 sparse.save_npz(outfilename, hp_sparse)
 print("Saved to {}".format(outfilename))