updated skymaps

include/Figures.h

@@ -20,8 +20,8 @@ bool E_sorter(ICEvent event_i, ICEvent event_j);
 void Likelihood_Table();
 
 // visualizations
-void SkyMap();
-void MW_Visualization();
+void IC_SkyMap();
+void MW_SkyMap();
 
 // testing
 void vMF_test();

py/SkyMap.py → py/IC_SkyMap.py

@@ -14,7 +14,7 @@
 
 print "Reading data..."
 # read in data
-dataf = open("data/SkyMap.txt", "r")
+dataf = open("data/IC_SkyMap.txt", "r")
 N_thetas, N_phis, N_Repeat = [int(x) for x in dataf.readline().split()]
 dts = ["count"]
 dt = [(d, "f") for d in dts]
@@ -67,42 +67,38 @@
 
 print "Drawing skymap lines..."
 # draw map with lines
-m = Basemap(projection = "hammer", lat_0 = 0, lon_0 = 0)
 m = Basemap(projection = "moll", lat_0 = 0, lon_0 = 0)
 m.drawmapboundary(fill_color = "none") # draw a line around the map region
 m.drawparallels(np.arange(-90., 120., 30.), labels = [1, 0, 0, 0], color = "gray", linewidth = 0.6, fontsize = 12) # draw parallels
 m.drawmeridians([0, 60, 120, 240, 300], fontsize = 12, color = "gray", linewidth = 0.6) # draw meridians
 
-lats = 90 - 180 * np.indices((N_thetas, N_phis))[0, :, :] / (N_thetas - 1)
-lons = 360 * np.indices((N_thetas, N_phis))[1, :, :] / (N_thetas - 1)
+lats = 90 - 180 * (np.indices((N_thetas, N_phis))[0, :, :] + 0.5) / N_thetas
+lons = 360 * (np.indices((N_thetas, N_phis))[1, :, :] + 0.5) / N_thetas
 lons = 180 - lons # reverse to agree with IC plots
 x, y = m(lons, lats)
 
 print "Drawing contours..."
 levels = np.linspace(dmin - 0.4, dmax, 100)
-print dmin, dmax
 c = m.contourf(x, y, data, levels = levels, cmap = plt.get_cmap("Oranges"), zorder = -9)
 plt.gca().set_rasterization_zorder(-1)
 print ""
-#print "Fixing contour edges..."
-#for collection in c.collections:
-#	collection.set_edgecolor("face")
 
 print "Drawing text..."
+ta = plt.gca().transAxes
 if gal:
-	plt.text(0.98, 0.02, r"${\rm Galactic}$", transform = plt.gca().transAxes, ha = "right", va = "bottom", fontsize = 12)
-	plt.text(0.24, 0.496, r"${\rm Galactic\ Plane}$", transform = plt.gca().transAxes, ha = "center", va = "top", fontsize = 7)
+	plt.text(0.98, 0.02, r"${\rm Galactic}$", transform = ta, ha = "right", va = "bottom", fontsize = 12)
+	plt.text(0.24, 0.496, r"${\rm Galactic\ Plane}$", transform = ta, ha = "center", va = "top", fontsize = 7)
+	plt.text(0.01, 0.52, r"$+180^\circ$", transform = ta, ha = "left", va = "bottom", fontsize = 7)
+	plt.text(0.99, 0.52, r"$-180^\circ$", transform = ta, ha = "right", va = "bottom", fontsize = 7)
 else:
-	plt.text(0.98, 0.02, r"${\rm Equatorial}$", transform = plt.gca().transAxes, ha = "right", va = "bottom", fontsize = 12)
-	plt.text(0.92, 0.496, r"${\rm Equatorial\ Plane}$", transform = plt.gca().transAxes, ha = "center", va = "top", fontsize = 7)
+	plt.text(0.98, 0.02, r"${\rm Equatorial}$", transform = ta, ha = "right", va = "bottom", fontsize = 12)
+	plt.text(0.92, 0.496, r"${\rm Equatorial\ Plane}$", transform = ta, ha = "center", va = "top", fontsize = 7)
 
 for i in xrange(len(event_ids)):
 	x, y = m(ls[i], bs[i])
 	plt.text(x, y, r"$%i$" % event_ids[i], ha = "center", va = "center", fontsize = 10, color = (0, 0, 0.6))
 
