Commit 0847394c57541fad7444abdfc8dc62c6c11cd418
1 parent
651d45a7
Enhance ctbin Python module test.
Showing
4 changed files
with
228 additions
and
129 deletions
cscripts/obsutils.py
| 1 | 1 | # ========================================================================== |
| 2 | -# This script provides a number of functions that are useful for handling | |
| 3 | -# CTA observations. | |
| 2 | +# CTA observation handling support functions. | |
| 4 | 3 | # |
| 5 | -# Copyright (C) 2011-2015 Juergen Knoedlseder | |
| 4 | +# Copyright (C) 2011-2016 Juergen Knoedlseder | |
| 6 | 5 | # |
| 7 | 6 | # This program is free software: you can redistribute it and/or modify |
| 8 | 7 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -31,21 +30,28 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, |
| 31 | 30 | """ |
| 32 | 31 | Simulate events for all observations in the container. |
| 33 | 32 | |
| 34 | - Parameters: | |
| 35 | - obs - Observation container | |
| 36 | - Keywords: | |
| 37 | - log - Create log file(s) | |
| 38 | - debug - Create console dump? | |
| 39 | - edisp - Apply energy dispersion? | |
| 40 | - seed - Seed value for simulations (default: 0) | |
| 41 | - nbins - Number of energy bins (default: 0=unbinned) | |
| 42 | - binsz - Pixel size for binned simulation (deg/pixel) | |
| 43 | - npix - Number of pixels in X and Y for binned simulation | |
| 33 | + Args: | |
| 34 | + obs: Observation container without events | |
| 35 | + | |
| 36 | + Kwargs: | |
| 37 | + log: Create log file(s) | |
| 38 | + debug: Create console dump? | |
| 39 | + chatter: Chatter level | |
| 40 | + edisp: Apply energy dispersion? | |
| 41 | + seed: Seed value for simulations (default: 0) | |
| 42 | + nbins: Number of energy bins (default: 0=unbinned) | |
| 43 | + binsz: Pixel size for binned simulation (deg/pixel) | |
| 44 | + npix: Number of pixels in X and Y for binned simulation | |
| 45 | + proj: Projection for binned simulation | |
| 46 | + coord: Coordinate system for binned simulation | |
| 47 | + | |
| 48 | + Returns: | |
| 49 | + Observation container, filled with simulated events. | |
| 44 | 50 | """ |
| 45 | 51 | |
| 46 | 52 | # Allocate ctobssim application and set parameters |
| 47 | 53 | sim = ctools.ctobssim(obs) |
| 48 | - sim["seed"] = seed | |
| 54 | + sim["seed"] = seed | |
| 49 | 55 | sim["edisp"] = edisp |
| 50 | 56 | |
| 51 | 57 | # Optionally open the log file |
| ... | ... | @@ -85,16 +91,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, |
| 85 | 91 | |
| 86 | 92 | # Allocate ctbin application and set parameters |
| 87 | 93 | bin = ctools.ctbin(sim.obs()) |
| 88 | - bin["ebinalg"] = "LOG" | |
| 89 | - bin["emin"] = emin | |
| 90 | - bin["emax"] = emax | |
| 94 | + bin["ebinalg"] = "LOG" | |
| 95 | + bin["emin"] = emin | |
| 96 | + bin["emax"] = emax | |
| 91 | 97 | bin["enumbins"] = nbins |
| 92 | - bin["usepnt"] = True # Use pointing for map centre | |
| 93 | - bin["nxpix"] = npix | |
| 94 | - bin["nypix"] = npix | |
| 95 | - bin["binsz"] = binsz | |
| 98 | + bin["usepnt"] = True # Use pointing for map centre | |
| 99 | + bin["nxpix"] = npix | |
| 100 | + bin["nypix"] = npix | |
| 101 | + bin["binsz"] = binsz | |
| 96 | 102 | bin["coordsys"] = coord |
| 97 | - bin["proj"] = proj | |
| 103 | + bin["proj"] = proj | |
| 98 | 104 | |
| 99 | 105 | # Optionally open the log file |
| 100 | 106 | if log: |
| ... | ... | @@ -350,28 +356,30 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ |
| 350 | 356 | # ======================= # |
| 351 | 357 | def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ |
| 352 | 358 | emin=0.1, emax=100.0, rad=5.0, \ |
| 353 | - irf="South_50h", caldb="prod2", id="000000", instrument="CTA"): | |
| 359 | + irf="South_50h", caldb="prod2", id="000000"): | |
| 354 | 360 | """ |
| 355 | - Returns a single CTA observation containing an empty CTA event list. | |
| 356 | - By looping over this function you can add CTA observations to the | |
| 357 | - observation container. | |
| 358 | - | |
| 359 | - Parameters: | |
| 360 | - pntdir - Pointing direction [GSkyDir] | |
| 361 | - Keywords: | |
| 362 | - tstart - Start time [seconds] (default: 0.0) | |
| 363 | - duration - Duration of observation [seconds] (default: 1800.0) | |
| 364 | - deadc - Deadtime correction factor (default: 0.95) | |
| 365 | - emin - Minimum event energy [TeV] (default: 0.1) | |
| 366 | - emax - Maximum event energy [TeV] (default: 100.0) | |
| 367 | - rad - ROI radius used for analysis [deg] (default: 5.0) | |
| 368 | - irf - Instrument response function (default: cta_dummy_irf) | |
| 369 | - caldb - Calibration database path (default: "dummy") | |
| 370 | - id - Run identifier (default: "000000") | |
| 371 | - instrument - Intrument (default: "CTA") | |
| 361 | + Set a single CTA observation. | |
| 362 | + | |
| 363 | + The function sets a single CTA observation containing an empty CTA | |
| 364 | + event list. By looping over this function you can add CTA observations | |
| 365 | + to the observation container. | |
| 366 | + | |
| 367 | + Args: | |
| 368 | + pntdir - Pointing direction [GSkyDir] | |
| 369 | + | |
| 370 | + Kwargs: | |
| 371 | + tstart: Start time (seconds) (default: 0.0) | |
| 372 | + duration: Duration of observation (seconds) (default: 1800.0) | |
| 373 | + deadc: Deadtime correction factor (default: 0.95) | |
| 374 | + emin: Minimum event energy (TeV) (default: 0.1) | |
| 375 | + emax: Maximum event energy (TeV) (default: 100.0) | |
| 376 | + rad: ROI radius used for analysis (deg) (default: 5.0) | |
| 377 | + irf: Instrument response function (default: "South_50h") | |
| 378 | + caldb: Calibration database path (default: "prod2") | |
| 379 | + id: Run identifier (default: "000000") | |
| 372 | 380 | """ |
| 373 | 381 | # Allocate CTA observation |
| 374 | - obs_cta = gammalib.GCTAObservation(instrument) | |
| 382 | + obs_cta = gammalib.GCTAObservation() | |
| 375 | 383 | |
| 376 | 384 | # Set calibration database |
| 377 | 385 | db = gammalib.GCaldb() |
| ... | ... | @@ -397,7 +405,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ |
| 397 | 405 | gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) |
| 398 | 406 | |
| 399 | 407 | # Set energy boundaries |
| 400 | - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ | |
| 408 | + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), | |
| 401 | 409 | gammalib.GEnergy(emax, "TeV")) |
| 402 | 410 | |
| 403 | 411 | # Allocate event list | ... | ... |
src/ctbin/ctbin.cpp
| ... | ... | @@ -52,6 +52,8 @@ const GEnergy g_energy_margin(1.0e-12, "TeV"); |
| 52 | 52 | |
| 53 | 53 | /***********************************************************************//** |
| 54 | 54 | * @brief Void constructor |
| 55 | + * | |
| 56 | + * Constructs an empty ctbin tool. | |
| 55 | 57 | ***************************************************************************/ |
| 56 | 58 | ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) |
| 57 | 59 | { |
| ... | ... | @@ -68,8 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) |
| 68 | 70 | * |
| 69 | 71 | * param[in] obs Observation container. |
| 70 | 72 | * |
| 71 | - * This method creates an instance of the class by copying an existing | |
| 72 | - * observations container. | |
| 73 | + * Constructs a ctbin tool from an observation container. | |
| 73 | 74 | ***************************************************************************/ |
| 74 | 75 | ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) |
| 75 | 76 | { |
| ... | ... | @@ -90,6 +91,8 @@ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) |
| 90 | 91 | * |
| 91 | 92 | * @param[in] argc Number of arguments in command line. |
| 92 | 93 | * @param[in] argv Array of command line arguments. |
| 94 | + * | |
| 95 | + * Constructs ctbin tool from command line arguments. | |
| 93 | 96 | ***************************************************************************/ |
| 94 | 97 | ctbin::ctbin(int argc, char *argv[]) : |
| 95 | 98 | ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv) |
| ... | ... | @@ -106,6 +109,8 @@ ctbin::ctbin(int argc, char *argv[]) : |
| 106 | 109 | * @brief Copy constructor |
| 107 | 110 | * |
| 108 | 111 | * @param[in] app Application. |
| 112 | + * | |
| 113 | + * Constructs ctbin tool from another ctbin instance. | |
| 109 | 114 | ***************************************************************************/ |
| 110 | 115 | ctbin::ctbin(const ctbin& app) : ctool(app) |
| 111 | 116 | { |
| ... | ... | @@ -122,6 +127,8 @@ ctbin::ctbin(const ctbin& app) : ctool(app) |
| 122 | 127 | |
| 123 | 128 | /***********************************************************************//** |
| 124 | 129 | * @brief Destructor |
| 130 | + * | |
| 131 | + * Destructs ctbin tool. | |
| 125 | 132 | ***************************************************************************/ |
| 126 | 133 | ctbin::~ctbin(void) |
| 127 | 134 | { |
| ... | ... | @@ -143,6 +150,9 @@ ctbin::~ctbin(void) |
| 143 | 150 | * @brief Assignment operator |
| 144 | 151 | * |
| 145 | 152 | * @param[in] app Application. |
| 153 | + * @return ctbin tool. | |
| 154 | + * | |
| 155 | + * Assigns ctbin tool. | |
| 146 | 156 | ***************************************************************************/ |
| 147 | 157 | ctbin& ctbin::operator=(const ctbin& app) |
| 148 | 158 | { |
| ... | ... | @@ -175,7 +185,9 @@ ctbin& ctbin::operator=(const ctbin& app) |
| 175 | 185 | ==========================================================================*/ |
| 176 | 186 | |
| 177 | 187 | /***********************************************************************//** |
| 178 | - * @brief Clear instance | |
| 188 | + * @brief Clear ctbin tool | |
| 189 | + * | |
| 190 | + * Clears ctbin tool. | |
| 179 | 191 | ***************************************************************************/ |
| 180 | 192 | void ctbin::clear(void) |
| 181 | 193 | { |
| ... | ... | @@ -195,12 +207,10 @@ void ctbin::clear(void) |
| 195 | 207 | |
| 196 | 208 | |
| 197 | 209 | /***********************************************************************//** |
| 198 | - * @brief Bin the event data | |
| 210 | + * @brief Run the ctbin tool | |
| 199 | 211 | * |
| 200 | - * This method loops over all observations found in the observation conatiner | |
| 201 | - * and bins all events from the event list(s) into counts map(s). Note that | |
| 202 | - * each event list is binned in a separate counts map, hence no summing of | |
| 203 | - * events is done. | |
| 212 | + * Gets the user parameters and loops over all CTA observations in the | |
| 213 | + * observation container to bin the events into a single counts cube. | |
| 204 | 214 | ***************************************************************************/ |
| 205 | 215 | void ctbin::run(void) |
| 206 | 216 | { |
| ... | ... | @@ -212,33 +222,17 @@ void ctbin::run(void) |
| 212 | 222 | // Get task parameters |
| 213 | 223 | get_parameters(); |
| 214 | 224 | |
| 215 | - // Write parameters into logger | |
| 216 | - if (logTerse()) { | |
| 217 | - log_parameters(); | |
| 218 | - log << std::endl; | |
| 219 | - } | |
| 220 | - | |
| 221 | 225 | // Write observation(s) into logger |
| 222 | 226 | if (logTerse()) { |
| 223 | 227 | log << std::endl; |
| 224 | - if (m_obs.size() > 1) { | |
| 225 | - log.header1("Observations"); | |
| 226 | - } | |
| 227 | - else { | |
| 228 | - log.header1("Observation"); | |
| 229 | - } | |
| 228 | + log.header1(gammalib::number("Observation", m_obs.size())); | |
| 230 | 229 | log << m_obs << std::endl; |
| 231 | 230 | } |
| 232 | 231 | |
| 233 | 232 | // Write header |
| 234 | 233 | if (logTerse()) { |
| 235 | 234 | log << std::endl; |
| 236 | - if (m_obs.size() > 1) { | |
| 237 | - log.header1("Bin observations"); | |
| 238 | - } | |
| 239 | - else { | |
| 240 | - log.header1("Bin observation"); | |
| 241 | - } | |
| 235 | + log.header1(gammalib::number("Bin observation", m_obs.size())); | |
| 242 | 236 | } |
| 243 | 237 | |
| 244 | 238 | // Loop over all observations in the container |
| ... | ... | @@ -247,10 +241,10 @@ void ctbin::run(void) |
| 247 | 241 | // Write header for observation |
| 248 | 242 | if (logTerse()) { |
| 249 | 243 | std::string header = m_obs[i]->instrument() + " observation"; |
| 250 | - if (m_obs[i]->name().length() > 1) { | |
| 244 | + if (!m_obs[i]->name().empty()) { | |
| 251 | 245 | header += " \"" + m_obs[i]->name() + "\""; |
| 252 | 246 | } |
| 253 | - if (m_obs[i]->id().length() > 1) { | |
| 247 | + if (!m_obs[i]->id().empty()) { | |
| 254 | 248 | header += " (id=" + m_obs[i]->id() +")"; |
| 255 | 249 | } |
| 256 | 250 | log.header3(header); |
| ... | ... | @@ -293,12 +287,7 @@ void ctbin::run(void) |
| 293 | 287 | // Write observation(s) into logger |
| 294 | 288 | if (logTerse()) { |
| 295 | 289 | log << std::endl; |
| 296 | - if (m_obs.size() > 1) { | |
| 297 | - log.header1("Binned observations"); | |
| 298 | - } | |
| 299 | - else { | |
| 300 | - log.header1("Binned observation"); | |
| 301 | - } | |
| 290 | + log.header1(gammalib::number("Binned observation", m_obs.size())); | |
| 302 | 291 | log << m_obs << std::endl; |
| 303 | 292 | } |
| 304 | 293 | |
| ... | ... | @@ -496,6 +485,12 @@ void ctbin::get_parameters(void) |
| 496 | 485 | m_outcube = (*this)["outcube"].filename(); |
| 497 | 486 | } |
| 498 | 487 | |
| 488 | + // Write parameters into logger | |
| 489 | + if (logTerse()) { | |
| 490 | + log_parameters(); | |
| 491 | + log << std::endl; | |
| 492 | + } | |
| 493 | + | |
| 499 | 494 | // Return |
| 500 | 495 | return; |
| 501 | 496 | } |
| ... | ... | @@ -521,8 +516,8 @@ void ctbin::fill_cube(GCTAObservation* obs) |
| 521 | 516 | const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs->events()); |
| 522 | 517 | if (events == NULL) { |
| 523 | 518 | std::string msg = "CTA Observation does not contain an event " |
| 524 | - "list. Event list information is needed to " | |
| 525 | - "fill the counts map."; | |
| 519 | + "list. An event list is needed to fill the " | |
| 520 | + "counts cube."; | |
| 526 | 521 | throw GException::invalid_value(G_FILL_CUBE, msg); |
| 527 | 522 | } |
| 528 | 523 | |
| ... | ... | @@ -534,7 +529,7 @@ void ctbin::fill_cube(GCTAObservation* obs) |
| 534 | 529 | |
| 535 | 530 | // Check for RoI sanity |
| 536 | 531 | if (!roi.is_valid()) { |
| 537 | - std::string msg = "No RoI information found in input observation " | |
| 532 | + std::string msg = "No valid RoI found in input observation " | |
| 538 | 533 | "\""+obs->name()+"\". Run ctselect to specify " |
| 539 | 534 | "an RoI for this observation before running " |
| 540 | 535 | "ctbin."; |
| ... | ... | @@ -624,10 +619,12 @@ void ctbin::fill_cube(GCTAObservation* obs) |
| 624 | 619 | * @brief Create output observation container. |
| 625 | 620 | * |
| 626 | 621 | * Creates an output observation container that combines all input CTA |
| 627 | - * observation into a single cube-style observation. All non CTA observations | |
| 628 | - * present in the observation container are kept. The method furthermore | |
| 629 | - * conserves any response information in case that a single CTA observation | |
| 630 | - * is provided. This supports the original binned analysis. | |
| 622 | + * observation into a single stacked observation. All non-CTA observations | |
| 623 | + * and all binned CTA observations that were present in the observation | |
| 624 | + * container are append to the observation container so that they can | |
| 625 | + * be used by other tools. The method furthermore conserves any response | |
| 626 | + * information in case that a single CTA observation is provided to support | |
| 627 | + * binned analysis. | |
| 631 | 628 | ***************************************************************************/ |
| 632 | 629 | void ctbin::obs_cube(void) |
| 633 | 630 | { |
| ... | ... | @@ -641,26 +638,37 @@ void ctbin::obs_cube(void) |
| 641 | 638 | GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]); |
| 642 | 639 | if (obs != NULL) { |
| 643 | 640 | |
| 644 | - // Only change the event type if we had an unbinned observation | |
| 641 | + // Change the event type if we had an unbinned observation | |
| 645 | 642 | if (obs->eventtype() == "EventList") { |
| 646 | 643 | |
| 644 | + // Assign cube to the observation | |
| 647 | 645 | obs->events(this->cube()); |
| 648 | - obs->eventfile(""); | |
| 649 | - } // endif: observation was unbinned | |
| 650 | - else { // input observation was binned and skipped | |
| 651 | 646 | |
| 652 | - // Create new and empty cube | |
| 653 | - GCTAEventCube cube = GCTAEventCube(m_cube, m_ebounds, obs->gti()); | |
| 647 | + } | |
| 648 | + | |
| 649 | + // ... otherwise the input observation was binned and hence | |
| 650 | + // skipped. In that case we simply append an empty counts | |
| 651 | + // cube | |
| 652 | + else { | |
| 653 | + | |
| 654 | + // Create empty counts cube | |
| 655 | + GCTAEventCube cube(m_cube, m_ebounds, obs->gti()); | |
| 654 | 656 | |
| 655 | - // Assign empty cube to have a new binned observation | |
| 657 | + // Assign empty cube to observation | |
| 656 | 658 | obs->events(cube); |
| 659 | + | |
| 657 | 660 | } |
| 658 | 661 | |
| 662 | + // Reset file name | |
| 663 | + obs->eventfile(""); | |
| 664 | + | |
| 659 | 665 | } // endif: obervation was valid |
| 660 | 666 | |
| 661 | 667 | } // endif: we only had one observation in the container |
| 662 | 668 | |
| 663 | - // ... otherwise put a single CTA observation in container | |
| 669 | + // ... otherwise put a single CTA observation in container and | |
| 670 | + // append all observations that have not been used in the binning | |
| 671 | + // to the container | |
| 664 | 672 | else { |
| 665 | 673 | |
| 666 | 674 | // Allocate observation container |
| ... | ... | @@ -672,10 +680,25 @@ void ctbin::obs_cube(void) |
| 672 | 680 | // Attach event cube to CTA observation |
| 673 | 681 | obs.events(this->cube()); |
| 674 | 682 | |
| 675 | - // Set map centre as pointing | |
| 676 | - GSkyPixel pixel(0.5*double(m_cube.nx()), 0.5*double(m_cube.ny())); | |
| 677 | - GSkyDir centre = m_cube.pix2dir(pixel); | |
| 678 | - GCTAPointing pointing(centre); | |
| 683 | + // Compute average pointing direction for all CTA event lists | |
| 684 | + double ra = 0.0; | |
| 685 | + double dec = 0.0; | |
| 686 | + double number = 0.0; | |
| 687 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 688 | + GCTAObservation* cta = dynamic_cast<GCTAObservation*>(m_obs[i]); | |
| 689 | + if ((cta != NULL) && (cta->eventtype() == "EventList")) { | |
| 690 | + ra += cta->pointing().dir().ra(); | |
| 691 | + dec += cta->pointing().dir().dec(); | |
| 692 | + number += 1.0; | |
| 693 | + } | |
| 694 | + } | |
| 695 | + if (number > 0.0) { | |
| 696 | + ra /= number; | |
| 697 | + dec /= number; | |
| 698 | + } | |
| 699 | + GSkyDir dir; | |
| 700 | + dir.radec(ra, dec); | |
| 701 | + GCTAPointing pointing(dir); | |
| 679 | 702 | |
| 680 | 703 | // Compute deadtime correction |
| 681 | 704 | double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; |
| ... | ... | @@ -683,8 +706,8 @@ void ctbin::obs_cube(void) |
| 683 | 706 | // Set CTA observation attributes |
| 684 | 707 | obs.pointing(pointing); |
| 685 | 708 | obs.obs_id(0); |
| 686 | - obs.ra_obj(centre.ra_deg()); //!< Dummy | |
| 687 | - obs.dec_obj(centre.dec_deg()); //!< Dummy | |
| 709 | + obs.ra_obj(dir.ra_deg()); //!< Dummy | |
| 710 | + obs.dec_obj(dir.dec_deg()); //!< Dummy | |
| 688 | 711 | obs.ontime(m_ontime); |
| 689 | 712 | obs.livetime(m_livetime); |
| 690 | 713 | obs.deadc(deadc); |
| ... | ... | @@ -701,6 +724,9 @@ void ctbin::obs_cube(void) |
| 701 | 724 | if (obs == NULL) { |
| 702 | 725 | container.append(*m_obs[i]); |
| 703 | 726 | } |
| 727 | + else if (obs->eventtype() != "EventList") { | |
| 728 | + container.append(*m_obs[i]); | |
| 729 | + } | |
| 704 | 730 | } |
| 705 | 731 | |
| 706 | 732 | // Set observation container | ... | ... |
src/ctbin/ctbin.par
| ... | ... | @@ -7,27 +7,27 @@ |
| 7 | 7 | # File: ctbin.par # |
| 8 | 8 | # Version: 1.1.0 # |
| 9 | 9 | # # |
| 10 | -# Author: Jurgen Knodlseder # | |
| 10 | +# Author: Juergen Knoedlseder # | |
| 11 | 11 | # Jurgen.Knodlseder@irap.omp.eu # |
| 12 | 12 | # IRAP # |
| 13 | 13 | # # |
| 14 | -# Purpose: Parameter file for the ctbin application # | |
| 14 | +# Purpose: Parameter file for the ctbin tool # | |
| 15 | 15 | # # |
| 16 | 16 | ############################################################################# |
| 17 | 17 | # |
| 18 | 18 | # General parameters |
| 19 | 19 | #=================== |
| 20 | -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" | |
| 21 | -outcube, f, a, "cntcube.fits",,, "Output counts cube file" | |
| 20 | +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" | |
| 21 | +outcube, f, a, cntcube.fits,,, "Output counts cube file" | |
| 22 | 22 | |
| 23 | 23 | # |
| 24 | 24 | # Energy binning parameters |
| 25 | 25 | #========================== |
| 26 | -ebinalg, s, a, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins" | |
| 26 | +ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins" | |
| 27 | 27 | emin, r, a, 0.1,,, "Start value for first energy bin in TeV" |
| 28 | 28 | emax, r, a, 100.0,,, "Stop value for last energy bin in TeV" |
| 29 | -enumbins, i, a, 20,,, "Number of energy bins" | |
| 30 | -ebinfile, f, a, "NONE",,, "Name of the file containing the energy bin definition" | |
| 29 | +enumbins, i, a, 20,1,200, "Number of energy bins" | |
| 30 | +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" | |
| 31 | 31 | |
| 32 | 32 | # |
| 33 | 33 | # Spatial binning parameters |
| ... | ... | @@ -36,10 +36,10 @@ usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" |
| 36 | 36 | nxpix, i, a, 200,,, "Size of the X axis in pixels" |
| 37 | 37 | nypix, i, a, 200,,, "Size of the Y axis in pixels" |
| 38 | 38 | binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" |
| 39 | -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" | |
| 39 | +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" | |
| 40 | 40 | xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" |
| 41 | 41 | yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" |
| 42 | -proj, s, a, "CAR",AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" | |
| 42 | +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" | |
| 43 | 43 | |
| 44 | 44 | # |
| 45 | 45 | # Standard parameters |
| ... | ... | @@ -48,5 +48,5 @@ publish, b, h, no,,, "Publish counts cube on VO Hub?" |
| 48 | 48 | chatter, i, h, 2,0,4, "Chattiness of output" |
| 49 | 49 | clobber, b, h, yes,,, "Overwrite existing output files with new output files?" |
| 50 | 50 | debug, b, h, no,,, "Debugging mode activated" |
| 51 | -mode, s, h, "ql",,, "Mode of automatic parameters" | |
| 52 | -logfile, f, h, "ctbin.log",,, "Log filename" | |
| 51 | +mode, s, h, ql,,, "Mode of automatic parameters" | |
| 52 | +logfile, f, h, ctbin.log,,, "Log filename" | ... | ... |
test/test_ctbin.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# This scripts performs unit tests for the ctbin tool. | |
| 3 | +# This scripts performs unit tests for the ctbin module. | |
| 4 | 4 | # |
| 5 | 5 | # Copyright (C) 2014-2016 Juergen Knoedlseder |
| 6 | 6 | # |
| ... | ... | @@ -27,8 +27,9 @@ import ctools |
| 27 | 27 | # ========================= # |
| 28 | 28 | class Test(gammalib.GPythonTestSuite): |
| 29 | 29 | """ |
| 30 | - Test class for ctbin tool. | |
| 30 | + Test class for ctbin module. | |
| 31 | 31 | """ |
| 32 | + | |
| 32 | 33 | # Constructor |
| 33 | 34 | def __init__(self): |
| 34 | 35 | """ |
| ... | ... | @@ -78,14 +79,8 @@ class Test(gammalib.GPythonTestSuite): |
| 78 | 79 | bin["xref"] = 83.63 |
| 79 | 80 | bin["yref"] = 22.01 |
| 80 | 81 | |
| 81 | - # Run tool | |
| 82 | - self.test_try("Run ctbin") | |
| 83 | - try: | |
| 84 | - bin.run() | |
| 85 | - self.test_try_success() | |
| 86 | - except: | |
| 87 | - msg = "Exception occured in ctbin." | |
| 88 | - self.test_try_failure(msg) | |
| 82 | + # Run ctbin tool | |
| 83 | + bin.run() | |
| 89 | 84 | |
| 90 | 85 | # Check content of observation and cube |
| 91 | 86 | self._test_observation(bin, 5542) |
| ... | ... | @@ -98,7 +93,7 @@ class Test(gammalib.GPythonTestSuite): |
| 98 | 93 | self._test_observation(cpy_bin, 5542) |
| 99 | 94 | self._test_cube(cpy_bin.cube(), 5542) |
| 100 | 95 | |
| 101 | - # Run copy of ctbin again | |
| 96 | + # Run copy of ctbin tool again | |
| 102 | 97 | cpy_bin.run() |
| 103 | 98 | |
| 104 | 99 | # Check content of observation and cube. We expect now an empty |
| ... | ... | @@ -128,46 +123,116 @@ class Test(gammalib.GPythonTestSuite): |
| 128 | 123 | evt = gammalib.GCTAEventCube("cntmap.fits") |
| 129 | 124 | self._test_cube(evt, 5542) |
| 130 | 125 | |
| 126 | + # Prepare observation container for stacked analysis | |
| 127 | + cta = gammalib.GCTAObservation(self.events_name) | |
| 128 | + obs = gammalib.GObservations() | |
| 129 | + cta.id("0001") | |
| 130 | + obs.append(cta) | |
| 131 | + cta.id("0002") | |
| 132 | + obs.append(cta) | |
| 133 | + cta.id("0003") | |
| 134 | + obs.append(cta) | |
| 135 | + | |
| 136 | + # Set-up ctbin using an observation container | |
| 137 | + bin = ctools.ctbin(obs) | |
| 138 | + bin["outcube"] = "cntmap.fits" | |
| 139 | + bin["ebinalg"] = "LOG" | |
| 140 | + bin["emin"] = 0.1 | |
| 141 | + bin["emax"] = 100.0 | |
| 142 | + bin["enumbins"] = 20 | |
| 143 | + bin["nxpix"] = 200 | |
| 144 | + bin["nypix"] = 200 | |
| 145 | + bin["binsz"] = 0.02 | |
| 146 | + bin["coordsys"] = "CEL" | |
| 147 | + bin["proj"] = "CAR" | |
| 148 | + bin["xref"] = 83.63 | |
| 149 | + bin["yref"] = 22.01 | |
| 150 | + | |
| 151 | + # Run ctbin tool | |
| 152 | + bin.run() | |
| 153 | + | |
| 154 | + # Check content of observation and cube (need multiplier=3 since | |
| 155 | + # three identical observations have been appended) | |
| 156 | + self._test_observation(bin, 5542, multiplier=3) | |
| 157 | + self._test_cube(bin.cube(), 5542, multiplier=3) | |
| 158 | + | |
| 159 | + # Set-up ctbin with invalid event file | |
| 160 | + bin = ctools.ctbin() | |
| 161 | + bin["inobs"] = "event_file_that_does_not_exist.fits" | |
| 162 | + bin["outcube"] = "cntmap.fits" | |
| 163 | + bin["ebinalg"] = "LOG" | |
| 164 | + bin["emin"] = 0.1 | |
| 165 | + bin["emax"] = 100.0 | |
| 166 | + bin["enumbins"] = 20 | |
| 167 | + bin["nxpix"] = 200 | |
| 168 | + bin["nypix"] = 200 | |
| 169 | + bin["binsz"] = 0.02 | |
| 170 | + bin["coordsys"] = "CEL" | |
| 171 | + bin["proj"] = "CAR" | |
| 172 | + bin["xref"] = 83.63 | |
| 173 | + bin["yref"] = 22.01 | |
| 174 | + | |
| 175 | + # Run ctbin tool | |
| 176 | + self.test_try("Run ctbin with invalid event file") | |
| 177 | + try: | |
| 178 | + bin.run() | |
| 179 | + self.test_try_failure() | |
| 180 | + except: | |
| 181 | + self.test_try_success() | |
| 182 | + | |
| 131 | 183 | # Return |
| 132 | 184 | return |
| 133 | 185 | |
| 134 | 186 | # Check observation |
| 135 | - def _test_observation(self, ctbin, nevents): | |
| 187 | + def _test_observation(self, ctbin, nevents, multiplier=1): | |
| 136 | 188 | """ |
| 137 | 189 | Test content of an observation. |
| 138 | 190 | |
| 139 | 191 | Args: |
| 140 | 192 | ctbin: ctbin instance |
| 141 | 193 | nevents: Expected number of events |
| 194 | + | |
| 195 | + Kwargs: | |
| 196 | + multiplier: Observation multiplier | |
| 142 | 197 | """ |
| 143 | 198 | # Test observation container |
| 144 | 199 | obs = gammalib.GCTAObservation(ctbin.obs()[0]) |
| 145 | 200 | pnt = obs.pointing() |
| 146 | - self.test_value(ctbin.obs().size(), 1, "There is one observation") | |
| 147 | - self.test_assert(obs.instrument() == "CTA", "Observation is CTA observation") | |
| 148 | - self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is 1800 sec") | |
| 149 | - self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is 1710 sec") | |
| 150 | - self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "Pointing Right Ascension is 83.63 deg") | |
| 151 | - self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "Pointing Declination is 22.01 deg") | |
| 201 | + self.test_value(ctbin.obs().size(), 1, | |
| 202 | + "There is one observation") | |
| 203 | + self.test_assert(obs.instrument() == "CTA", | |
| 204 | + "Observation is CTA observation") | |
| 205 | + self.test_value(obs.ontime(), 1800.0*multiplier, 1.0e-6, | |
| 206 | + "Ontime is 1800 sec") | |
| 207 | + self.test_value(obs.livetime(), 1710.0*multiplier, 1.0e-6, | |
| 208 | + "Livetime is 1710 sec") | |
| 209 | + self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, | |
| 210 | + "Pointing Right Ascension is 83.63 deg") | |
| 211 | + self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, | |
| 212 | + "Pointing Declination is 22.01 deg") | |
| 152 | 213 | |
| 153 | 214 | # Test event cube |
| 154 | - self._test_cube(obs.events(), nevents) | |
| 215 | + self._test_cube(obs.events(), nevents, multiplier=multiplier) | |
| 155 | 216 | |
| 156 | 217 | # Return |
| 157 | 218 | return |
| 158 | 219 | |
| 159 | 220 | # Check event cube |
| 160 | - def _test_cube(self, cube, nevents): | |
| 221 | + def _test_cube(self, cube, nevents, multiplier=1): | |
| 161 | 222 | """ |
| 162 | 223 | Test content of event cube." |
| 163 | 224 | |
| 164 | 225 | Args: |
| 165 | 226 | cube: Event cube |
| 166 | 227 | nevents: Expected number of events. |
| 228 | + | |
| 229 | + Kwargs: | |
| 230 | + multiplier: Event number multiplier | |
| 167 | 231 | """ |
| 168 | 232 | # Test event cube |
| 169 | 233 | self.test_value(cube.size(), 800000, "800000 event bins") |
| 170 | - self.test_value(cube.number(), nevents, str(nevents)+" events") | |
| 234 | + self.test_value(cube.number(), nevents*multiplier, | |
| 235 | + str(nevents*multiplier)+" events") | |
| 171 | 236 | |
| 172 | 237 | # Return |
| 173 | 238 | return | ... | ... |