Skip to content

mapmaker_README.md

Latest commit

 

History

History
196 lines (173 loc) · 9.52 KB

mapmaker_README.md

File metadata and controls

196 lines (173 loc) · 9.52 KB

Code to calculate the probability of seeing a EM counterpart to a LIGO event starting from the LIGO probability maps.

At present, calculates the limiting magnitude of a 30 second i-band Blanco image, all sky.

Try: ./setup.sh

import matplotlib.pyplot as plt; plt.ion() ; figure = plt.figure(1) import numpy as np; import os import mags

get ligo map

import hp2np import healpy as hp skymap = "lalinference.fits.gz"; resolution = 128 ra,dec,ligo = hp2np.hp2np(skymap) j,j,ligo_dist = hp2np.hp2np(skymap, degrade=resolution, field=1) j,j,ligo_dist_var =hp2np.hp2np(skymap, degrade=resolution, field=2) j,j,ligo_dist_norm = hp2np.hp2np(skymap, degrade=resolution, field=3)

compute limiting mag

mjd = 55435.60552; obs = mags.observed(ra,dec,ligo, mjd); omags = obs.limitMag("i")

since this turns out to be during the day, change the time.

obs.resetTime(mjd+(10./24.)); obs.limitMag("i")

perhaps make the exposures longer than 30 sec, say 90 seconds

obs.limitMag("i",exposure=90.)

plot as ra,dec map (all coordinates are stored as radians)

plt.clf(); plt.hexbin( obs.ra360./2/np.pi, obs.dec360./2/np.pi, obs.maglim, vmin=15); plt.colorbar()

plot in an equal area projection- x,y are projected into a mcbryde-thomas projection

plt.axes().set_aspect('equal') plt.clf(); plt.hexbin( obs.x, obs.y, obs.maglim, vmin=15); plt.colorbar()

work with the source probabilty

import sourceProb sm=sourceProb.map(obs,type="NS"); sm.calculateProb(ligo, ligo_dist, ligo_dist_sigma) plt.clf();plt.hexbin(obs.x,obs.y,sm.probMap); plt.colorbar() plt.clf();plt.hexbin(obs.x,obs.y,sm.probMap*ligo); plt.colorbar()

show the ligo probability contours on the limiting magnitude map

plt.hexbin( obs.x, obs.y, obs.maglim, vmin=15); plt.colorbar() sm.plotLigoContours(obs)

show the ligo probability contours on the source detection probability map

plt.hexbin( obs.x, obs.y, sm.probMap); plt.colorbar() sm.plotLigoContours(obs)

show the ligo*source probability contours on the limiting magnitude map

plt.hexbin( obs.x, obs.y, obs.maglim, vmin=17); plt.colorbar() sm.plotLigoContours(obs, type="ls")

do it better

from equalArea import mcplot plt.clf(); mcplot.plot(obs.x,obs.y,obs.maglim,300,vmin=17); sm.plotLigoContours(obs,type="ls") plt.savefig("thefig.pdf", bbox_inches="tight")

get the limiting distance map, assuming a fixed absolute magnitude

sm=sourceProb.map(obs, type="NS"); sm.calculateProb() plt.clf(); plt.hexbin( obs.x, obs.y, sm.prob); plt.colorbar()

lower the resolution of the maps for speed

ra,dec,ligo = hp2np.hp2np(os.environ['DESGW_DATA_DIR'] + "/bayestar.fits.gz") ra2,dec2,ligo2 =hp2np.map2np(ligo,resolution=32) obs = mags.observed(ra2,dec2,ligo2 mjd, degradeRes=True); obs.resetTime(mjd+(10./24.)); obs.limitMag("i", exposure=180); sm=sourceProb.map(obs, type="NS");sm.calculateProb();

plot using hour angle instead of RA- area always in center

mcplot.plot(obs.hx,obs.hy,sm.probMap,cmap="gray"); sm.plotLigoContours(obs,type="ligo",hourangle=True);

More to come. Jim Annis Dec 5, 2014

=====================================================

Now let us industrialize this.

First for a simulation, sim 13681

import getHexObservations; import mapsAtTimeT; import allMaps; import modelRead

=== we have to get the simulation information

simNum = 13681 data_dir = "/data/des30.a/data/annis/des-gw/ligo/sims/" sims, mjds, distances, ligoMap = allMaps.selectSingleSim(simNum, data_dir)

=== let's reset the time of the burst to make it useable

mjds[0]=57295. - 0.5

=== now we are ready to do our standard work

ligo = hp.read_map(ligoMap) ra,dec,ligo = hp2np.map2np(ligo,256, fluxConservation=True) obs = mags.observed(ra,dec,ligo, mjds[0]) obs.limitMag("i",exposure=180.)

==== get the neutron star explosion models

models = modelRead.getModels()

==== calculate maps during a full night of observing

probs,times = mapsAtTimeT.oneDayOfTotalProbability(obs,mjds[0],distances[0],models) mapsAtTimeT.probabilityMapSaver (obs, sims[0], mjds[0], distances[0], models, times, probs, "jack/")

==== figure out what to observe during that night

n_slots, map_zero = getHexObservations.contemplateTheDivisionsOfTime( probs, times, hoursAvailable=6) best_slot = getHexObservations.now( n_slots,
mapDirectory="jack/", simNumber=13681, mapZero=map_zero)

==== show what we learned

plt.clf();getHexObservations.observingPlot(figure,13681,5,"jack/",12)

recall that "figure" is from "figure = plt.figure(1)"