-dpi = 400
+dpi = 250
 print "Saving .eps..."
-plt.savefig("fig/SkyMap.eps", bbox_inches = "tight", dpi = dpi)
-#print "Saving .pdf..."
-#plt.savefig("fig/SkyMap.pdf", bbox_inches = "tight", dpi = dpi)
+plt.savefig("fig/IC_SkyMap.eps", bbox_inches = "tight", dpi = dpi)
 

py/MW_SkyMap.py

@@ -0,0 +1,64 @@
+"""
+This code is free to use, copy, distribute, and modify.
+If you use this code or any modification of this code, we request that you reference both this code https://zenodo.org/record/x and the paper https://arxiv.org/abs/17xx.xxxxx.
+"""
+
+print "Importing..."
+from mpl_toolkits.basemap import Basemap
+import matplotlib.pyplot as plt
+import numpy as np
+
+print "Reading data..."
+# read in data
+dataf = open("data/MW_SkyMap.txt", "r")
+N_thetas, N_phis, N_Repeat = [int(x) for x in dataf.readline().split()]
+dts = ["count"]
+dt = [(d, "f") for d in dts]
+data = np.loadtxt(dataf, dtype = dt)
+dataf.close()
+data = np.array(data["count"])
+
+# get ranges
+dmax = np.nanmax(data)
+
+for i in xrange(len(data)):
+	if np.isinf(data[i]):
+		data[i] = dmax + 1
+
+dmin = np.nanmin(data)
+
+for i in xrange(len(data)):
+	if data[i] == dmax + 1:
+		data[i] = dmin - 0.3
+
+# reshape
+data.resize(N_thetas, N_phis)
+
+print "Drawing skymap lines..."
+# draw map with lines
+m = Basemap(projection = "moll", lat_0 = 0, lon_0 = 0)
+m.drawmapboundary(fill_color = "none") # draw a line around the map region
+m.drawparallels(np.arange(-90., 120., 30.), labels = [1, 0, 0, 0], color = "gray", linewidth = 0.6, fontsize = 12) # draw parallels
+m.drawmeridians([0, 60, 120, 240, 300], fontsize = 12, color = "gray", linewidth = 0.6) # draw meridians
+
+lats = 90 - 180 * (np.indices((N_thetas, N_phis))[0, :, :] + 0.5) / N_thetas
+lons = 360 * (np.indices((N_thetas, N_phis))[1, :, :] + 0.5) / N_thetas
+lons = 180 - lons # reverse to agree with IC plots
+x, y = m(lons, lats)
+
+print "Drawing contours..."
+levels = np.linspace(dmin - 0.4, dmax, 100)
+c = m.contourf(x, y, data, levels = levels, cmap = plt.get_cmap("Blues"), zorder = -9)
+plt.gca().set_rasterization_zorder(-1)
+print ""
+
+print "Drawing text..."
+ta = plt.gca().transAxes
+plt.text(0.98, 0.02, r"${\rm Galactic}$", transform = ta, ha = "right", va = "bottom", fontsize = 12)
+plt.text(0.01, 0.52, r"$+180^\circ$", transform = ta, ha = "left", va = "bottom", fontsize = 7)
+plt.text(0.99, 0.52, r"$-180^\circ$", transform = ta, ha = "right", va = "bottom", fontsize = 7)
+
+dpi = 250
+print "Saving .eps..."
+plt.savefig("fig/MW_SkyMap.eps", bbox_inches = "tight", dpi = dpi)
+

src/Figures.cpp

@@ -113,11 +113,11 @@ void Likelihood_Table()
 	data.close();
 }
 
