1# (c) Stefan Countryman 2016-2018 

2 

3""" 

4FileHandlers for the MATLAB coincident analysis code. Sees whether neutrinos 

5and gravitational waves overlap and plots join skymaps. Deprecated. 

6""" 

7 

8from abc import abstractproperty 

9import logging 

10from llama.filehandler import JSONFile 

11from llama.files import matlab 

12from llama.files.skymap_info import SkymapInfo 

13from llama.files.lvc_skymap_mat import LvcSkymapMat 

14from llama.files.i3.json import IceCubeNeutrinoList 

15from llama.files.team_receipts import TeamEmailReceipt 

16from llama.files.gracedb import GraceDBReceipt 

17from llama.files.gwastro import GWAstroReceipt 

18LOGGER = logging.getLogger(__name__) 

19 

20 

21# The following class is abstract, so we shouldn't need to define abstract 

22# methods. 

23# pylint: disable=abstract-method 

24@matlab.MATLABCallingFileHandler.set_class_attributes 

25class CoincAnalysis(matlab.MATLABCallingFileHandler): 

26 """ 

27 An abstract FileHandler class for files generated as part of the 

28 SingleGWHENAnalysis pipeline. 

29 """ 

30 

31 @property 

32 def matlab_command(self): 

33 """ 

34 Use MATLAB to run Single_GWHEN_Analysis for this event. 

35 """ 

36 # get the event corresponding to this file handler. 

37 gpstime = SkymapInfo(self).event_time_gps 

38 far = SkymapInfo(self).far 

39 cmdformat = "Single_GWHEN_Analysis('{}', 'icecube', {}, {}, '{}')" 

40 return cmdformat.format(self.eventid, gpstime, far, self.rundir) 

41 

42 # use '-nojvm' option for a faster MATLAB startup time 

43 matlab_options = ['-nojvm'] 

44 

45 DEPENDENCIES = (SkymapInfo, LvcSkymapMat, IceCubeNeutrinoList) 

46 

47 

48@CoincAnalysis.sync_shared_manifests 

49class IceCubeNeutrinoListMat(CoincAnalysis): 

50 """ 

51 A MATLAB-formatted file containing Neutrino data. Used for running 

52 analysis and for generating coincidence plots. 

53 """ 

54 

55 FILENAME = "IceCubeNeutrinoList.mat" 

56 

57 

58@GWAstroReceipt.upload_this() 

59@TeamEmailReceipt.upload_this() 

60@CoincAnalysis.sync_shared_manifests 

61class CoincAnalysisInitialIcecubeJson(CoincAnalysis, JSONFile): 

62 """ 

63 A JSON document containing the quantitative results of the LLAMA 

64 coincident analysis, as documented in arXiv:1112.1140. Uses the old O2 

65 MATLAB implementation. 

66 """ 

67 

68 FILENAME = "coinc_analysis_initial_icecube.json" 

69 

70 @property 

71 def reconstructed_neutrinos(self): 

72 """ 

73 Read a list of reconstructed neutrinos from the JSON file output by 

74 the pipeline. If there are no coincident neutrinos, the list will be 

75 empty. 

76 """ 

77 data = self.read_json()['Output'] 

78 # because MATLAB saves 1x1 matrices as scalars, we have to explicitly 

79 # check whether we are dealing with a list of reconstructed neutrinos 

80 # or just a single scalar value. We will explicitly convert scalars 

81 # to 1x1 lists. 

82 if data['L'] is None: 

83 return [] 

84 elif not isinstance(data['L'], list): 

85 likelihoods = [data['L']] 

86 p_combined = [data['p_combined']] 

87 reco_dec = [data['reconstructed']['DEC']] 

88 reco_ra = [data['reconstructed']['RA']] 

89 energy_gev = [data['HEN']['E']] 

90 sigma_deg = [data['HEN']['sigma']] 

91 gps_time = [data['HEN']['t']] 

92 deltat = [gps_time[0] - data['GW']['time']] 

93 else: 

94 likelihoods = data['L'] 

95 p_combined = data['p_combined'] 