And let's try an example that works for Dylan Brout's MainInjector

have a devel version of maininjector around, say at:

/data/des41.a/data/annis/maininjector

=== the recycler object knows about skymap, mjd, trigger_id, and outfolder

=== it has exposure_length as a config variable and distance from the event

import numpy as np; import os; import getHexObservations out_dir = "jack/"; data_dir = out_dir web_dir = "jack/web/" skymap = "/data/des41.a/data/desgw/maininjector/test-triggers/M168556/M168556_bayestar.fits.gz" trigger_id = "/M168556" ; trigger_mjd = 55472.8245; distance = 91.24 hours_available=6.; exposure_length=180.;

==== calculate maps during a full night of observing

probs,times = getHexObservations.prepare(skymap, trigger_mjd, trigger_id,
data_dir, exposure_length=exposure_length, distance=distance) n_slots, map_zero = getHexObservations.contemplateTheDivisionsOfTime( probs, times, hoursAvailable=hours_available);

==== figure out what to observe during that night

best_slot = getHexObservations.now(
n_slots, mapDirectory=data_dir, simNumber=trigger_id, mapZero=map_zero) getHexObservations.makeObservingPlots( n_slots, trigger_id, best_slot, data_dir)

======== debug hexalation

obs, trigger_id, mjd, distance, models, times, probs,data_dir =
getHexObservations.prepare(skymap, trigger_mjd, trigger_id,
data_dir, exposure_length=exposure_length, distance=distance) gw_data_dir = os.environ["DESGW_DATA_DIR"];hexFile = gw_data_dir + "all-sky-hexCenters.txt" print probs,times obs,sm = mapsAtTimeT.probabilityMaps( obs, mjd, times[10], distance, models) raHexen, decHexen, hexVals, rank = hexalate.cutAndHexalate (
obs, sm, hexFile)

======== real event

import numpy as np; import os; import getHexObservations out_dir = "jack2/"; data_dir = out_dir skymap = "/data/des41.a/data/desgw/maininjector_devel/real-triggers/G184098/G184098_bayestar.fits" skymap = "/data/des41.a/data/desgw/maininjector_devel/real-triggers/G184098/skyprobcc_cWB_complete.fits" trigger_id = "G184098" ; trigger_mjd = 57279.4102; distance = 10. hours_available=3.; exposure_length=180.;

==== calculate maps during a full night of observing

probs,times = getHexObservations.prepare(skymap, trigger_mjd, trigger_id,
data_dir, exposure_length=exposure_length, distance=distance) n_slots, map_zero = getHexObservations.contemplateTheDivisionsOfTime( probs, times, hoursAvailable=hours_available);

==== figure out what to observe during that night

best_slot = getHexObservations.now(
n_slots, mapDirectory=data_dir, simNumber=trigger_id, mapZero=map_zero, skipJson=True) getHexObservations.makeObservingPlots( n_slots, trigger_id, best_slot, data_dir)

test out the new start_mjd and the rationalization of slot duration etc

30 minute per slot, 5 hexes/slot

thursday mid night 57283 -0.5 call it 57282.6

probs,times = getHexObservations.prepare(skymap, trigger_mjd, trigger_id,
data_dir, exposure_length=exposure_length, distance=distance,
start_mjd=57285.0, deltaTime =0.02083) n_slots, map_zero = getHexObservations.contemplateTheDivisionsOfTime( probs, times, hoursAvailable=hours_available, slotDuration=30.); best_slot = getHexObservations.now(
n_slots, mapDirectory=data_dir, simNumber=trigger_id,
mapZero=map_zero, maxHexesPerSlot=5, skipJson=True) getHexObservations.makeObservingPlots( n_slots, trigger_id, best_slot, data_dir) reload(getHexObservations); reload(mapsAtTimeT)

rahex,dechex = np.genfromtxt(hexFile, unpack=True) ra2, dec2 = np.genfromtxt("python/jack2/G184098-ra-dec-prob-mjd-slot.txt", unpack=True) ra_lmc, dec_lmc = np.genfromtxt("python/jack2/lmc.radec.txt", unpack=True) tra, tdec, tval = getHexObservations.tester(ra2, dec2, rahex, dechex, hexVals) print " in 3 slots", tval.sum() tra, tdec, tval = getHexObservations.tester(ra_lmc, dec_lmc, ra, dec, vals) print " in lmc ", tval.sum() ix = (ra2< 140); print ra2.size, ra2[ix].size tra, tdec, tval = getHexObservations.tester(ra2[ix], dec2[ix], ra, dec, vals) print " in 2.4 slots", tval.sum()

ras, decs = np.genfromtxt("jack2/lasilla.txt", unpack=True, usecols=(1,2))

night 3 : 57285.0 data_dir="night3/"; start_mjd=57285.0 night 4 : 57286.0 data_dir="night4/"; start_mjd=57286.0 probs,times = getHexObservations.prepare(skymap, trigger_mjd, trigger_id,data_dir, exposure_length=exposure_length, distance=distance,start_mjd=57286.0, deltaTime =0.02083,onlyHexesAlreadyDone="../records/G184098-ra-dec-prob-mjd-slot.txt") n_slots, map_zero = getHexObservations.contemplateTheDivisionsOfTime(probs, times, hoursAvailable=hours_available, slotDuration=30.); best_slot = getHexObservations.now(n_slots, mapDirectory=data_dir, simNumber=trigger_id,mapZero=map_zero, maxHexesPerSlot=5, skipJson=False)

More to come. Jim Annis Sept 11, 2015