diff --git a/flarestack/core/results.py b/flarestack/core/results.py index a0795c17..7927bce6 100644 --- a/flarestack/core/results.py +++ b/flarestack/core/results.py @@ -265,7 +265,7 @@ def nu_astronomy(self, flux, e_pdf_dict): :param e_pdf_dict: Dictionary containing energy PDF information :return: Value for the neutrino luminosity """ - return calculate_astronomy(flux, e_pdf_dict, self.sources) + return calculate_astronomy(flux, e_pdf_dict) def clean_merged_data(self): """Function to clear cache of all data""" diff --git a/flarestack/core/unblinding.py b/flarestack/core/unblinding.py index 1a9a1b75..a5ddc5f8 100644 --- a/flarestack/core/unblinding.py +++ b/flarestack/core/unblinding.py @@ -237,11 +237,13 @@ def calculate_upper_limits(self): flux_uls = [] fluence_uls = [] - e_per_source_uls = [] energy_flux = [] x_axis = [] + # e_per_source_uls = [] + for subdir in sorted(os.listdir(self.pickle_dir)): + # this is a loop over different `gamma` values! root = os.path.join(self.unblind_dict["background_ts"], subdir) @@ -316,9 +318,10 @@ def calculate_upper_limits(self): fluence_uls.append(astro_dict[key] * inj_time) - e_per_source_uls.append( - astro_dict["Mean Luminosity (erg/s)"] * inj_time - ) + # this functionality has been retired + # e_per_source_uls.append( + # astro_dict["Mean Luminosity (erg/s)"] * inj_time + # ) energy_flux.append(astro_dict[key]) x_axis.append(float(subdir)) @@ -345,10 +348,10 @@ def calculate_upper_limits(self): plt.figure() ax1 = plt.subplot(111) - ax2 = ax1.twinx() + # ax2 = ax1.twinx() ax1.plot(x_axis, fluence_uls, marker="o") - ax2.plot(x_axis, e_per_source_uls, marker="o") + # ax2.plot(x_axis, e_per_source_uls, marker="o") if self.mock_unblind: ax1.text( @@ -360,16 +363,20 @@ def calculate_upper_limits(self): transform=ax1.transAxes, ) - ax2.grid(True, which="both") + ax1.grid(True, which="both") + # ax2.grid(True, which="both") ax1.set_xlabel("Gamma") - ax2.set_xlabel("Gamma") + # ax2.set_xlabel("Gamma") ax1.set_ylabel(r"Total Fluence [GeV cm$^{-2}$ s$^{-1}$]") - ax2.set_ylabel(r"Isotropic-Equivalent Luminosity $L_{\nu}$ (erg/s)") + # ax2.set_ylabel(r"Isotropic-Equivalent Luminosity $L_{\nu}$ (erg/s)") ax1.set_yscale("log") - ax2.set_yscale("log") + # ax2.set_yscale("log") + + axes = [ax1] # was: axes = [ax1, ax2] + ordinates = [fluence_uls] # was: [fluence_uls, e_per_source_uls] - for k, ax in enumerate([ax1, ax2]): - y = [fluence_uls, e_per_source_uls][k] + for k, ax in enumerate(axes): + y = ordinates[k] ax.set_ylim(0.95 * min(y), 1.1 * max(y)) plt.tight_layout() @@ -387,7 +394,7 @@ def calculate_upper_limits(self): "flux": flux_uls, "energy_flux": energy_flux, "fluence": fluence_uls, - "energy": e_per_source_uls, + # "energy": e_per_source_uls, } with open(self.limit_path, "wb") as f: diff --git a/flarestack/utils/neutrino_astronomy.py b/flarestack/utils/neutrino_astronomy.py index 8a3ac3f4..6b51eb6a 100644 --- a/flarestack/utils/neutrino_astronomy.py +++ b/flarestack/utils/neutrino_astronomy.py @@ -5,32 +5,11 @@ import math from astropy.coordinates import Distance from flarestack.core.energy_pdf import EnergyPDF -from flarestack.utils.catalogue_loader import calculate_source_weight logger = logging.getLogger(__name__) -e_0 = 1 * u.GeV -# Set parameters for conversion from CR luminosity to nu luminosity -f_pi = 0.1 -waxmann_bachall = (3.0 / 8.0) * f_pi - -f_cr_to_nu = 0.05 - - -def find_zfactor(distance): - """For a given distance, converts this distance to a redshift, and then - returns the 1+z factor for that distance - - :param distance: Astropy distance (with units) - :return: Corresponding 1+z factor - """ - redshift = astropy.coordinates.Distance(distance).compute_z() - zfactor = 1 + redshift - return zfactor - - -def calculate_astronomy(flux, e_pdf_dict, catalogue): +def calculate_astronomy(flux, e_pdf_dict): flux /= u.GeV * u.cm**2 * u.s @@ -38,7 +17,7 @@ def calculate_astronomy(flux, e_pdf_dict, catalogue): astro_res = dict() - phi_integral = energy_PDF.flux_integral() * u.GeV + # phi_integral = energy_PDF.flux_integral() * u.GeV # unused e_integral = energy_PDF.fluence_integral() * u.GeV**2 @@ -50,150 +29,12 @@ def calculate_astronomy(flux, e_pdf_dict, catalogue): logger.debug("Energy Flux:{0}".format(tot_fluence)) - src_1 = np.sort(catalogue, order="distance_mpc")[0] - - frac = calculate_source_weight(src_1) / calculate_source_weight(catalogue) - - si = flux * frac - - astro_res["Flux from nearest source"] = si - logger.debug("Total flux: {0}".format(flux)) - logger.debug("Fraction from nearest source: {0}".format(frac)) - logger.debug("Flux from nearest source: {0}".format(flux * frac)) - - lumdist = src_1["distance_mpc"] * u.Mpc - - area = 4 * math.pi * (lumdist.to(u.cm)) ** 2 - - dNdA = (si * phi_integral).to(u.s**-1 * u.cm**-2) - - # int_dNdA += dNdA - N = dNdA * area - - logger.debug("There would be {:.3g} neutrinos emitted.".format(N)) logger.debug( "The energy range was assumed to be between {0} and {1}".format( energy_PDF.integral_e_min, energy_PDF.integral_e_max ) ) - # Energy requires a 1/(1+z) factor - - zfactor = find_zfactor(lumdist) - etot = (si * area * e_integral).to(u.erg / u.s) * zfactor - - astro_res["Mean Luminosity (erg/s)"] = etot.value - - logger.debug("The required neutrino luminosity was {0}".format(etot)) - - cr_e = etot / f_cr_to_nu - - logger.debug( - "Assuming {0:.3g}% was transferred from CR to neutrinos, we would require a total CR luminosity of {1}".format( - 100 * f_cr_to_nu, cr_e - ) - ) return astro_res - - -# def calculate_neutrinos(source, season, inj_kwargs): -# -# print source -# inj = Injector(season, [source], **inj_kwargs) -# -# print "\n" -# print source["Name"], season["Name"] -# print "\n" -# -# energy_pdf = inj_kwargs["Injection Energy PDF"] -# -# energy = energy_pdf["Energy Flux"] * inj.sig_time_pdf.effective_injection_time( -# source) -# print "Neutrino Energy is", energy -# -# lumdist = source["Distance (Mpc)"] * u.Mpc -# -# print "Source is", lumdist, "away" -# -# area = 4 * math.pi * (lumdist.to(u.cm)) ** 2 -# -# nu_flu = energy.to(u.GeV)/area -# -# print "Neutrino Fluence is", nu_flu -# -# # if "E Min" in energy_pdf.keys(): -# # e_min = energy_pdf["E Min"] * u.GeV -# # else: -# # e_min = (100 * u.GeV) -# # -# # if "E Max" in energy_pdf.keys(): -# # e_max = energy_pdf["E Max"] * u.GeV -# # else: -# # e_max = (10 * u.PeV).to(u.GeV) -# -# e_int = energy_pdf.fluence_integral() -# -# flux_1GeV = nu_flu/e_int -# -# print "Flux at 1GeV would be", flux_1GeV, "\n" -# -# source_mc, omega, band_mask = inj.select_mc_band(inj._mc, source) -# -# source_mc["ow"] = flux_1GeV * (inj.mc_weights[band_mask] / omega) -# n_inj = np.sum(source_mc["ow"]) -# -# print "" -# -# print "This corresponds to", n_inj, "neutrinos" -# -# return n_inj - - -# def calculate_neutrinos(source, season, inj_kwargs): -# -# inj = season.make_injector([source], **inj_kwargs) -# energy_pdf = inj_kwargs["Injection Energy PDF"] -# -# print("\n") -# print(source["Name"], season["Name"]) -# print("\n") -# -# if "Flux at 1GeV" not in list(energy_pdf.keys()): -# # if "Source Energy (erg)" in energy_pdf.keys(): -# energy_pdf["Flux at 1GeV"] = \ -# energy_pdf["Source Energy (erg)"] / energy_pdf.fluence_integral() -# -# flux_1_gev = energy_pdf["Flux at 1GeV"] * \ -# inj.time_pdf.effective_injection_time(source) * u.s -# -# print(energy_pdf["Flux at 1GeV"]) -# -# print("Flux at 1GeV would be", flux_1_gev, "\n") -# print( -# "Time is {0} years".format( -# inj.time_pdf.effective_injection_time(source) / (60*60*24*365)) -# ) -# print("Raw Flux is", energy_pdf["Flux at 1GeV"]) -# -# source_mc, omega, band_mask = inj.select_mc_band(inj._mc, source) -# -# # northern_mask = inj._mc["sinDec"] > 0.0 -# -# # print "OneWights:", np.sum(inj._mc["ow"]) -# # print np.sum(inj._mc["ow"] * inj._mc["trueE"]**-energy_pdf["Gamma"]) -# # print np.sum(inj.mc_weights * flux_1_gev) -# # -# # print "Now Ludwig-style:", -# # -# # print np.sum(flux_1_gev * inj.mc_weights[northern_mask]) -# -# source_mc["ow"] = flux_1_gev * inj.mc_weights[band_mask] / omega -# n_inj = np.sum(source_mc["ow"]) -# -# print("") -# -# print("This corresponds to", n_inj, "neutrinos") -# -# return n_inj diff --git a/tests/test_util_astronomy.py b/tests/test_util_astronomy.py index e4a37fa3..a72da018 100644 --- a/tests/test_util_astronomy.py +++ b/tests/test_util_astronomy.py @@ -10,8 +10,6 @@ true_res_astro = { "Energy Flux (GeV cm^{-2} s^{-1})": 1.151292546497023e-08, - "Flux from nearest source": 8.61854306789315e-10 / (u.cm**2 * u.GeV * u.s), - "Mean Luminosity (erg/s)": 9.384215679708581e45, } catalogue = tde_catalogue_name("jetted") @@ -29,23 +27,17 @@ def test_neutrino_astronomy(self): cat = load_catalogue(catalogue) - res_astro = calculate_astronomy(1.0e-9, injection_energy_pdf, cat) + res_astro = calculate_astronomy(1.0e-9, injection_energy_pdf) logging.info("Calculated values {0}".format(res_astro)) logging.info("Reference values {0}".format(true_res_astro)) - for key in ["Energy Flux (GeV cm^{-2} s^{-1})", "Mean Luminosity (erg/s)"]: + for key in ["Energy Flux (GeV cm^{-2} s^{-1})"]: self.assertAlmostEqual( np.log(res_astro[key]), np.log(true_res_astro[key]), places=2 ) - self.assertAlmostEqual( - np.log(res_astro["Flux from nearest source"].value), - np.log(true_res_astro["Flux from nearest source"].value), - places=2, - ) - if __name__ == "__main__": unittest.main()