96 reco_dec = data['reconstructed']['DEC'] 

97 reco_ra = data['reconstructed']['RA'] 

98 energy_gev = data['HEN']['E'] 

99 sigma_deg = data['HEN']['sigma'] 

100 gps_time = data['HEN']['t'] 

101 deltat = [gps_time[i] - data['GW']['time'][i] 

102 for i in range(len(gps_time))] 

103 coincident_neutrino_inds = [] 

104 for i, val in enumerate(likelihoods): 

105 if val != 0: 

106 coincident_neutrino_inds.append(i) 

107 return [{ 

108 'L': likelihoods[i], 

109 'p_combined': p_combined[i], 

110 'DEC': reco_dec[i], 

111 'RA': reco_ra[i], 

112 'energy_gev': energy_gev[i], 

113 'sigma_deg': sigma_deg[i], 

114 'gps_time': gps_time[i], 

115 'dt': deltat[i] 

116 } for i in coincident_neutrino_inds] 

117 

118 @property 

119 def best_neutrino(self): 

120 """ 

121 Read the list of reconstructed neutrinos from the JSON file output 

122 by the pipeline and pick the best reconstructed neutrino. If there are 

123 no neutrinos, this property will equal None. 

124 """ 

125 neutrinos = self.reconstructed_neutrinos 

126 if not neutrinos: 

127 return None 

128 # sort on joint likelihood 

129 return sorted(neutrinos, key=lambda a: a['L'], reverse=True)[0] 

130 

131 

132@CoincAnalysis.sync_shared_manifests 

133class HENlistONSOURCEIceCubeMat(CoincAnalysis): 

134 """ 

135 A PDF plot of the LVC GW skymap with neutrinos provided by ICECUBE 

136 superimposed as point-sources. Uses the environmental variable 

137 MATLAB_EXECUTABLE_PATH to find a usable MATLAB binary for that part 

138 of the analysis. Default path subject to change. 

139 """ 

140 

141 FILENAME = "HENlist_ONSOURCE_icecube.mat" 

142 

143 

144# The following class is abstract, so we shouldn't need to define abstract 

145# methods. 

146# pylint: disable=abstract-method 

147class CoincSkymapInitialIcecube(matlab.MATLABCallingFileHandler): 

148 """ 

149 Plots of the LVC GW skymap with neutrinos provided by ICECUBE 

150 superimposed as point-sources. You can manually override the latitude and 

151 longitude limits by changing those variables' values before generating. 

152 """ 

153 DEPENDENCIES = (IceCubeNeutrinoListMat, LvcSkymapMat) 

154 

155 # specify the latitude and longitude limits for the plot in degrees. A 

156 # full skymap plot will have latitude limits of [-90, 90] and longitude 

157 # limits of [-180, 180]. if this is modified to be a programmatically 

158 # calculated property in the future, make sure to update this class's 

159 # docstring. 

160 latitude_limits = [-90, 90] 

161 longitude_limits = [-180, 180] 

162 

163 @property 

164 def matlab_command(self): 

165 """ 

166 Use MATLAB to run PlotGWHEN_5 for this event. 

167 """ 

168 neutrinodata = IceCubeNeutrinoListMat(self).fullpath 

169 gwdata = LvcSkymapMat(self).fullpath 

170 outfile = self.fullpath 

171 cmdformat = "PlotGWHEN_5('{}', '{}', '{}', {}, {})" 

172 return cmdformat.format(neutrinodata, gwdata, outfile, 

173 self.latitude_limits, self.longitude_limits) 

174 

175 

176@GWAstroReceipt.upload_this() 

177@GraceDBReceipt.upload_this(lambda s: ("LLAMA ICECUBE neutrino coincidence " 

178 "followup analysis using " 

179 f"{SkymapInfo(s).skymap_filename}")) 

180@CoincSkymapInitialIcecube.set_class_attributes 

181class CoincSkymapO2Large(CoincSkymapInitialIcecube): 

182 """ 

183 A larger version of the MATLAB-generated GW/HEN Coincident skymap. 

184 """ 

185 

186 FILENAME = "coinc_o2.png"