Commit dd4efb4666383c0dbbb51fa371c30c18dcf595c4
1 parent
d8a27288
move to examples directory
Showing
2 changed files
with
0 additions
and
845 deletions
test/example_binned.py deleted
| 1 | -#! /usr/bin/env python | ||
| 2 | -# ========================================================================== | ||
| 3 | -# This script illustrates how to build a binned analysis pipeline using the | ||
| 4 | -# ctatools. If matplotlib is installed, a counts spectrum and an offset | ||
| 5 | -# histogram will be displayed. | ||
| 6 | -# | ||
| 7 | -# Copyright (C) 2011 Jurgen Knodlseder | ||
| 8 | -# | ||
| 9 | -# This program is free software: you can redistribute it and/or modify | ||
| 10 | -# it under the terms of the GNU General Public License as published by | ||
| 11 | -# the Free Software Foundation, either version 3 of the License, or | ||
| 12 | -# (at your option) any later version. | ||
| 13 | -# | ||
| 14 | -# This program is distributed in the hope that it will be useful, | ||
| 15 | -# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 16 | -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 17 | -# GNU General Public License for more details. | ||
| 18 | -# | ||
| 19 | -# You should have received a copy of the GNU General Public License | ||
| 20 | -# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 21 | -# | ||
| 22 | -# ========================================================================== | ||
| 23 | -from ctatools import * | ||
| 24 | -from gammalib import * | ||
| 25 | -from math import * | ||
| 26 | -import os | ||
| 27 | -import glob | ||
| 28 | - | ||
| 29 | - | ||
| 30 | -# ============================= # | ||
| 31 | -# Analysis pipeline - version 1 # | ||
| 32 | -# ============================= # | ||
| 33 | -def pipeline_v1(): | ||
| 34 | - """ | ||
| 35 | - Binned analysis pipeline - save intermediate results. | ||
| 36 | - | ||
| 37 | - This function implements an analysis pipeline that successively calls | ||
| 38 | - ctobssim, ctbin and ctlike by saving the intermediate results as FITS | ||
| 39 | - files on disk. This replicates the way how the analysis would be done | ||
| 40 | - in a ftools-like approach. | ||
| 41 | - """ | ||
| 42 | - # Set script parameters | ||
| 43 | - model_name = "data/crab.xml" | ||
| 44 | - events_name = "events.fits" | ||
| 45 | - cntmap_name = "cntmap.fits" | ||
| 46 | - result_name = "results.xml" | ||
| 47 | - caldb = "irf" | ||
| 48 | - irf = "kb_E_50h_v3" | ||
| 49 | - ra = 83.63 | ||
| 50 | - dec = 22.01 | ||
| 51 | - rad_sim = 10.0 | ||
| 52 | - tstart = 0.0 | ||
| 53 | - tstop = 1800.0 | ||
| 54 | - emin = 0.1 | ||
| 55 | - emax = 100.0 | ||
| 56 | - enumbins = 20 | ||
| 57 | - nxpix = 200 | ||
| 58 | - nypix = 200 | ||
| 59 | - binsz = 0.02 | ||
| 60 | - coordsys = "CEL" | ||
| 61 | - proj = "CAR" | ||
| 62 | - | ||
| 63 | - # Initialise timing | ||
| 64 | - wall_seconds = 0.0 | ||
| 65 | - cpu_seconds = 0.0 | ||
| 66 | - | ||
| 67 | - # Simulate events | ||
| 68 | - sim = ctobssim() | ||
| 69 | - sim.logFileOpen() # We need this to explicitely open the log file in Python mode | ||
| 70 | - sim["infile"].filename(model_name) | ||
| 71 | - sim["outfile"].filename(events_name) | ||
| 72 | - sim["caldb"].string(caldb) | ||
| 73 | - sim["irf"].string(irf) | ||
| 74 | - sim["ra"].real(ra) | ||
| 75 | - sim["dec"].real(dec) | ||
| 76 | - sim["rad"].real(rad_sim) | ||
| 77 | - sim["tmin"].real(tstart) | ||
| 78 | - sim["tmax"].real(tstop) | ||
| 79 | - sim["emin"].real(emin) | ||
| 80 | - sim["emax"].real(emax) | ||
| 81 | - sim.execute() | ||
| 82 | - print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 83 | - | ||
| 84 | - # Update timing | ||
| 85 | - wall_seconds += sim.telapse() | ||
| 86 | - cpu_seconds += sim.celapse() | ||
| 87 | - | ||
| 88 | - # Bin events into counts map | ||
| 89 | - bin = ctbin() | ||
| 90 | - bin.logFileOpen() # We need this to explicitely open the log file in Python mode | ||
| 91 | - bin["evfile"].filename(events_name) | ||
| 92 | - bin["outfile"].filename(cntmap_name) | ||
| 93 | - bin["emin"].real(emin) | ||
| 94 | - bin["emax"].real(emax) | ||
| 95 | - bin["enumbins"].integer(enumbins) | ||
| 96 | - bin["nxpix"].integer(nxpix) | ||
| 97 | - bin["nypix"].integer(nypix) | ||
| 98 | - bin["binsz"].real(binsz) | ||
| 99 | - bin["coordsys"].string(coordsys) | ||
| 100 | - bin["xref"].real(ra) | ||
| 101 | - bin["yref"].real(dec) | ||
| 102 | - bin["proj"].string(proj) | ||
| 103 | - bin.execute() | ||
| 104 | - print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)" | ||
| 105 | - | ||
| 106 | - # Update timing | ||
| 107 | - wall_seconds += bin.telapse() | ||
| 108 | - cpu_seconds += bin.celapse() | ||
| 109 | - | ||
| 110 | - # Perform maximum likelihood fitting | ||
| 111 | - like = ctlike() | ||
| 112 | - like.logFileOpen() # We need this to explicitely open the log file in Python mode | ||
| 113 | - like["cntmap"].filename(cntmap_name) | ||
| 114 | - like["srcmdl"].filename(model_name) | ||
| 115 | - like["outmdl"].filename(result_name) | ||
| 116 | - like["method"].string("BINNED") | ||
| 117 | - like["caldb"].string(caldb) | ||
| 118 | - like["irf"].string(irf) | ||
| 119 | - like.execute() | ||
| 120 | - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 121 | - | ||
| 122 | - # Update timing | ||
| 123 | - wall_seconds += like.telapse() | ||
| 124 | - cpu_seconds += like.celapse() | ||
| 125 | - | ||
| 126 | - # Show total times | ||
| 127 | - print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 128 | - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 129 | - | ||
| 130 | - # Return | ||
| 131 | - return | ||
| 132 | - | ||
| 133 | - | ||
| 134 | -# ============================= # | ||
| 135 | -# Analysis pipeline - version 2 # | ||
| 136 | -# ============================= # | ||
| 137 | -def pipeline_v2(): | ||
| 138 | - """ | ||
| 139 | - Binned analysis pipeline - keep intermediate results in memory. | ||
| 140 | - | ||
| 141 | - This function implements an analysis pipeline that successively calls | ||
| 142 | - ctobssim, ctbin and ctlike without saving the intermediate results as | ||
| 143 | - FITS files on disk. All data are hold in memory. | ||
| 144 | - | ||
| 145 | - At the end, results are plotted (if matplotlib is installed) | ||
| 146 | - """ | ||
| 147 | - # Set script parameters | ||
| 148 | - model_name = "data/crab.xml" | ||
| 149 | - caldb = "irf" | ||
| 150 | - irf = "kb_E_50h_v3" | ||
| 151 | - ra = 83.63 | ||
| 152 | - dec = 22.01 | ||
| 153 | - rad_sim = 10.0 | ||
| 154 | - tstart = 0.0 | ||
| 155 | - tstop = 1800.0 | ||
| 156 | - emin = 0.1 | ||
| 157 | - emax = 100.0 | ||
| 158 | - enumbins = 20 | ||
| 159 | - nxpix = 200 | ||
| 160 | - nypix = 200 | ||
| 161 | - binsz = 0.02 | ||
| 162 | - coordsys = "CEL" | ||
| 163 | - proj = "CAR" | ||
| 164 | - | ||
| 165 | - # Initialise timing | ||
| 166 | - wall_seconds = 0.0 | ||
| 167 | - cpu_seconds = 0.0 | ||
| 168 | - | ||
| 169 | - # Simulate events | ||
| 170 | - sim = ctobssim() | ||
| 171 | - sim["infile"].filename(model_name) | ||
| 172 | - sim["caldb"].string(caldb) | ||
| 173 | - sim["irf"].string(irf) | ||
| 174 | - sim["ra"].real(ra) | ||
| 175 | - sim["dec"].real(dec) | ||
| 176 | - sim["rad"].real(rad_sim) | ||
| 177 | - sim["tmin"].real(tstart) | ||
| 178 | - sim["tmax"].real(tstop) | ||
| 179 | - sim["emin"].real(emin) | ||
| 180 | - sim["emax"].real(emax) | ||
| 181 | - sim.run() | ||
| 182 | - print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 183 | - | ||
| 184 | - # Update timing | ||
| 185 | - wall_seconds += sim.telapse() | ||
| 186 | - cpu_seconds += sim.celapse() | ||
| 187 | - | ||
| 188 | - # Bin events into counts map | ||
| 189 | - bin = ctbin(sim.obs()) | ||
| 190 | - bin["emin"].real(emin) | ||
| 191 | - bin["emax"].real(emax) | ||
| 192 | - bin["enumbins"].integer(enumbins) | ||
| 193 | - bin["nxpix"].integer(nxpix) | ||
| 194 | - bin["nypix"].integer(nypix) | ||
| 195 | - bin["binsz"].real(binsz) | ||
| 196 | - bin["coordsys"].string(coordsys) | ||
| 197 | - bin["xref"].real(ra) | ||
| 198 | - bin["yref"].real(dec) | ||
| 199 | - bin["proj"].string(proj) | ||
| 200 | - bin.run() | ||
| 201 | - print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)" | ||
| 202 | - | ||
| 203 | - # Update timing | ||
| 204 | - wall_seconds += bin.telapse() | ||
| 205 | - cpu_seconds += bin.celapse() | ||
| 206 | - | ||
| 207 | - # Perform maximum likelihood fitting | ||
| 208 | - like = ctlike(bin.obs()) | ||
| 209 | - like.run() | ||
| 210 | - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 211 | - | ||
| 212 | - # Update timing | ||
| 213 | - wall_seconds += like.telapse() | ||
| 214 | - cpu_seconds += like.celapse() | ||
| 215 | - | ||
| 216 | - # Show total times | ||
| 217 | - print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 218 | - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 219 | - | ||
| 220 | - # Show model fitting results | ||
| 221 | - #print like.obs().models() | ||
| 222 | - | ||
| 223 | - # Plot counts | ||
| 224 | - plot_counts(bin.obs()) | ||
| 225 | - | ||
| 226 | - # Return | ||
| 227 | - return | ||
| 228 | - | ||
| 229 | - | ||
| 230 | -# =========== # | ||
| 231 | -# Plot counts # | ||
| 232 | -# =========== # | ||
| 233 | -def plot_counts(observations): | ||
| 234 | - """ | ||
| 235 | - Plot counts. | ||
| 236 | - """ | ||
| 237 | - # Only proceed if matplotlib is available | ||
| 238 | - try: | ||
| 239 | - # Import matplotlib | ||
| 240 | - import matplotlib.pyplot as plt | ||
| 241 | - | ||
| 242 | - # Set legend fontsize | ||
| 243 | - params = {'legend.fontsize': 10} | ||
| 244 | - plt.rcParams.update(params) | ||
| 245 | - | ||
| 246 | - # Set plot styles | ||
| 247 | - styles = ['b-', 'g-', 'y-', 'n-'] | ||
| 248 | - | ||
| 249 | - # Dump header | ||
| 250 | - print "" | ||
| 251 | - print "Make plots (using matplotlib):" | ||
| 252 | - print "==============================" | ||
| 253 | - | ||
| 254 | - # Create figure 1 | ||
| 255 | - plt.figure(1,figsize=(12,6)) | ||
| 256 | - plt.subplots_adjust(hspace=.7) | ||
| 257 | - | ||
| 258 | - # Create subplot 1 | ||
| 259 | - plt.subplot(121) | ||
| 260 | - plt.title("Spectrum (summed over all pixels)") | ||
| 261 | - | ||
| 262 | - # Loop over observations | ||
| 263 | - for obs in observations: | ||
| 264 | - | ||
| 265 | - # Get event cube | ||
| 266 | - cube = cast_GCTAEventCube(obs.events()) | ||
| 267 | - | ||
| 268 | - # Create energy axis | ||
| 269 | - energy = [] | ||
| 270 | - ebounds = cube.ebounds() | ||
| 271 | - for i in range(ebounds.size()): | ||
| 272 | - energy.append(ebounds.elogmean(i).TeV()) | ||
| 273 | - | ||
| 274 | - # Create spectrum | ||
| 275 | - print "Extract data:" | ||
| 276 | - counts = [0.0 for i in range(ebounds.size())] | ||
| 277 | - for bin in cube: | ||
| 278 | - index = ebounds.index(bin.energy()) | ||
| 279 | - counts[index] = counts[index] + bin.counts() | ||
| 280 | - | ||
| 281 | - # Create error bars | ||
| 282 | - error = [sqrt(c) for c in counts] | ||
| 283 | - | ||
| 284 | - # Plot spectrum | ||
| 285 | - plt.loglog(energy, counts, 'ro', label='data') | ||
| 286 | - plt.errorbar(energy, counts, error, fmt=None, ecolor='r') | ||
| 287 | - | ||
| 288 | - # Extract models | ||
| 289 | - print "Extract models:" | ||
| 290 | - sum_model = [0.0 for i in range(ebounds.size())] | ||
| 291 | - for k, m in enumerate(observations.models()): | ||
| 292 | - print "- "+m.name() | ||
| 293 | - model = [0.0 for i in range(ebounds.size())] | ||
| 294 | - for bin in cube: | ||
| 295 | - index = ebounds.index(bin.energy()) | ||
| 296 | - prob = m.eval(bin, obs) | ||
| 297 | - model[index] = model[index] + prob * bin.size() | ||
| 298 | - for i in range(ebounds.size()): | ||
| 299 | - sum_model[i] = sum_model[i] + model[i] | ||
| 300 | - plt.loglog(energy, model, styles[k], label=m.name()) | ||
| 301 | - plt.loglog(energy, sum_model, 'r-', label='total') | ||
| 302 | - | ||
| 303 | - # Put labels | ||
| 304 | - plt.xlabel("Energy (TeV)") | ||
| 305 | - plt.ylabel("Counts") | ||
| 306 | - plt.legend(loc="lower left") | ||
| 307 | - | ||
| 308 | - # Create subplot 2 | ||
| 309 | - plt.subplot(122) | ||
| 310 | - plt.title("Offset (summed over all energies)") | ||
| 311 | - | ||
| 312 | - # Set Crab direction | ||
| 313 | - crab = GSkyDir() | ||
| 314 | - crab.radec_deg(83.63, 22.01) | ||
| 315 | - | ||
| 316 | - # Loop over observations | ||
| 317 | - for obs in observations: | ||
| 318 | - | ||
| 319 | - # Get event cube | ||
| 320 | - cube = cast_GCTAEventCube(obs.events()) | ||
| 321 | - | ||
| 322 | - # Create offset histogram | ||
| 323 | - print "Extract data:" | ||
| 324 | - nx = 30 | ||
| 325 | - doffset2 = 0.01 | ||
| 326 | - offset2 = [(i+0.5)*doffset2 for i in range(nx)] | ||
| 327 | - counts = [0.0 for i in range(nx)] | ||
| 328 | - for bin in cube: | ||
| 329 | - off = bin.dir().dist_deg(crab) | ||
| 330 | - off2 = off*off | ||
| 331 | - index = int(off2/doffset2) | ||
| 332 | - if index < nx: | ||
| 333 | - counts[index] = counts[index] + bin.counts() | ||
| 334 | - | ||
| 335 | - # Create error bars | ||
| 336 | - error = [sqrt(c) for c in counts] | ||
| 337 | - | ||
| 338 | - # Plot distribution | ||
| 339 | - plt.semilogy(offset2, counts, 'ro', label='data') | ||
| 340 | - | ||
| 341 | - # Extract models | ||
| 342 | - print "Extract models:" | ||
| 343 | - sum_model = [0.0 for i in range(nx)] | ||
| 344 | - for k, m in enumerate(observations.models()): | ||
| 345 | - print "- "+m.name() | ||
| 346 | - model = [0.0 for i in range(nx)] | ||
| 347 | - for bin in cube: | ||
| 348 | - off = bin.dir().dist_deg(crab) | ||
| 349 | - off2 = off*off | ||
| 350 | - index = int(off2/doffset2) | ||
| 351 | - if index < nx: | ||
| 352 | - prob = m.eval(bin, obs) | ||
| 353 | - model[index] = model[index] + prob * bin.size() | ||
| 354 | - for i in range(nx): | ||
| 355 | - sum_model[i] = sum_model[i] + model[i] | ||
| 356 | - plt.plot(offset2, model, styles[k], label=m.name()) | ||
| 357 | - plt.plot(offset2, sum_model, 'r-', label='total') | ||
| 358 | - plt.ylim(ymin=0.1) | ||
| 359 | - | ||
| 360 | - # Put labels | ||
| 361 | - plt.xlabel("Offset (deg^2)") | ||
| 362 | - plt.ylabel("Counts") | ||
| 363 | - plt.legend(loc="upper right") | ||
| 364 | - | ||
| 365 | - # Show counts spectra | ||
| 366 | - plt.show() | ||
| 367 | - | ||
| 368 | - except: | ||
| 369 | - print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown." | ||
| 370 | - | ||
| 371 | - # Return | ||
| 372 | - return | ||
| 373 | - | ||
| 374 | - | ||
| 375 | -#==========================# | ||
| 376 | -# Main routine entry point # | ||
| 377 | -#==========================# | ||
| 378 | -if __name__ == '__main__': | ||
| 379 | - """ | ||
| 380 | - CTA binned analysis script. | ||
| 381 | - | ||
| 382 | - This script implements the analysis steps to perform a binned analysis. | ||
| 383 | - The steps include: | ||
| 384 | - - Event simulation | ||
| 385 | - - Event binning | ||
| 386 | - - Binned maximum likelihood fitting | ||
| 387 | - | ||
| 388 | - Two variations of the script are implemented. One that calls the classes | ||
| 389 | - like executables, storing the intermediate results on disk, and a second | ||
| 390 | - that uses the results of the precedent class directly without storing | ||
| 391 | - any data on disk. | ||
| 392 | - | ||
| 393 | - The second variation will also show a counts spectrum of the data after | ||
| 394 | - fitting (only if matplotlib is available on your system). Note that the | ||
| 395 | - counts spectrum is summed over the entire data space, i.e. no event selection | ||
| 396 | - around the Crab has been performed. | ||
| 397 | - """ | ||
| 398 | - # Dump header | ||
| 399 | - print "*************************************" | ||
| 400 | - print "* CTA binned analysis scripts *" | ||
| 401 | - print "*************************************" | ||
| 402 | - print "... this script will take 1-2 minutes" | ||
| 403 | - | ||
| 404 | - # Remove any existing result files | ||
| 405 | - list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] | ||
| 406 | - for files in list: | ||
| 407 | - for file in files: | ||
| 408 | - os.remove(file) | ||
| 409 | - | ||
| 410 | - # Save intermediate results on disk | ||
| 411 | - print "" | ||
| 412 | - print "Executable analysis pipeline:" | ||
| 413 | - print "=============================" | ||
| 414 | - pipeline_v1() | ||
| 415 | - | ||
| 416 | - # Do analysis in memory | ||
| 417 | - print "" | ||
| 418 | - print "In memory analysis pipeline:" | ||
| 419 | - print "============================" | ||
| 420 | - pipeline_v2() | ||
| 421 | - print "" |
test/example_unbinned.py deleted
| 1 | -#! /usr/bin/env python | ||
| 2 | -# ========================================================================== | ||
| 3 | -# This script illustrates how to build an unbinned analysis pipeline with | ||
| 4 | -# the ctatools. If matplotlib is installed, the spectrum will be displayed | ||
| 5 | -# on the screen. | ||
| 6 | -# | ||
| 7 | -# Copyright (C) 2011 Jurgen Knodlseder | ||
| 8 | -# | ||
| 9 | -# This program is free software: you can redistribute it and/or modify | ||
| 10 | -# it under the terms of the GNU General Public License as published by | ||
| 11 | -# the Free Software Foundation, either version 3 of the License, or | ||
| 12 | -# (at your option) any later version. | ||
| 13 | -# | ||
| 14 | -# This program is distributed in the hope that it will be useful, | ||
| 15 | -# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 16 | -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 17 | -# GNU General Public License for more details. | ||
| 18 | -# | ||
| 19 | -# You should have received a copy of the GNU General Public License | ||
| 20 | -# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 21 | -# | ||
| 22 | -# ========================================================================== | ||
| 23 | -from ctatools import * | ||
| 24 | -from gammalib import * | ||
| 25 | -from math import * | ||
| 26 | -import os | ||
| 27 | -import glob | ||
| 28 | -import sys | ||
| 29 | - | ||
| 30 | - | ||
| 31 | -# ============================= # | ||
| 32 | -# Analysis pipeline - version 1 # | ||
| 33 | -# ============================= # | ||
| 34 | -def pipeline_v1(): | ||
| 35 | - """ | ||
| 36 | - Unbinned analysis pipeline - save intermediate results. | ||
| 37 | - | ||
| 38 | - This function implements an analysis pipeline that successively calls | ||
| 39 | - ctobssim, ctselect and ctlike by saving the intermediate results as FITS | ||
| 40 | - files on disk. This replicates the way how the analysis would be done | ||
| 41 | - in a ftools-like approach. | ||
| 42 | - """ | ||
| 43 | - # Set script parameters | ||
| 44 | - model_name = "data/crab.xml" | ||
| 45 | - events_name = "events.fits" | ||
| 46 | - selected_events_name = "selected_events.fits" | ||
| 47 | - result_name = "results.xml" | ||
| 48 | - caldb = "irf" | ||
| 49 | - irf = "kb_E_50h_v3" | ||
| 50 | - ra = 83.63 | ||
| 51 | - dec = 22.01 | ||
| 52 | - rad_sim = 10.0 | ||
| 53 | - tstart = 0.0 | ||
| 54 | - tstop = 1800.0 | ||
| 55 | - emin = 0.1 | ||
| 56 | - emax = 100.0 | ||
| 57 | - rad_select = 3.0 | ||
| 58 | - | ||
| 59 | - # Initialise timing | ||
| 60 | - wall_seconds = 0.0 | ||
| 61 | - cpu_seconds = 0.0 | ||
| 62 | - | ||
| 63 | - # Simulate events | ||
| 64 | - sim = ctobssim() | ||
| 65 | - sim.logFileOpen() | ||
| 66 | - sim["infile"].filename(model_name) | ||
| 67 | - sim["outfile"].filename(events_name) | ||
| 68 | - sim["caldb"].string(caldb) | ||
| 69 | - sim["irf"].string(irf) | ||
| 70 | - sim["ra"].real(ra) | ||
| 71 | - sim["dec"].real(dec) | ||
| 72 | - sim["rad"].real(rad_sim) | ||
| 73 | - sim["tmin"].real(tstart) | ||
| 74 | - sim["tmax"].real(tstop) | ||
| 75 | - sim["emin"].real(emin) | ||
| 76 | - sim["emax"].real(emax) | ||
| 77 | - sim.execute() | ||
| 78 | - print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 79 | - | ||
| 80 | - # Update timing | ||
| 81 | - wall_seconds += sim.telapse() | ||
| 82 | - cpu_seconds += sim.celapse() | ||
| 83 | - | ||
| 84 | - # Select events | ||
| 85 | - select = ctselect() | ||
| 86 | - select.logFileOpen() | ||
| 87 | - select["infile"].filename(events_name) | ||
| 88 | - select["outfile"].filename(selected_events_name) | ||
| 89 | - select["ra"].real(ra) | ||
| 90 | - select["dec"].real(dec) | ||
| 91 | - select["rad"].real(rad_select) | ||
| 92 | - select["tmin"].real(tstart) | ||
| 93 | - select["tmax"].real(tstop) | ||
| 94 | - select["emin"].real(emin) | ||
| 95 | - select["emax"].real(emax) | ||
| 96 | - select.execute() | ||
| 97 | - print "Select event ("+str(select.celapse())+" CPU seconds)" | ||
| 98 | - | ||
| 99 | - # Update timing | ||
| 100 | - wall_seconds += select.telapse() | ||
| 101 | - cpu_seconds += select.celapse() | ||
| 102 | - | ||
| 103 | - # Perform maximum likelihood fitting | ||
| 104 | - like = ctlike() | ||
| 105 | - like.logFileOpen() | ||
| 106 | - like["evfile"].filename(selected_events_name) | ||
| 107 | - like["srcmdl"].filename(model_name) | ||
| 108 | - like["outmdl"].filename(result_name) | ||
| 109 | - like["method"].string("UNBINNED") | ||
| 110 | - like["caldb"].string(caldb) | ||
| 111 | - like["irf"].string(irf) | ||
| 112 | - like.execute() | ||
| 113 | - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 114 | - | ||
| 115 | - # Update timing | ||
| 116 | - wall_seconds += like.telapse() | ||
| 117 | - cpu_seconds += like.celapse() | ||
| 118 | - | ||
| 119 | - # Show total times | ||
| 120 | - print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 121 | - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 122 | - | ||
| 123 | - # Return | ||
| 124 | - return | ||
| 125 | - | ||
| 126 | - | ||
| 127 | -# ============================= # | ||
| 128 | -# Analysis pipeline - version 2 # | ||
| 129 | -# ============================= # | ||
| 130 | -def pipeline_v2(show_data): | ||
| 131 | - """ | ||
| 132 | - Unbinned analysis pipeline - keep intermediate results in memory. | ||
| 133 | - | ||
| 134 | - This function implements an analysis pipeline that successively calls | ||
| 135 | - ctobssim, ctselect and ctlike without saving the intermediate results as | ||
| 136 | - FITS files on disk. All data is only hold in memory. | ||
| 137 | - """ | ||
| 138 | - # Set script parameters | ||
| 139 | - model_name = "data/crab.xml" | ||
| 140 | - caldb = "irf" | ||
| 141 | - irf = "kb_E_50h_v3" | ||
| 142 | - ra = 83.63 | ||
| 143 | - dec = 22.01 | ||
| 144 | - rad_sim = 10.0 | ||
| 145 | - tstart = 0.0 | ||
| 146 | - tstop = 1800.0 | ||
| 147 | - emin = 0.1 | ||
| 148 | - emax = 100.0 | ||
| 149 | - rad_select = 3.0 | ||
| 150 | - | ||
| 151 | - # Initialise timing | ||
| 152 | - wall_seconds = 0.0 | ||
| 153 | - cpu_seconds = 0.0 | ||
| 154 | - | ||
| 155 | - # Simulate events | ||
| 156 | - sim = ctobssim() | ||
| 157 | - sim["infile"].filename(model_name) | ||
| 158 | - sim["caldb"].string(caldb) | ||
| 159 | - sim["irf"].string(irf) | ||
| 160 | - sim["ra"].real(ra) | ||
| 161 | - sim["dec"].real(dec) | ||
| 162 | - sim["rad"].real(rad_sim) | ||
| 163 | - sim["tmin"].real(tstart) | ||
| 164 | - sim["tmax"].real(tstop) | ||
| 165 | - sim["emin"].real(emin) | ||
| 166 | - sim["emax"].real(emax) | ||
| 167 | - sim.run() | ||
| 168 | - print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 169 | - | ||
| 170 | - # Update timing | ||
| 171 | - wall_seconds += sim.telapse() | ||
| 172 | - cpu_seconds += sim.celapse() | ||
| 173 | - | ||
| 174 | - # Select events | ||
| 175 | - select = ctselect(sim.obs()) | ||
| 176 | - select["ra"].real(ra) | ||
| 177 | - select["dec"].real(dec) | ||
| 178 | - select["rad"].real(rad_select) | ||
| 179 | - select["tmin"].real(tstart) | ||
| 180 | - select["tmax"].real(tstop) | ||
| 181 | - select["emin"].real(emin) | ||
| 182 | - select["emax"].real(emax) | ||
| 183 | - select.run() | ||
| 184 | - print "Select event ("+str(select.celapse())+" CPU seconds)" | ||
| 185 | - | ||
| 186 | - # Update timing | ||
| 187 | - wall_seconds += select.telapse() | ||
| 188 | - cpu_seconds += select.celapse() | ||
| 189 | - | ||
| 190 | - # Perform maximum likelihood fitting | ||
| 191 | - like = ctlike(select.obs()) | ||
| 192 | - like.run() | ||
| 193 | - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 194 | - | ||
| 195 | - # Update timing | ||
| 196 | - wall_seconds += like.telapse() | ||
| 197 | - cpu_seconds += like.celapse() | ||
| 198 | - | ||
| 199 | - # Show total times | ||
| 200 | - print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 201 | - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 202 | - | ||
| 203 | - # Show model fitting results | ||
| 204 | - #print like.obs().models() | ||
| 205 | - | ||
| 206 | - # Optionally plot counts | ||
| 207 | - if show_data: | ||
| 208 | - plot_counts(like.obs()) | ||
| 209 | - | ||
| 210 | - # Return | ||
| 211 | - return | ||
| 212 | - | ||
| 213 | - | ||
| 214 | -# =========== # | ||
| 215 | -# Plot counts # | ||
| 216 | -# =========== # | ||
| 217 | -def plot_counts(observations): | ||
| 218 | - """ | ||
| 219 | - Plot counts. | ||
| 220 | - """ | ||
| 221 | - # Only proceed if matplotlib is available | ||
| 222 | - try: | ||
| 223 | - # Import matplotlib | ||
| 224 | - import matplotlib.pyplot as plt | ||
| 225 | - | ||
| 226 | - # Set legend fontsize | ||
| 227 | - params = {'legend.fontsize': 10} | ||
| 228 | - plt.rcParams.update(params) | ||
| 229 | - | ||
| 230 | - # Set plot styles | ||
| 231 | - styles = ['b-', 'g-', 'y-', 'n-'] | ||
| 232 | - | ||
| 233 | - # Dump header | ||
| 234 | - print "" | ||
| 235 | - print "Make plots (using matplotlib):" | ||
| 236 | - print "==============================" | ||
| 237 | - | ||
| 238 | - # Create figure 1 | ||
| 239 | - plt.figure(1,figsize=(12,6)) | ||
| 240 | - plt.subplots_adjust(hspace=.7) | ||
| 241 | - | ||
| 242 | - # Create subplot 1 | ||
| 243 | - plt.subplot(121) | ||
| 244 | - plt.title("Spectrum (summed over all pixels)") | ||
| 245 | - | ||
| 246 | - # Loop over observations | ||
| 247 | - for obs in observations: | ||
| 248 | - | ||
| 249 | - # Get event list | ||
| 250 | - list = cast_GCTAEventList(obs.events()) | ||
| 251 | - | ||
| 252 | - # Create energy axis | ||
| 253 | - ebounds = GEbounds() | ||
| 254 | - emin = GEnergy() | ||
| 255 | - emax = GEnergy() | ||
| 256 | - emin.TeV(0.1) | ||
| 257 | - emax.TeV(100.0) | ||
| 258 | - ebounds.setlog(emin, emax, 10) | ||
| 259 | - energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())] | ||
| 260 | - | ||
| 261 | - # Create spectrum | ||
| 262 | - print "Extract data:" | ||
| 263 | - counts = [0.0 for i in range(ebounds.size())] | ||
| 264 | - print list.size() | ||
| 265 | - for atom in list: | ||
| 266 | - index = ebounds.index(atom.energy()) | ||
| 267 | - counts[index] = counts[index] + 1.0 | ||
| 268 | - | ||
| 269 | - # Create error bars | ||
| 270 | - error = [sqrt(c) for c in counts] | ||
| 271 | - | ||
| 272 | - # Plot spectrum | ||
| 273 | - plt.loglog(energy, counts, 'ro', label='data') | ||
| 274 | - #plt.errorbar(energy, counts, error, fmt=None, ecolor='r') | ||
| 275 | - | ||
| 276 | - # Extract models | ||
| 277 | - print "Extract models:" | ||
| 278 | - sum_model = [0.0 for i in range(ebounds.size())] | ||
| 279 | - for k, m in enumerate(observations.models()): | ||
| 280 | - print "- "+m.name() | ||
| 281 | - model = [0.0 for i in range(ebounds.size())] | ||
| 282 | - for atom in list: | ||
| 283 | - index = ebounds.index(atom.energy()) | ||
| 284 | - prob = m.eval(atom, obs) | ||
| 285 | - model[index] = model[index] + prob * atom.size() | ||
| 286 | - for i in range(ebounds.size()): | ||
| 287 | - sum_model[i] = sum_model[i] + model[i] | ||
| 288 | - #plt.loglog(energy, model, styles[k], label=m.name()) | ||
| 289 | - #plt.loglog(energy, sum_model, 'r-', label='total') | ||
| 290 | - | ||
| 291 | - # Put labels | ||
| 292 | - plt.xlabel("Energy (TeV)") | ||
| 293 | - plt.ylabel("Counts") | ||
| 294 | - plt.legend(loc="lower left") | ||
| 295 | - | ||
| 296 | - # Create subplot 2 | ||
| 297 | - plt.subplot(122) | ||
| 298 | - plt.title("Offset (summed over all energies)") | ||
| 299 | - | ||
| 300 | - # Set Crab direction | ||
| 301 | - crab = GSkyDir() | ||
| 302 | - crab.radec_deg(83.63, 22.01) | ||
| 303 | - | ||
| 304 | - # Loop over observations | ||
| 305 | - for obs in observations: | ||
| 306 | - | ||
| 307 | - # Get event cube | ||
| 308 | - list = cast_GCTAEventList(obs.events()) | ||
| 309 | - | ||
| 310 | - # Create offset histogram | ||
| 311 | - print "Extract data:" | ||
| 312 | - nx = 30 | ||
| 313 | - doffset2 = 0.01 | ||
| 314 | - offset2 = [(i+0.5)*doffset2 for i in range(nx)] | ||
| 315 | - counts = [0.0 for i in range(nx)] | ||
| 316 | - for atom in list: | ||
| 317 | - off = atom.dir().dist_deg(crab) | ||
| 318 | - off2 = off*off | ||
| 319 | - index = int(off2/doffset2) | ||
| 320 | - if index < nx: | ||
| 321 | - counts[index] = counts[index] + 1.0 | ||
| 322 | - | ||
| 323 | - # Create error bars | ||
| 324 | - error = [sqrt(c) for c in counts] | ||
| 325 | - | ||
| 326 | - # Plot distribution | ||
| 327 | - plt.semilogy(offset2, counts, 'ro', label='data') | ||
| 328 | - | ||
| 329 | - # Extract models | ||
| 330 | - print "Extract models:" | ||
| 331 | - sum_model = [0.0 for i in range(nx)] | ||
| 332 | - for k, m in enumerate(observations.models()): | ||
| 333 | - print "- "+m.name() | ||
| 334 | - model = [0.0 for i in range(nx)] | ||
| 335 | - for atom in list: | ||
| 336 | - off = atom.dir().dist_deg(crab) | ||
| 337 | - off2 = off*off | ||
| 338 | - index = int(off2/doffset2) | ||
| 339 | - if index < nx: | ||
| 340 | - prob = m.eval(atom, obs) | ||
| 341 | - model[index] = model[index] + prob * atom.size() | ||
| 342 | - for i in range(nx): | ||
| 343 | - sum_model[i] = sum_model[i] + model[i] | ||
| 344 | - #plt.plot(offset2, model, styles[k], label=m.name()) | ||
| 345 | - #plt.plot(offset2, sum_model, 'r-', label='total') | ||
| 346 | - #plt.ylim(ymin=0.1) | ||
| 347 | - | ||
| 348 | - # Put labels | ||
| 349 | - plt.xlabel("Offset (deg^2)") | ||
| 350 | - plt.ylabel("Counts") | ||
| 351 | - plt.legend(loc="upper right") | ||
| 352 | - | ||
| 353 | - # Show counts spectra | ||
| 354 | - plt.show() | ||
| 355 | - | ||
| 356 | - except ImportError: | ||
| 357 | - print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown." | ||
| 358 | - | ||
| 359 | - # Return | ||
| 360 | - return | ||
| 361 | - | ||
| 362 | - | ||
| 363 | -#==========================# | ||
| 364 | -# Main routine entry point # | ||
| 365 | -#==========================# | ||
| 366 | -if __name__ == '__main__': | ||
| 367 | - """ | ||
| 368 | - CTA unbinned analysis script. | ||
| 369 | - | ||
| 370 | - This script implements the analysis steps to perform an unbinned analysis. | ||
| 371 | - The steps include: | ||
| 372 | - - Event simulation | ||
| 373 | - - Event selection | ||
| 374 | - - Binned maximum likelihood fitting | ||
| 375 | - | ||
| 376 | - Two variations of the script are implemented. One that calls the classes | ||
| 377 | - like executables, storing the intermediate results on disk, and a second | ||
| 378 | - that uses the results of the precedent class directly without storing | ||
| 379 | - any data on disk. | ||
| 380 | - """ | ||
| 381 | - # Initialise flags | ||
| 382 | - need_help = False | ||
| 383 | - show_data = True | ||
| 384 | - | ||
| 385 | - # Test for command line arguments | ||
| 386 | - print sys.argv[0] | ||
| 387 | - if (len(sys.argv) > 1): | ||
| 388 | - if sys.argv[1] == "-h": | ||
| 389 | - need_help = True | ||
| 390 | - elif sys.argv[1] == "-noshow": | ||
| 391 | - show_data = False | ||
| 392 | - else: | ||
| 393 | - need_help = True | ||
| 394 | - | ||
| 395 | - # Print help if needed and exit | ||
| 396 | - if need_help: | ||
| 397 | - print "Usage: example_binned.py [OPTIONS]" | ||
| 398 | - print " -h Display this usage message" | ||
| 399 | - print " -noshow Do not show data using matplotlib" | ||
| 400 | - sys.exit() | ||
| 401 | - | ||
| 402 | - # Dump header | ||
| 403 | - print "*********************************" | ||
| 404 | - print "* CTA unbinned analysis scripts *" | ||
| 405 | - print "*********************************" | ||
| 406 | - | ||
| 407 | - # Remove any existing result files | ||
| 408 | - list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] | ||
| 409 | - for files in list: | ||
| 410 | - for file in files: | ||
| 411 | - os.remove(file) | ||
| 412 | - | ||
| 413 | - # Save intermediate results on disk | ||
| 414 | - print "" | ||
| 415 | - print "Executable analysis pipeline:" | ||
| 416 | - print "=============================" | ||
| 417 | - pipeline_v1() | ||
| 418 | - | ||
| 419 | - # Do analysis in memory | ||
| 420 | - print "" | ||
| 421 | - print "In memory analysis pipeline:" | ||
| 422 | - print "============================" | ||
| 423 | - pipeline_v2(show_data) | ||
| 424 | - print "" |