Commit 2b7d1365a6b41379316f63f16e1d46179fcf3abc
1 parent
09afe671
initial version
Showing
3 changed files
with
537 additions
and
0 deletions
test/example_survey.jpg
0 → 100644
245 KB
test/example_survey.png
0 → 100644
251 KB
test/example_survey.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ===========================================================================================# | ||
| 3 | +# This script simulates CTA surveys. | ||
| 4 | +# | ||
| 5 | +# If matplotlib is installed, the spectrum will be displayed on the screen. | ||
| 6 | +# ===========================================================================================# | ||
| 7 | +from ctatools import * | ||
| 8 | +from gammalib import * | ||
| 9 | +from math import * | ||
| 10 | +import os | ||
| 11 | +import glob | ||
| 12 | +import sys | ||
| 13 | + | ||
| 14 | + | ||
| 15 | +# =========== # | ||
| 16 | +# Plot counts # | ||
| 17 | +# =========== # | ||
| 18 | +def plot_counts(obs): | ||
| 19 | + """ | ||
| 20 | + Plot counts. | ||
| 21 | + """ | ||
| 22 | + # Only proceed if matplotlib is available | ||
| 23 | + if 1 == 1: | ||
| 24 | + #try: | ||
| 25 | + # Import matplotlib | ||
| 26 | + import matplotlib.pyplot as plt | ||
| 27 | + | ||
| 28 | + # Set legend fontsize | ||
| 29 | + params = {'legend.fontsize': 10} | ||
| 30 | + plt.rcParams.update(params) | ||
| 31 | + | ||
| 32 | + # Set plot styles | ||
| 33 | + styles = ['b-', 'g-', 'y-', 'n-'] | ||
| 34 | + | ||
| 35 | + # Dump header | ||
| 36 | + print "" | ||
| 37 | + print "Make plots (using matplotlib):" | ||
| 38 | + print "==============================" | ||
| 39 | + | ||
| 40 | + # Create figure 1 | ||
| 41 | + plt.figure(1,figsize=(12,6)) | ||
| 42 | + plt.subplots_adjust(hspace=.7) | ||
| 43 | + | ||
| 44 | + # Create subplot 1 | ||
| 45 | + plt.subplot(121) | ||
| 46 | + plt.title("Spectrum (summed over all pixels)") | ||
| 47 | + | ||
| 48 | + # Loop over observations | ||
| 49 | + for run in obs: | ||
| 50 | + | ||
| 51 | + # Get event list | ||
| 52 | + list = cast_GCTAEventList(run.events()) | ||
| 53 | + | ||
| 54 | + # Create energy axis | ||
| 55 | + ebounds = GEbounds() | ||
| 56 | + emin = GEnergy() | ||
| 57 | + emax = GEnergy() | ||
| 58 | + emin.TeV(0.1) | ||
| 59 | + emax.TeV(100.0) | ||
| 60 | + ebounds.setlog(emin, emax, 10) | ||
| 61 | + energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())] | ||
| 62 | + | ||
| 63 | + # Create spectrum | ||
| 64 | + print "Extract data:" | ||
| 65 | + counts = [0.0 for i in range(ebounds.size())] | ||
| 66 | + print list.size() | ||
| 67 | + for atom in list: | ||
| 68 | + index = ebounds.index(atom.energy()) | ||
| 69 | + counts[index] = counts[index] + 1.0 | ||
| 70 | + | ||
| 71 | + # Create error bars | ||
| 72 | + error = [sqrt(c) for c in counts] | ||
| 73 | + | ||
| 74 | + # Plot spectrum | ||
| 75 | + plt.loglog(energy, counts, 'ro', label='data') | ||
| 76 | + #plt.errorbar(energy, counts, error, fmt=None, ecolor='r') | ||
| 77 | + | ||
| 78 | + # Extract models | ||
| 79 | + print "Extract models:" | ||
| 80 | + sum_model = [0.0 for i in range(ebounds.size())] | ||
| 81 | + for k, m in enumerate(obs.models()): | ||
| 82 | + print "- "+m.name() | ||
| 83 | + model = [0.0 for i in range(ebounds.size())] | ||
| 84 | + for atom in list: | ||
| 85 | + index = ebounds.index(atom.energy()) | ||
| 86 | + prob = m.eval(atom, run) | ||
| 87 | + model[index] = model[index] + prob * atom.size() | ||
| 88 | + for i in range(ebounds.size()): | ||
| 89 | + sum_model[i] = sum_model[i] + model[i] | ||
| 90 | + #plt.loglog(energy, model, styles[k], label=m.name()) | ||
| 91 | + #plt.loglog(energy, sum_model, 'r-', label='total') | ||
| 92 | + | ||
| 93 | + # Put labels | ||
| 94 | + plt.xlabel("Energy (TeV)") | ||
| 95 | + plt.ylabel("Counts") | ||
| 96 | + plt.legend(loc="lower left") | ||
| 97 | + | ||
| 98 | + # Create subplot 2 | ||
| 99 | + plt.subplot(122) | ||
| 100 | + plt.title("Offset (summed over all energies)") | ||
| 101 | + | ||
| 102 | + # Set Crab direction | ||
| 103 | + crab = GSkyDir() | ||
| 104 | + crab.radec_deg(83.63, 22.01) | ||
| 105 | + | ||
| 106 | + # Loop over observations | ||
| 107 | + for run in obs: | ||
| 108 | + | ||
| 109 | + # Get event cube | ||
| 110 | + list = cast_GCTAEventList(run.events()) | ||
| 111 | + | ||
| 112 | + # Create offset histogram | ||
| 113 | + print "Extract data:" | ||
| 114 | + nx = 30 | ||
| 115 | + doffset2 = 0.01 | ||
| 116 | + offset2 = [(i+0.5)*doffset2 for i in range(nx)] | ||
| 117 | + counts = [0.0 for i in range(nx)] | ||
| 118 | + for atom in list: | ||
| 119 | + off = atom.dir().dist_deg(crab) | ||
| 120 | + off2 = off*off | ||
| 121 | + index = int(off2/doffset2) | ||
| 122 | + if index < nx: | ||
| 123 | + counts[index] = counts[index] + 1.0 | ||
| 124 | + | ||
| 125 | + # Create error bars | ||
| 126 | + error = [sqrt(c) for c in counts] | ||
| 127 | + | ||
| 128 | + # Plot distribution | ||
| 129 | + #plt.semilogy(offset2, counts, 'ro', label='data') | ||
| 130 | + | ||
| 131 | + # Extract models | ||
| 132 | + print "Extract models:" | ||
| 133 | + sum_model = [0.0 for i in range(nx)] | ||
| 134 | + for k, m in enumerate(obs.models()): | ||
| 135 | + print "- "+m.name() | ||
| 136 | + model = [0.0 for i in range(nx)] | ||
| 137 | + for atom in list: | ||
| 138 | + off = atom.dir().dist_deg(crab) | ||
| 139 | + off2 = off*off | ||
| 140 | + index = int(off2/doffset2) | ||
| 141 | + if index < nx: | ||
| 142 | + prob = m.eval(atom, run) | ||
| 143 | + model[index] = model[index] + prob * atom.size() | ||
| 144 | + for i in range(nx): | ||
| 145 | + sum_model[i] = sum_model[i] + model[i] | ||
| 146 | + #plt.plot(offset2, model, styles[k], label=m.name()) | ||
| 147 | + #plt.plot(offset2, sum_model, 'r-', label='total') | ||
| 148 | + #plt.ylim(ymin=0.1) | ||
| 149 | + | ||
| 150 | + # Put labels | ||
| 151 | + plt.xlabel("Offset (deg^2)") | ||
| 152 | + plt.ylabel("Counts") | ||
| 153 | + #plt.legend(loc="upper right") | ||
| 154 | + | ||
| 155 | + # Show counts spectra | ||
| 156 | + plt.show() | ||
| 157 | + | ||
| 158 | + #except: | ||
| 159 | + # print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown." | ||
| 160 | + | ||
| 161 | + # Return | ||
| 162 | + return | ||
| 163 | + | ||
| 164 | + | ||
| 165 | +# ===================== # | ||
| 166 | +# Simulate observations # | ||
| 167 | +# ===================== # | ||
| 168 | +def sim_obs(obs): | ||
| 169 | + """ | ||
| 170 | + Simulate events for all observations in the container. | ||
| 171 | + """ | ||
| 172 | + # Simulate events | ||
| 173 | + sim = ctobssim(obs) | ||
| 174 | + sim.logFileOpen() # We need this to explicitely open the log file in Python mode | ||
| 175 | + sim.run() | ||
| 176 | + | ||
| 177 | + # Extract copy of observations (we need a copy here as sim goes out of | ||
| 178 | + # scope once we leave the function) | ||
| 179 | + obs = sim.obs().copy() | ||
| 180 | + | ||
| 181 | + # Return observations | ||
| 182 | + return obs | ||
| 183 | + | ||
| 184 | + | ||
| 185 | +# ================ # | ||
| 186 | +# Fit observations # | ||
| 187 | +# ================ # | ||
| 188 | +def fit_obs(obs): | ||
| 189 | + """ | ||
| 190 | + Perform maximum likelihood fitting of observations. | ||
| 191 | + """ | ||
| 192 | + # Perform maximum likelihood fitting | ||
| 193 | + like = ctlike(obs) | ||
| 194 | + like.logFileOpen() # We need this to explicitely open the log file in Python mode | ||
| 195 | + like.run() | ||
| 196 | + | ||
| 197 | + # Extract copy of observations (we need a copy here as like goes out of | ||
| 198 | + # scope once we leave the function) | ||
| 199 | + #obs = like.obs().copy() | ||
| 200 | + | ||
| 201 | + #print like.opt() | ||
| 202 | + | ||
| 203 | + # Return observations | ||
| 204 | + return like | ||
| 205 | + | ||
| 206 | + | ||
| 207 | +# ================= # | ||
| 208 | +# Create counts map # | ||
| 209 | +# ================= # | ||
| 210 | +def make_counts_map(obs): | ||
| 211 | + """ | ||
| 212 | + Make counts map by combining the events of all observations. | ||
| 213 | + """ | ||
| 214 | + # Allocate counts map | ||
| 215 | + nxpix = 600 | ||
| 216 | + nypix = 200 | ||
| 217 | + cntmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1) | ||
| 218 | + | ||
| 219 | + # Fill all observations | ||
| 220 | + for run in obs: | ||
| 221 | + | ||
| 222 | + # Loop over all events | ||
| 223 | + for event in run.events(): | ||
| 224 | + | ||
| 225 | + # Cast to CTA event | ||
| 226 | + event = cast_GCTAEventAtom(event) | ||
| 227 | + | ||
| 228 | + # Determine sky pixel | ||
| 229 | + skydir = event.dir().skydir() | ||
| 230 | + pixel = cntmap.dir2pix(skydir) | ||
| 231 | + | ||
| 232 | + # Set pixel | ||
| 233 | + cntmap[pixel] += 1.0 | ||
| 234 | + | ||
| 235 | + # Save sky map | ||
| 236 | + cntmap.save("survey_cntmap.fits", True) | ||
| 237 | + | ||
| 238 | + # Return counts map | ||
| 239 | + return cntmap | ||
| 240 | + | ||
| 241 | + | ||
| 242 | +# ================ # | ||
| 243 | +# Create model map # | ||
| 244 | +# ================ # | ||
| 245 | +def make_model_map(obs): | ||
| 246 | + """ | ||
| 247 | + Make model map by combining all observations. | ||
| 248 | + """ | ||
| 249 | + # Allocate model map | ||
| 250 | + nxpix = 600 | ||
| 251 | + nypix = 200 | ||
| 252 | + modmap = GSkymap("CAR", "GAL", 0.0, 0.0, -0.05, 0.05, nxpix, nypix, 1) | ||
| 253 | + | ||
| 254 | + # Set energy and time | ||
| 255 | + energy = GEnergy() | ||
| 256 | + time = GTime() | ||
| 257 | + instdir = GCTAInstDir() | ||
| 258 | + energy.TeV(0.1) | ||
| 259 | + | ||
| 260 | + # Loop over all map pixels | ||
| 261 | + for pixel in range(modmap.npix()): | ||
| 262 | + | ||
| 263 | + # Get sky direction | ||
| 264 | + skydir = modmap.pix2dir(pixel) | ||
| 265 | + instdir.skydir(skydir) | ||
| 266 | + | ||
| 267 | + # Create event atom for map pixel | ||
| 268 | + atom = GCTAEventAtom() | ||
| 269 | + atom.dir(instdir) | ||
| 270 | + atom.energy(energy) | ||
| 271 | + atom.time(time) | ||
| 272 | + | ||
| 273 | + # Initialise model value | ||
| 274 | + value = 0.0 | ||
| 275 | + | ||
| 276 | + # Loop over all observations | ||
| 277 | + for run in obs: | ||
| 278 | + value += obs.models().eval(atom, run) | ||
| 279 | + | ||
| 280 | + # Set map value | ||
| 281 | + modmap[pixel] = value | ||
| 282 | + | ||
| 283 | + # Save sky map | ||
| 284 | + modmap.save("survey_modmap.fits", True) | ||
| 285 | + | ||
| 286 | + # Return model map | ||
| 287 | + return modmap | ||
| 288 | + | ||
| 289 | + | ||
| 290 | +# ======================= # | ||
| 291 | +# Set one CTA observation # | ||
| 292 | +# ======================= # | ||
| 293 | +def set_one_obs(pntdir, tstart=0.0, duration=1800.0, emin=0.1, emax=100.0, rad=3.0, \ | ||
| 294 | + irf="kb_E_50h_v3", caldb="irf"): | ||
| 295 | + """ | ||
| 296 | + Parameters: | ||
| 297 | + pntdir - Pointing direction | ||
| 298 | + Keywords: | ||
| 299 | + tstart - Start time (seconds) | ||
| 300 | + duration - Duration of observation (seconds) | ||
| 301 | + emin - Minimum event energy (TeV) | ||
| 302 | + emax - Maximum event energy (TeV) | ||
| 303 | + rad - ROI radius (deg) | ||
| 304 | + irf - Instrument response function | ||
| 305 | + caldb - Calibration database path | ||
| 306 | + """ | ||
| 307 | + # Allocate CTA observation | ||
| 308 | + obs = GCTAObservation() | ||
| 309 | + | ||
| 310 | + # Set pointing direction | ||
| 311 | + pnt = GCTAPointing() | ||
| 312 | + pnt.dir(pntdir) | ||
| 313 | + obs.pointing(pnt) | ||
| 314 | + | ||
| 315 | + # Set ROI | ||
| 316 | + roi = GCTARoi() | ||
| 317 | + instdir = GCTAInstDir() | ||
| 318 | + instdir.skydir(pntdir) | ||
| 319 | + roi.centre(instdir) | ||
| 320 | + roi.radius(rad) | ||
| 321 | + | ||
| 322 | + # Set GTI | ||
| 323 | + gti = GGti() | ||
| 324 | + start = GTime() | ||
| 325 | + stop = GTime() | ||
| 326 | + start.met(tstart) | ||
| 327 | + stop.met(tstart+duration) | ||
| 328 | + gti.append(start, stop) | ||
| 329 | + | ||
| 330 | + # Set energy boundaries | ||
| 331 | + ebounds = GEbounds() | ||
| 332 | + e_min = GEnergy() | ||
| 333 | + e_max = GEnergy() | ||
| 334 | + e_min.TeV(emin) | ||
| 335 | + e_max.TeV(emax) | ||
| 336 | + ebounds.append(e_min, e_max) | ||
| 337 | + | ||
| 338 | + # Allocate event list | ||
| 339 | + events = GCTAEventList() | ||
| 340 | + events.roi(roi) | ||
| 341 | + events.gti(gti) | ||
| 342 | + events.ebounds(ebounds) | ||
| 343 | + obs.events(events) | ||
| 344 | + | ||
| 345 | + # Set instrument response | ||
| 346 | + obs.response(irf, caldb) | ||
| 347 | + | ||
| 348 | + # Return observation | ||
| 349 | + return obs | ||
| 350 | + | ||
| 351 | + | ||
| 352 | +# ======================================== # | ||
| 353 | +# Add CTA background model to observations # | ||
| 354 | +# ======================================== # | ||
| 355 | +def add_background_model(obs): | ||
| 356 | + """ | ||
| 357 | + Add standard CTA background model to observations container. | ||
| 358 | + | ||
| 359 | + We use a simple power law here, scaled to Konrad's E configuration | ||
| 360 | + performance table. The model needs still to be validated. | ||
| 361 | + """ | ||
| 362 | + # Recover models from observation | ||
| 363 | + models = obs.models() | ||
| 364 | + | ||
| 365 | + # Define background model | ||
| 366 | + bgd_radial = GCTAModelRadialGauss(3.0) | ||
| 367 | + bgd_spectrum = GModelSpectralPlaw(61.8, -1.85) | ||
| 368 | + bgd_spectrum["Prefactor"].scale(1.0e-6) | ||
| 369 | + bgd_spectrum["PivotEnergy"].value(1.0) | ||
| 370 | + bgd_spectrum["PivotEnergy"].scale(1.0e6) | ||
| 371 | + bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum) | ||
| 372 | + bgd_model.name("Background") | ||
| 373 | + bgd_model.instruments("CTA") | ||
| 374 | + | ||
| 375 | + # Add background model to container | ||
| 376 | + models.append(bgd_model) | ||
| 377 | + | ||
| 378 | + # Put container back in observation container | ||
| 379 | + obs.models(models) | ||
| 380 | + | ||
| 381 | + # Return observation container | ||
| 382 | + return obs | ||
| 383 | + | ||
| 384 | + | ||
| 385 | +# ================= # | ||
| 386 | +# Set Crab spectrum # | ||
| 387 | +# ================= # | ||
| 388 | +def crab_spec(): | ||
| 389 | + """ | ||
| 390 | + Set Crab spectrum based on MAGIC observations | ||
| 391 | + (Albert et al. 2008, ApJ, 674, 1037) | ||
| 392 | + """ | ||
| 393 | + # Set parameters | ||
| 394 | + spectrum = GModelSpectralPlaw(5.7, -2.48) | ||
| 395 | + spectrum["Prefactor"].scale(1.0e-16) | ||
| 396 | + spectrum["PivotEnergy"].value(0.3) | ||
| 397 | + spectrum["PivotEnergy"].scale(1.0e6) | ||
| 398 | + | ||
| 399 | + # Return spectrum | ||
| 400 | + return spectrum | ||
| 401 | + | ||
| 402 | + | ||
| 403 | +# =============================== # | ||
| 404 | +# Setup single observation survey # | ||
| 405 | +# =============================== # | ||
| 406 | +def survey_single(): | ||
| 407 | + """ | ||
| 408 | + Creates a single observation survey for test purposes. | ||
| 409 | + """ | ||
| 410 | + # Allocate observation container | ||
| 411 | + obs = GObservations() | ||
| 412 | + | ||
| 413 | + # Set single pointing at galactic centre | ||
| 414 | + pntdir = GSkyDir() | ||
| 415 | + pntdir.lb_deg(0.0, 0.0) | ||
| 416 | + run = set_one_obs(pntdir) | ||
| 417 | + obs.append(run) | ||
| 418 | + | ||
| 419 | + # Define single point source with Crab flux at galactic centre | ||
| 420 | + center = GSkyDir() | ||
| 421 | + center.lb_deg(0.0, 0.0) | ||
| 422 | + point_spatial = GModelSpatialPtsrc(center) | ||
| 423 | + point_spectrum = crab_spec() | ||
| 424 | + point = GModelPointSource(point_spatial, point_spectrum) | ||
| 425 | + point.name('GC source') | ||
| 426 | + | ||
| 427 | + # Create model container | ||
| 428 | + models = GModels() | ||
| 429 | + models.append(point) | ||
| 430 | + obs.models(models) | ||
| 431 | + | ||
| 432 | + # Return observation container | ||
| 433 | + return obs | ||
| 434 | + | ||
| 435 | + | ||
| 436 | +# =========================== # | ||
| 437 | +# Setup Galactic plane survey # | ||
| 438 | +# =========================== # | ||
| 439 | +def survey_gplane(lrange=10, lstep=2): | ||
| 440 | + """ | ||
| 441 | + Creates a single observation survey for test purposes. | ||
| 442 | + | ||
| 443 | + Keywords: | ||
| 444 | + lrange - Longitude range (integer deg) | ||
| 445 | + lstep - Longitude step size (integer deg) | ||
| 446 | + """ | ||
| 447 | + # Allocate observation container | ||
| 448 | + obs = GObservations() | ||
| 449 | + | ||
| 450 | + # Loop over longitudes | ||
| 451 | + for l in range(-lrange,lrange+lstep,lstep): | ||
| 452 | + | ||
| 453 | + # Set pointing | ||
| 454 | + pntdir = GSkyDir() | ||
| 455 | + pntdir.lb_deg(l, 0.0) | ||
| 456 | + run = set_one_obs(pntdir) | ||
| 457 | + obs.append(run) | ||
| 458 | + | ||
| 459 | + # Define single point source with Crab flux at galactic centre | ||
| 460 | + center = GSkyDir() | ||
| 461 | + center.lb_deg(0.0, 0.0) | ||
| 462 | + point_spatial = GModelSpatialPtsrc(center) | ||
| 463 | + point_spectrum = crab_spec() | ||
| 464 | + point = GModelPointSource(point_spatial, point_spectrum) | ||
| 465 | + point.name('GC source') | ||
| 466 | + | ||
| 467 | + # Create model container | ||
| 468 | + models = GModels() | ||
| 469 | + models.append(point) | ||
| 470 | + obs.models(models) | ||
| 471 | + | ||
| 472 | + # Return observation container | ||
| 473 | + return obs | ||
| 474 | + | ||
| 475 | + | ||
| 476 | +#==========================# | ||
| 477 | +# Main routine entry point # | ||
| 478 | +#==========================# | ||
| 479 | +if __name__ == '__main__': | ||
| 480 | + """ | ||
| 481 | + CTA survey simulation and analysis script. | ||
| 482 | + """ | ||
| 483 | + # Initialise flags | ||
| 484 | + need_help = False | ||
| 485 | + | ||
| 486 | + # Test for command line arguments | ||
| 487 | + print sys.argv[0] | ||
| 488 | + if (len(sys.argv) > 1): | ||
| 489 | + if sys.argv[1] == "-h": | ||
| 490 | + need_help = True | ||
| 491 | + else: | ||
| 492 | + need_help = True | ||
| 493 | + | ||
| 494 | + # Print help if needed and exit | ||
| 495 | + if need_help: | ||
| 496 | + print "Usage: example_survey.py [OPTIONS]" | ||
| 497 | + print " -h Display this usage message" | ||
| 498 | + sys.exit() | ||
| 499 | + | ||
| 500 | + # Dump header | ||
| 501 | + print "**********************" | ||
| 502 | + print "* Simulate CTA survey *" | ||
| 503 | + print "**********************" | ||
| 504 | + | ||
| 505 | + # Remove any existing result files | ||
| 506 | + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] | ||
| 507 | + for files in list: | ||
| 508 | + for file in files: | ||
| 509 | + os.remove(file) | ||
| 510 | + | ||
| 511 | + # Setup single observation survey | ||
| 512 | + obs = survey_single() | ||
| 513 | + obs = survey_gplane() | ||
| 514 | + | ||
| 515 | + # Add background model | ||
| 516 | + obs = add_background_model(obs) | ||
| 517 | + | ||
| 518 | + # Simulate events | ||
| 519 | + print "Simulate events" | ||
| 520 | + obs = sim_obs(obs) | ||
| 521 | + | ||
| 522 | + # Make counts map | ||
| 523 | + print "Make counts map" | ||
| 524 | + cntmap = make_counts_map(obs) | ||
| 525 | + | ||
| 526 | + # Fit observations | ||
| 527 | + print "Fit observatins" | ||
| 528 | + like = fit_obs(obs) | ||
| 529 | + #print like.opt() | ||
| 530 | + #print like.obs().models() | ||
| 531 | + | ||
| 532 | + # Make model map | ||
| 533 | + print "Make model map (this step will take some time)" | ||
| 534 | + modmap = make_model_map(obs) | ||
| 535 | + | ||
| 536 | + # Show fit results | ||
| 537 | + #plot_counts(like.obs()) |