Commit 5986e22e9c64a1ecd6c8d9634cc2d94acd0b6d6a
1 parent
ee3661aa
initial version
Showing
2 changed files
with
445 additions
and
0 deletions
test/Makefile.am
0 → 100644
| 1 | +## Process this file with automake to produce Makefile.in | ||
| 2 | + | ||
| 3 | +# we use srcdir/src for "private" includes | ||
| 4 | +# we use srcdir/include for the static global include files | ||
| 5 | +# we use builddir/src for the generated include files to get found | ||
| 6 | + | ||
| 7 | +# Optional python bindings | ||
| 8 | +if ENABLE_PYTHON_BINDING | ||
| 9 | +if HAVE_PYTHON | ||
| 10 | +if WITH_PYTHON_BINDING | ||
| 11 | + TEST_PYTHON_SCRIPT = test_python.py | ||
| 12 | + TEST_PYTHON_ENV = PYTHONPATH=$(top_builddir)/pyext/build:$(PYTHONPATH) | ||
| 13 | +endif | ||
| 14 | +endif | ||
| 15 | +endif | ||
| 16 | + | ||
| 17 | +# Include directories | ||
| 18 | +INCLUDES = -I$(top_srcdir)/include \ | ||
| 19 | + -I$(top_builddir)/src | ||
| 20 | + | ||
| 21 | +# Test programs (those will be compiled) | ||
| 22 | +check_PROGRAMS = | ||
| 23 | + | ||
| 24 | +# Set test environment | ||
| 25 | +TESTS_ENVIRONMENT = LD_LIBRARY_PATH=$(top_builddir)/src/.libs:$(LD_LIBRARY_PATH) \ | ||
| 26 | + $(TEST_PYTHON_ENV) | ||
| 27 | + | ||
| 28 | +# Tests (these are the tests to be executed) | ||
| 29 | +TESTS = test_ctatools.sh \ | ||
| 30 | + $(TEST_PYTHON_SCRIPT) | ||
| 31 | + | ||
| 32 | +# Standard test sources and links | ||
| 33 | +#test_GVector_SOURCES = test_GVector.cpp | ||
| 34 | +#test_GVector_LDFLAGS = @LDFLAGS@ | ||
| 35 | +#test_GVector_CPPFLAGS = @CPPFLAGS@ | ||
| 36 | +#test_GVector_LDADD = $(top_srcdir)/src/libgamma.la |
test/example_unbinned.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ===========================================================================================# | ||
| 3 | +# This script illustrates how to build a unbinned analysis chain using CTAtools. | ||
| 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 | +# Analysis pipeline - version 1 # | ||
| 17 | +# ============================= # | ||
| 18 | +def pipeline_v1(): | ||
| 19 | + """ | ||
| 20 | + Unbinned analysis pipeline - save intermediate results. | ||
| 21 | + | ||
| 22 | + This function implements an analysis pipeline that successively calls | ||
| 23 | + ctobssim, ctselect and ctlike by saving the intermediate results as FITS | ||
| 24 | + files on disk. This replicates the way how the analysis would be done | ||
| 25 | + in a ftools-like approach. | ||
| 26 | + """ | ||
| 27 | + # Set script parameters | ||
| 28 | + model_name = "data/crab.xml" | ||
| 29 | + events_name = "events.fits" | ||
| 30 | + selected_events_name = "selected_events.fits" | ||
| 31 | + result_name = "results.xml" | ||
| 32 | + caldb = "irf" | ||
| 33 | + irf = "kb_E_50h_v3" | ||
| 34 | + ra = 83.63 | ||
| 35 | + dec = 22.01 | ||
| 36 | + rad_sim = 10.0 | ||
| 37 | + tstart = 0.0 | ||
| 38 | + tstop = 1800.0 | ||
| 39 | + emin = 0.1 | ||
| 40 | + emax = 100.0 | ||
| 41 | + rad_select = 3.0 | ||
| 42 | + | ||
| 43 | + # Initialise timing | ||
| 44 | + wall_seconds = 0.0 | ||
| 45 | + cpu_seconds = 0.0 | ||
| 46 | + | ||
| 47 | + # Simulate events | ||
| 48 | + sim = ctobssim() | ||
| 49 | + sim.logFileOpen() | ||
| 50 | + sim["infile"].filename(model_name) | ||
| 51 | + sim["outfile"].filename(events_name) | ||
| 52 | + sim["caldb"].string(caldb) | ||
| 53 | + sim["irf"].string(irf) | ||
| 54 | + sim["ra"].real(ra) | ||
| 55 | + sim["dec"].real(dec) | ||
| 56 | + sim["rad"].real(rad_sim) | ||
| 57 | + sim["tmin"].real(tstart) | ||
| 58 | + sim["tmax"].real(tstop) | ||
| 59 | + sim["emin"].real(emin) | ||
| 60 | + sim["emax"].real(emax) | ||
| 61 | + sim.execute() | ||
| 62 | + print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 63 | + | ||
| 64 | + # Update timing | ||
| 65 | + wall_seconds += sim.telapse() | ||
| 66 | + cpu_seconds += sim.celapse() | ||
| 67 | + | ||
| 68 | + # Select events | ||
| 69 | + select = ctselect() | ||
| 70 | + select.logFileOpen() | ||
| 71 | + select["infile"].filename(events_name) | ||
| 72 | + select["outfile"].filename(selected_events_name) | ||
| 73 | + select["ra"].real(ra) | ||
| 74 | + select["dec"].real(dec) | ||
| 75 | + select["rad"].real(rad_select) | ||
| 76 | + select["tmin"].real(tstart) | ||
| 77 | + select["tmax"].real(tstop) | ||
| 78 | + select["emin"].real(emin) | ||
| 79 | + select["emax"].real(emax) | ||
| 80 | + select.execute() | ||
| 81 | + print "Select event ("+str(select.celapse())+" CPU seconds)" | ||
| 82 | + | ||
| 83 | + # Update timing | ||
| 84 | + wall_seconds += select.telapse() | ||
| 85 | + cpu_seconds += select.celapse() | ||
| 86 | + | ||
| 87 | + # Perform maximum likelihood fitting | ||
| 88 | + like = ctlike() | ||
| 89 | + like.logFileOpen() | ||
| 90 | + like["evfile"].filename(selected_events_name) | ||
| 91 | + like["srcmdl"].filename(model_name) | ||
| 92 | + like["outmdl"].filename(result_name) | ||
| 93 | + like["method"].string("UNBINNED") | ||
| 94 | + like["caldb"].string(caldb) | ||
| 95 | + like["irf"].string(irf) | ||
| 96 | + like.execute() | ||
| 97 | + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 98 | + | ||
| 99 | + # Update timing | ||
| 100 | + wall_seconds += like.telapse() | ||
| 101 | + cpu_seconds += like.celapse() | ||
| 102 | + | ||
| 103 | + # Show total times | ||
| 104 | + print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 105 | + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 106 | + | ||
| 107 | + # Return | ||
| 108 | + return | ||
| 109 | + | ||
| 110 | + | ||
| 111 | +# ============================= # | ||
| 112 | +# Analysis pipeline - version 2 # | ||
| 113 | +# ============================= # | ||
| 114 | +def pipeline_v2(show_data): | ||
| 115 | + """ | ||
| 116 | + Unbinned analysis pipeline - keep intermediate results in memory. | ||
| 117 | + | ||
| 118 | + This function implements an analysis pipeline that successively calls | ||
| 119 | + ctobssim, ctselect and ctlike without saving the intermediate results as | ||
| 120 | + FITS files on disk. All data is only hold in memory. | ||
| 121 | + """ | ||
| 122 | + # Set script parameters | ||
| 123 | + model_name = "data/crab.xml" | ||
| 124 | + caldb = "irf" | ||
| 125 | + irf = "kb_E_50h_v3" | ||
| 126 | + ra = 83.63 | ||
| 127 | + dec = 22.01 | ||
| 128 | + rad_sim = 10.0 | ||
| 129 | + tstart = 0.0 | ||
| 130 | + tstop = 1800.0 | ||
| 131 | + emin = 0.1 | ||
| 132 | + emax = 100.0 | ||
| 133 | + rad_select = 3.0 | ||
| 134 | + | ||
| 135 | + # Initialise timing | ||
| 136 | + wall_seconds = 0.0 | ||
| 137 | + cpu_seconds = 0.0 | ||
| 138 | + | ||
| 139 | + # Simulate events | ||
| 140 | + sim = ctobssim() | ||
| 141 | + sim["infile"].filename(model_name) | ||
| 142 | + sim["caldb"].string(caldb) | ||
| 143 | + sim["irf"].string(irf) | ||
| 144 | + sim["ra"].real(ra) | ||
| 145 | + sim["dec"].real(dec) | ||
| 146 | + sim["rad"].real(rad_sim) | ||
| 147 | + sim["tmin"].real(tstart) | ||
| 148 | + sim["tmax"].real(tstop) | ||
| 149 | + sim["emin"].real(emin) | ||
| 150 | + sim["emax"].real(emax) | ||
| 151 | + sim.run() | ||
| 152 | + print "Simulated events ("+str(sim.celapse())+" CPU seconds)" | ||
| 153 | + | ||
| 154 | + # Update timing | ||
| 155 | + wall_seconds += sim.telapse() | ||
| 156 | + cpu_seconds += sim.celapse() | ||
| 157 | + | ||
| 158 | + # Select events | ||
| 159 | + select = ctselect(sim.obs()) | ||
| 160 | + select["ra"].real(ra) | ||
| 161 | + select["dec"].real(dec) | ||
| 162 | + select["rad"].real(rad_select) | ||
| 163 | + select["tmin"].real(tstart) | ||
| 164 | + select["tmax"].real(tstop) | ||
| 165 | + select["emin"].real(emin) | ||
| 166 | + select["emax"].real(emax) | ||
| 167 | + select.run() | ||
| 168 | + print "Select event ("+str(select.celapse())+" CPU seconds)" | ||
| 169 | + | ||
| 170 | + # Update timing | ||
| 171 | + wall_seconds += select.telapse() | ||
| 172 | + cpu_seconds += select.celapse() | ||
| 173 | + | ||
| 174 | + # Perform maximum likelihood fitting | ||
| 175 | + like = ctlike(select.obs()) | ||
| 176 | + like.run() | ||
| 177 | + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)" | ||
| 178 | + | ||
| 179 | + # Update timing | ||
| 180 | + wall_seconds += like.telapse() | ||
| 181 | + cpu_seconds += like.celapse() | ||
| 182 | + | ||
| 183 | + # Show total times | ||
| 184 | + print "Total wall time elapsed: "+str(wall_seconds)+" seconds" | ||
| 185 | + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds" | ||
| 186 | + | ||
| 187 | + # Show model fitting results | ||
| 188 | + #print like.obs().models() | ||
| 189 | + | ||
| 190 | + # Optionally plot counts | ||
| 191 | + if show_data: | ||
| 192 | + plot_counts(like.obs()) | ||
| 193 | + | ||
| 194 | + # Return | ||
| 195 | + return | ||
| 196 | + | ||
| 197 | + | ||
| 198 | +# =========== # | ||
| 199 | +# Plot counts # | ||
| 200 | +# =========== # | ||
| 201 | +def plot_counts(observations): | ||
| 202 | + """ | ||
| 203 | + Plot counts. | ||
| 204 | + """ | ||
| 205 | + # Only proceed if matplotlib is available | ||
| 206 | + if 0 == 0: | ||
| 207 | + #try: | ||
| 208 | + # Import matplotlib | ||
| 209 | + import matplotlib.pyplot as plt | ||
| 210 | + | ||
| 211 | + # Set legend fontsize | ||
| 212 | + params = {'legend.fontsize': 10} | ||
| 213 | + plt.rcParams.update(params) | ||
| 214 | + | ||
| 215 | + # Set plot styles | ||
| 216 | + styles = ['b-', 'g-', 'y-', 'n-'] | ||
| 217 | + | ||
| 218 | + # Dump header | ||
| 219 | + print "" | ||
| 220 | + print "Make plots (using matplotlib):" | ||
| 221 | + print "==============================" | ||
| 222 | + | ||
| 223 | + # Create figure 1 | ||
| 224 | + plt.figure(1,figsize=(12,6)) | ||
| 225 | + plt.subplots_adjust(hspace=.7) | ||
| 226 | + | ||
| 227 | + # Create subplot 1 | ||
| 228 | + plt.subplot(121) | ||
| 229 | + plt.title("Spectrum (summed over all pixels)") | ||
| 230 | + | ||
| 231 | + # Loop over observations | ||
| 232 | + for obs in observations: | ||
| 233 | + | ||
| 234 | + # Get event list | ||
| 235 | + list = cast_GCTAEventList(obs.events()) | ||
| 236 | + | ||
| 237 | + # Create energy axis | ||
| 238 | + ebounds = GEbounds() | ||
| 239 | + emin = GEnergy() | ||
| 240 | + emax = GEnergy() | ||
| 241 | + emin.TeV(0.1) | ||
| 242 | + emax.TeV(100.0) | ||
| 243 | + ebounds.setlog(emin, emax, 10) | ||
| 244 | + energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())] | ||
| 245 | + | ||
| 246 | + # Create spectrum | ||
| 247 | + print "Extract data:" | ||
| 248 | + counts = [0.0 for i in range(ebounds.size())] | ||
| 249 | + print list.size() | ||
| 250 | + for atom in list: | ||
| 251 | + index = ebounds.index(atom.energy()) | ||
| 252 | + counts[index] = counts[index] + 1.0 | ||
| 253 | + | ||
| 254 | + # Create error bars | ||
| 255 | + error = [sqrt(c) for c in counts] | ||
| 256 | + | ||
| 257 | + # Plot spectrum | ||
| 258 | + plt.loglog(energy, counts, 'ro', label='data') | ||
| 259 | + #plt.errorbar(energy, counts, error, fmt=None, ecolor='r') | ||
| 260 | + | ||
| 261 | + # Extract models | ||
| 262 | + print "Extract models:" | ||
| 263 | + sum_model = [0.0 for i in range(ebounds.size())] | ||
| 264 | + for k, m in enumerate(observations.models()): | ||
| 265 | + print "- "+m.name() | ||
| 266 | + model = [0.0 for i in range(ebounds.size())] | ||
| 267 | + for atom in list: | ||
| 268 | + index = ebounds.index(atom.energy()) | ||
| 269 | + prob = m.eval(atom, obs) | ||
| 270 | + model[index] = model[index] + prob * atom.size() | ||
| 271 | + for i in range(ebounds.size()): | ||
| 272 | + sum_model[i] = sum_model[i] + model[i] | ||
| 273 | + #plt.loglog(energy, model, styles[k], label=m.name()) | ||
| 274 | + #plt.loglog(energy, sum_model, 'r-', label='total') | ||
| 275 | + | ||
| 276 | + # Put labels | ||
| 277 | + plt.xlabel("Energy (TeV)") | ||
| 278 | + plt.ylabel("Counts") | ||
| 279 | + plt.legend(loc="lower left") | ||
| 280 | + | ||
| 281 | + # Create subplot 2 | ||
| 282 | + plt.subplot(122) | ||
| 283 | + plt.title("Offset (summed over all energies)") | ||
| 284 | + | ||
| 285 | + # Set Crab direction | ||
| 286 | + crab = GSkyDir() | ||
| 287 | + crab.radec_deg(83.63, 22.01) | ||
| 288 | + | ||
| 289 | + # Loop over observations | ||
| 290 | + for obs in observations: | ||
| 291 | + | ||
| 292 | + # Get event cube | ||
| 293 | + list = cast_GCTAEventList(obs.events()) | ||
| 294 | + | ||
| 295 | + # Create offset histogram | ||
| 296 | + print "Extract data:" | ||
| 297 | + nx = 30 | ||
| 298 | + doffset2 = 0.01 | ||
| 299 | + offset2 = [(i+0.5)*doffset2 for i in range(nx)] | ||
| 300 | + counts = [0.0 for i in range(nx)] | ||
| 301 | + for atom in list: | ||
| 302 | + off = atom.dir().dist_deg(crab) | ||
| 303 | + off2 = off*off | ||
| 304 | + index = int(off2/doffset2) | ||
| 305 | + if index < nx: | ||
| 306 | + counts[index] = counts[index] + 1.0 | ||
| 307 | + | ||
| 308 | + # Create error bars | ||
| 309 | + error = [sqrt(c) for c in counts] | ||
| 310 | + | ||
| 311 | + # Plot distribution | ||
| 312 | + plt.semilogy(offset2, counts, 'ro', label='data') | ||
| 313 | + | ||
| 314 | + # Extract models | ||
| 315 | + print "Extract models:" | ||
| 316 | + sum_model = [0.0 for i in range(nx)] | ||
| 317 | + for k, m in enumerate(observations.models()): | ||
| 318 | + print "- "+m.name() | ||
| 319 | + model = [0.0 for i in range(nx)] | ||
| 320 | + for atom in list: | ||
| 321 | + off = atom.dir().dist_deg(crab) | ||
| 322 | + off2 = off*off | ||
| 323 | + index = int(off2/doffset2) | ||
| 324 | + if index < nx: | ||
| 325 | + prob = m.eval(atom, obs) | ||
| 326 | + model[index] = model[index] + prob * atom.size() | ||
| 327 | + for i in range(nx): | ||
| 328 | + sum_model[i] = sum_model[i] + model[i] | ||
| 329 | + #plt.plot(offset2, model, styles[k], label=m.name()) | ||
| 330 | + #plt.plot(offset2, sum_model, 'r-', label='total') | ||
| 331 | + #plt.ylim(ymin=0.1) | ||
| 332 | + | ||
| 333 | + # Put labels | ||
| 334 | + plt.xlabel("Offset (deg^2)") | ||
| 335 | + plt.ylabel("Counts") | ||
| 336 | + plt.legend(loc="upper right") | ||
| 337 | + | ||
| 338 | + # Show counts spectra | ||
| 339 | + plt.show() | ||
| 340 | + | ||
| 341 | + #except: | ||
| 342 | + # print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown." | ||
| 343 | + | ||
| 344 | + # Return | ||
| 345 | + return | ||
| 346 | + | ||
| 347 | + | ||
| 348 | +#==========================# | ||
| 349 | +# Main routine entry point # | ||
| 350 | +#==========================# | ||
| 351 | +if __name__ == '__main__': | ||
| 352 | + """ | ||
| 353 | + CTA unbinned analysis script. | ||
| 354 | + | ||
| 355 | + This script implements the analysis steps to perform an unbinned analysis. | ||
| 356 | + The steps include: | ||
| 357 | + - Event simulation | ||
| 358 | + - Event selection | ||
| 359 | + - Binned maximum likelihood fitting | ||
| 360 | + | ||
| 361 | + Two variations of the script are implemented. One that calls the classes | ||
| 362 | + like executables, storing the intermediate results on disk, and a second | ||
| 363 | + that uses the results of the precedent class directly without storing | ||
| 364 | + any data on disk. | ||
| 365 | + """ | ||
| 366 | + # Initialise flags | ||
| 367 | + need_help = False | ||
| 368 | + show_data = True | ||
| 369 | + | ||
| 370 | + # Test for command line arguments | ||
| 371 | + print sys.argv[0] | ||
| 372 | + if (len(sys.argv) > 1): | ||
| 373 | + if sys.argv[1] == "-h": | ||
| 374 | + need_help = True | ||
| 375 | + elif sys.argv[1] == "-noshow": | ||
| 376 | + show_data = False | ||
| 377 | + else: | ||
| 378 | + need_help = True | ||
| 379 | + | ||
| 380 | + # Print help if needed and exit | ||
| 381 | + if need_help: | ||
| 382 | + print "Usage: example_binned.py [OPTIONS]" | ||
| 383 | + print " -h Display this usage message" | ||
| 384 | + print " -noshow Do not show data using matplotlib" | ||
| 385 | + sys.exit() | ||
| 386 | + | ||
| 387 | + # Dump header | ||
| 388 | + print "*********************************" | ||
| 389 | + print "* CTA unbinned analysis scripts *" | ||
| 390 | + print "*********************************" | ||
| 391 | + | ||
| 392 | + # Remove any existing result files | ||
| 393 | + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] | ||
| 394 | + for files in list: | ||
| 395 | + for file in files: | ||
| 396 | + os.remove(file) | ||
| 397 | + | ||
| 398 | + # Save intermediate results on disk | ||
| 399 | + print "" | ||
| 400 | + print "Executable analysis pipeline:" | ||
| 401 | + print "=============================" | ||
| 402 | + pipeline_v1() | ||
| 403 | + | ||
| 404 | + # Do analysis in memory | ||
| 405 | + print "" | ||
| 406 | + print "In memory analysis pipeline:" | ||
| 407 | + print "============================" | ||
| 408 | + pipeline_v2(show_data) | ||
| 409 | + print "" |