-void SkyMap()
+void IC_SkyMap()
 {
-	std::cout << "Generating SkyMap..." << std::endl;
+	std::cout << "Generating IC SkyMap..." << std::endl;
 
-	std::ofstream data("data/SkyMap.txt");
+	std::ofstream data("data/IC_SkyMap.txt");
 	coord2D coord_gal, coord_gal_smeared;
 
 	int k, l; // indices for the grid
@@ -129,6 +129,7 @@ void SkyMap()
 	data << N_thetas << " " << N_phis << " " << N_Repeat << std::endl;
 
 	int *grid = new int[N_thetas * N_phis];
+	for (int i = 0; i < N_thetas * N_phis; i++) grid[i] = 0; // initialize the grid to zero
 
 	Progress_Bar *pbar = new Progress_Bar();
 	pbar->update(0);
@@ -139,8 +140,8 @@ void SkyMap()
 		for (int j = 0; j < N_Repeat; j++)
 		{
 			coord_gal_smeared = vMF_smear(coord_gal, events[i].sigma_direction);
-			k = (N_thetas - 1) * fmod(coord_gal_smeared.theta, M_PI) / M_PI;
-			l = (N_phis - 1) * fmod(11 * M_PI + coord_gal_smeared.phi, 2 * M_PI) / (2 * M_PI);
+			k = N_thetas * fmod(coord_gal_smeared.theta, M_PI) / M_PI;
+			l = N_phis * fmod(11 * M_PI + coord_gal_smeared.phi, 2 * M_PI) / (2 * M_PI);
 			grid[k * N_phis + l]++;
 
 			if (j % 10000 == 0)
@@ -153,11 +154,8 @@ void SkyMap()
 	for (int i = 0; i < N_thetas * N_phis; i++)
 	{
 		k = int(1.0 * i / (N_phis));
-		theta = M_PI * k / (N_phis - 1);
-		if (k == 0 or k == N_phis - 1)
-			data << log(grid[i] / N_Repeat) << std::endl;
-		else
-			data << log(grid[i] / sin(theta) / N_Repeat) << std::endl;
+		theta = M_PI * (k + 0.5) / N_phis;
+		data << log(grid[i] / sin(theta) / N_Repeat) << std::endl;
 	}
 	data.close();
 
@@ -166,9 +164,56 @@ void SkyMap()
 	std::cout << "Done." << std::endl;
 }
 
-void MW_Visualization()
+void MW_SkyMap()
 {
-	std::ofstream data("data/MW_Visualization.txt");
+	std::cout << "Generating MW SkyMap..." << std::endl;
+
+	std::ofstream data("data/MW_SkyMap.txt");
+	coord2D coord_gal;
+
+	int k, l; // indices for the grid
+	const int N_thetas = 5e2;
+	const int N_phis = 5e2;
+	int N_Repeat = 1e9;
+	double theta;
+
+	data << N_thetas << " " << N_phis << " " << N_Repeat << std::endl;
+
+	int *grid = new int[N_thetas * N_phis];
+	for (int i = 0; i < N_thetas * N_phis; i++) grid[i] = 0; // initialize the grid to zero
+
+	Progress_Bar *pbar = new Progress_Bar();
+	pbar->update(0);
+	// populates the grid
+	for (int i = 0; i < N_Repeat; i++)
+	{
+		coord_gal = MW();
+		k = N_thetas * (fmod(coord_gal.theta, M_PI) / M_PI);
+		l = N_phis * (fmod(11 * M_PI + coord_gal.phi, 2 * M_PI) / (2 * M_PI));
+		grid[k * N_phis + l]++;
+
+		if (i % 10000 == 0)
+			pbar->update(0, N_Repeat, i, true);
+	} // i, N_Repeat
+	delete pbar;
+
+	// writes the grid to file
+	for (int i = 0; i < N_thetas * N_phis; i++)
+	{
+		k = int(1.0 * i / (N_phis));
+		theta = M_PI * (k + 0.5) / N_phis;
+		data << log(grid[i] / sin(theta) / N_Repeat) << std::endl;
+	}
+	data.close();
+
+	delete[] grid;
+
+	std::cout << "Done." << std::endl;
+}
+
+void MW_SkyMap_old()
+{
+	std::ofstream data("data/MW_SkyMap.txt");
 	int N_Repeat;
 	coord_cart coord_c;
 	N_Repeat = 1e7;

src/main.cpp

@@ -24,8 +24,8 @@ int main()
 	// generate the data for figures and tables
 	Likelihood();
 	Likelihood_Table();
-//	SkyMap();
-//	MW_Visualization();
+//	IC_SkyMap();
+//	MW_SkyMap();
 	vMF_test();
 
 	return 0;