diff --git a/ChangeLog b/ChangeLog index 00d517d..a99b94a 100755 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,7 @@ ======================== Apply event bin weighting in ctmodel for stacked computation (#2721) + Enable ctbin to produce counts cubes for a joint binned analysis (#2714) Add Jupyter notebooks for H.E.S.S. DR1 (#2712) Add code generator (#2709) Add csbkgmodel script (#2694) diff --git a/NEWS b/NEWS index b131636..86388de 100755 --- a/NEWS +++ b/NEWS @@ -61,6 +61,14 @@ of the filename is now done by the GApplicationPar class. ctbin - CTA event binning ------------------------- +Added the option to produce an individual counts cube for each event list in +the input observation definition XML file. The option is enabled by setting +the new hidden parameter "stack" to "no". The cube() method now has an optional +index argument, providing access to the i'th counts cube that was produced. +The cubes() method specifies the number of counts cubes that were produced. +If no counts cube was produced, the cube() method now throws an exception +(#2714). + ctbin now accouts for varying energy thresholds when computing the weight array (#2685). diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 4a67755..4178fbf 100755 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -164,6 +164,7 @@ testscripts = $(srcdir)/tests/__init__.py \ # Test data to be distributed in $(pyexecdir)/cscripts/tests/data testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/obs_unbinned.xml \ + $(top_srcdir)/test/data/obs_unbinned_two.xml \ $(top_srcdir)/test/data/obs_stacked.xml \ $(top_srcdir)/test/data/obs_stacked_two.xml \ $(top_srcdir)/test/data/crab_events.fits \ diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 5a8927a..80ef31a 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -43,6 +43,9 @@ Bug fixes Improvements ------------ +* [`2714 `_] - + Enable :ref:`ctbin` to produce counts cubes for a joint binned analysis. + Rename ``outcube`` to ``outobs`` parameter. * [`2712 `_] - Add Jupyter notebooks for H.E.S.S. DR1 * [`2685 `_] - diff --git a/doc/source/users/reference_manual/ctbin.rst b/doc/source/users/reference_manual/ctbin.rst index d10e07d..2c523e2 100755 --- a/doc/source/users/reference_manual/ctbin.rst +++ b/doc/source/users/reference_manual/ctbin.rst @@ -14,17 +14,27 @@ a 3-dimensional data cube spanned by Right Ascension or Galactic longitude, Declination or Galactic latitude, and energy. The energy binning may be either linear, logarithmic, or custom defined using an input file. The events are either taken from a single event list file or from the event lists that are -specified in an observation definition file. In case that multiple event -lists are given in an observation definition file, the tool will loop over -all event lists and fill all events into a single counts cube. +specified in an observation definition file. -:ref:`ctbin` generates a counts cube FITS file comprising three extensions. The +In case that multiple event lists are given in an observation definition file +there are two options. By default, the tool will loop over all event lists and +fill all events into a single counts cube. If ``stack=no`` is specified, the +tool will bin each event list into a separate counts cube. + +:ref:`ctbin` generates counts cube FITS file(s) comprising three extensions. The primary extension contains a 3-dimensional image that contains the counts cube values. The next extension named ``EBOUNDS`` contains a binary table that defines the energy boundaries of the counts cube. The last extension named ``GTI`` contains a binary table that defines the Good Time Intervals of all event lists that have been filled into the counts cube. +By default, the tool will write out a single, stacked, counts cube. If ``stack=no`` +is specified, the tool will write all counts cubes into a separate FITS file, +prefixing the input event list file name with the string that is specified by +the hidden ``prefix`` parameter. In addition, the tool will create an observation +definition XML file that lists all produced counts cubes, and that can be used +as input for any further analysis step. + General parameters ------------------ @@ -32,8 +42,18 @@ General parameters ``inobs [file]`` Input event list or observation definition XML file. -``outcube [file]`` - Output counts cube file. +``outobs [file]`` + Output counts cube file or observation definition XML file. + +``(stack = yes) [boolean]`` + Specify whether stacking of the events in a single counts cube is requested. + By default this hidden parameter is set to yes, and the tool produces a + single counts cube FITS file on output. If ``no`` is specified, the tool + will produce a counts cube for each event list that is found in the input + observation, and will also write an observation definition XML file. + +``(prefix = cntcube_) [string]`` + Prefix for output counts cube in observation definition XML file. ``ebinalg [string]`` Algorithm for defining energy bins. For ``FILE``, the energy bins are defined diff --git a/doc/source/users/reference_manual/ctselect.rst b/doc/source/users/reference_manual/ctselect.rst index 01e4e71..f85f487 100755 --- a/doc/source/users/reference_manual/ctselect.rst +++ b/doc/source/users/reference_manual/ctselect.rst @@ -52,7 +52,7 @@ General parameters Output event list or observation definition XML file. ``(prefix = "selected_") [string]`` - Prefix for output event lists in observation definition file. + Prefix for output event lists in observation definition XML file. ``(usepnt = yes) [boolean]`` Use pointing instead of RA/DEC parameters? diff --git a/pyext/ctbin.i b/pyext/ctbin.i index e6c564c..e6ac9d8 100755 --- a/pyext/ctbin.i +++ b/pyext/ctbin.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -46,9 +46,10 @@ public: // Methods void clear(void); void run(void); + int cubes(void) const; + const GCTAEventCube& cube(const int& index = 0) const; void save(void); void publish(const std::string& name = ""); - const GCTAEventCube& cube(void) const; }; diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index dfb5c3a..56ec036 100755 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -38,6 +38,7 @@ #endif /* __ Method name definitions ____________________________________________ */ +#define G_CUBE "ctbin::cube(int&)" #define G_GET_PARAMETERS "ctbin::get_parameters()" #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)" #define G_SET_WEIGHTS "ctbin::set_weights(GCTAObservation*)" @@ -234,6 +235,9 @@ void ctbin::run(void) // Get task parameters get_parameters(); + // Initialise sky direction cache for stacking + init_sky_dir_cache(); + // Write input observation container into logger log_observations(NORMAL, m_obs, "Input observation"); @@ -253,64 +257,171 @@ void ctbin::run(void) std::string msg = " Including unbinned "+obs->instrument()+ " observation"; log_string(NORMAL, msg); + + } // endfor: looped over all unbinned observations + + // Set number of relevant observations + int nobs = obs_list.size(); + + // Initialise arrays + if (nobs > 0) { + m_counts = std::vector(nobs); + m_weights = std::vector(nobs); + m_cubes = std::vector(nobs); + } + else { + m_counts.clear(); + m_weights.clear(); + m_cubes.clear(); } // Write header into logger - log_header1(TERSE, gammalib::number("Bin observation", m_obs.size())); + log_header1(TERSE, gammalib::number("Bin observation", nobs)); // Loop over all unbinned CTA observations in the container #pragma omp parallel for - for (int i = 0; i < obs_list.size(); ++i) { + for (int i = 0; i < nobs; ++i) { // Get pointer to observation GCTAObservation* obs = obs_list[i]; // Fill the cube - fill_cube(obs); + GSkyMap counts = fill_cube(obs); // Set the counts cube weights - set_weights(obs); + GSkyMap weights = set_weights(obs); // Dispose events to free memory obs->dispose_events(); + // Post process cubes + #pragma omp critical(ctbin_fill_cube) + { + // If stacking of cubes is requested then add up counts cube + // and set weights in the first slot of the m_counts and + // m_weights arrays + if (m_stack) { + if (i == 0) { + m_counts[0] = counts; + m_weights[0] = weights; + } + else { + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { + for (int pixel = 0; pixel < counts.npix(); ++pixel) { + m_counts[0](pixel, iebin) += counts(pixel, iebin); + if (weights(pixel, iebin) == 1.0) { + m_weights[0](pixel, iebin) = 1.0; + } + } + } + } + } + + // ... otherwise store counts and weights + else { + m_counts[i] = counts; + m_weights[i] = weights; + } + + } // end: omp critical section + } // endfor: looped over observations - // Compute total livetime in all observations - double livetime = 0.0; - for (int i = 0; i < obs_list.size(); ++i) { - livetime += obs_list[i]->livetime(); - } + // Combine GTIs of all observations and compute total ontime and livetime + for (int i = 0; i < nobs; ++i) { + + // Get pointer to observation + GCTAObservation* obs = obs_list[i]; - // Compute livetime fraction per energy bin so that varying energy - // thresholds are correctly taken into account in the weighting - // array. Computational speed is not optimum, we have to see later - // whether this is an issue. - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - double livetime_ebin = 0.0; - for (int i = 0; i < obs_list.size(); ++i) { - const GCTAEventList* events = dynamic_cast(obs_list[i]->events()); - std::vector usage = cube_layer_usage(m_ebounds, events->ebounds()); - if (usage[iebin]) { - livetime_ebin += obs_list[i]->livetime(); + // Get Good Time Intervals for observation + GGti gti = obs->gti(); + + // If resulting GTI is not filled then copy over the reference from + // the observation and use it for the resulting GTI + if (m_gti.is_empty()) { + m_gti.reference(gti.reference()); + } + + // Append GTIs + m_gti.extend(gti); + + // Update ontime and livetime + m_ontime += obs->ontime(); + m_livetime += obs->livetime(); + + } // endfor: looped over all observations + + // Post process data for a stacked cube + if (m_stack && nobs > 0) { + + // Compute livetime fraction per energy bin so that varying energy + // thresholds are correctly taken into account in the weighting + // array. Computational speed is not optimum, we have to see later + // whether this is an issue. + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { + double livetime_ebin = 0.0; + for (int i = 0; i < obs_list.size(); ++i) { + std::vector usage = + cube_layer_usage(m_ebounds, obs_list[i]->ebounds()); + if (usage[iebin]) { + livetime_ebin += obs_list[i]->livetime(); + } } - } - if (livetime > 0.0) { - double fraction = livetime_ebin / livetime; - for (int pixel = 0; pixel < m_counts.npix(); ++pixel) { - if (m_weights(pixel, iebin) > 0.0) { - m_weights(pixel, iebin) *= fraction; + if (m_livetime > 0.0) { + double fraction = livetime_ebin / m_livetime; + for (int pixel = 0; pixel < m_counts[0].npix(); ++pixel) { + if (m_weights[0](pixel, iebin) > 0.0) { + m_weights[0](pixel, iebin) *= fraction; + } } } } - } - // Build event cube (needs to come before obs_cube() since this method - // relies on correct setting of m_cube) - m_cube = GCTAEventCube(m_counts, m_weights, m_ebounds, m_gti); + // Set a single cube in the observation container + obs_cube_stacked(); + + } // endif: stacked cube post processing + + // Post process data for non-stacked cubes + else { + + // Loop over all relevant observations + for (int i = 0; i < nobs; ++i) { + + // Get pointer to observation + GCTAObservation* obs = obs_list[i]; + + // Build counts cube + m_cubes[i] = GCTAEventCube(m_counts[i], m_weights[i], m_ebounds, + obs->gti()); + + // Attach event cube in the correct slot of the observation + // container + obs->events(m_cubes[i]); + + // Construct unique filename for possible saving of the counts + // cube into FITS files. Each observation has a unique pair of + // instrument() and id() attributes, hence we construct a lower + // case file name from these attributes, making sure that all + // whitespaces are converted into "_" characters. Note that the + // id() attribute is allowed to be blank, which can only happen + // for a single observation. + std::string inst = gammalib::tolower(obs->instrument()); + std::string id = gammalib::tolower(obs->id()); + inst = gammalib::replace_segment(inst, " ", "_"); + id = gammalib::replace_segment(id, " ", "_"); + std::string filename = m_prefix + inst; + if (!id.empty()) { + filename += "_" + id; + } + filename += ".fits"; - // Set a single cube in the observation container - obs_cube(); + // Set filename of observation + obs->eventfile(filename); + + } // endfor: looped over relevant observations + + } // endelse: non-stacked cube post processing // Write resulting observation container into logger log_observations(NORMAL, m_obs, "Binned observation"); @@ -326,37 +437,106 @@ void ctbin::run(void) /***********************************************************************//** + * @brief Return event cube at index + * + * @param[in] index Cube index (0,...,cubes()-1). + * @return Reference to event cube. + * + * @exception GException::out_of_range + * Event cube index out of range. + * + * Returns a reference to the event cube at the given index. + ***************************************************************************/ +const GCTAEventCube& ctbin::cube(const int& index) const +{ + // Compile option: raise an exception if index is out of range + if (index < 0 || index >= cubes()) { + throw GException::out_of_range(G_CUBE, "Cube index", index, cubes()); + } + + // Return cube + return m_cubes[index]; +} + + +/***********************************************************************//** * @brief Save counts cube * - * Saves the counts cube into a FITS file. + * Saves the counts cube. If the counts cube was stacked single FITS file + * will be written. Otherwise each generated counts cube will be written + * into a FITS file, and an observation definition XML file will be written + * out that contains an observation for each counts cube. ***************************************************************************/ void ctbin::save(void) { // Write header log_header1(TERSE, "Save counts cube"); - // Save only if filename is valid and if there is at least one - // observation - if ((*this)["outcube"].is_valid() && m_obs.size() > 0) { + // If counts cube was stacked then save first observation + if (m_stack) { - // Get counts cube filename - m_outcube = (*this)["outcube"].filename(); + // Save only if filename is valid and if there is at least one + // observation + if ((*this)["outobs"].is_valid() && m_obs.size() > 0) { - // Get CTA observation from observation container - GCTAObservation* obs = dynamic_cast(m_obs[0]); + // Get counts cube filename + GFilename outobs = (*this)["outobs"].filename(); - // Save only if observation is valid - if (obs != NULL) { + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(m_obs[0]); + + // Save only if observation is valid + if (obs != NULL) { + + // Log counts cube file name + log_value(NORMAL, "Counts cube file", outobs.url()); - // Log counts cube file name - log_value(NORMAL, "Counts cube file", m_outcube.url()); + // Save cube + obs->save(outobs, clobber()); - // Save cube - obs->save(m_outcube, clobber()); + } // endif: observation was valid + + } // endif: outobs file was valid + + } // endif: counts cube was stacked + + // ... otherwise write out all relevant counts cubes and the observation + // definition XML file + else { + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(m_obs[i]); - } // endif: observation was valid + // Continue only if observation is valid + if (obs != NULL) { - } // endif: outcube file was valid + // Retrieve filename + std::string filename = obs->eventfile(); + + // Save counts cube + obs->save(filename, clobber()); + + } // endif: observation was valid + + } // endfor: looped over all observations + + // Write XML file + if ((*this)["outobs"].is_valid()) { + + // Get observation definition XML filename + std::string outobs = (*this)["outobs"].filename(); + + // Issue warning if output filename has no .xml suffix + log_string(TERSE, warn_xml_suffix(outobs)); + + // Save XML file + m_obs.save(outobs); + } + + } // endelse: counts cube were not stacked // Return return; @@ -382,8 +562,10 @@ void ctbin::publish(const std::string& name) // Write counts cube name into logger log_value(NORMAL, "Counts cube name", user_name); - // Publish counts cube - m_counts.publish(user_name); + // Publish first counts cube if it exists + if (cubes() > 0) { + m_counts[0].publish(user_name); + } // Return return; @@ -402,17 +584,18 @@ void ctbin::publish(const std::string& name) void ctbin::init_members(void) { // Initialise members - m_outcube.clear(); m_usepnt = false; + m_stack = true; + m_prefix.clear(); m_publish = false; m_chatter = static_cast(2); // Initialise protected members + m_cubes.clear(); m_counts.clear(); m_weights.clear(); m_ebounds.clear(); m_gti.clear(); - m_cube.clear(); m_ontime = 0.0; m_livetime = 0.0; @@ -432,17 +615,18 @@ void ctbin::init_members(void) void ctbin::copy_members(const ctbin& app) { // Copy attributes - m_outcube = app.m_outcube; m_usepnt = app.m_usepnt; + m_stack = app.m_stack; + m_prefix = app.m_prefix; m_publish = app.m_publish; m_chatter = app.m_chatter; // Copy protected members + m_cubes = app.m_cubes; m_counts = app.m_counts; m_weights = app.m_weights; m_ebounds = app.m_ebounds; m_gti = app.m_gti; - m_cube = app.m_cube; m_ontime = app.m_ontime; m_livetime = app.m_livetime; @@ -479,32 +663,23 @@ void ctbin::get_parameters(void) // information and do not accept counts cubes. setup_observations(m_obs, false, true, false); - // Create an event cube based on task parameters - GCTAEventCube cube = create_cube(m_obs); - - // Get the skymap from the cube and initialise all counts cube bins and - // weights to zero - m_counts = cube.counts(); - m_counts = 0.0; - m_weights = m_counts; + // Create an event cube to query task parameters + GCTAEventCube cube = this->ctool::create_cube(m_obs); // Get energy boundaries m_ebounds = cube.ebounds(); + // Get stack and prefix parameters + m_stack = (*this)["stack"].boolean(); + m_prefix = (*this)["prefix"].string(); + // Get remaining parameters m_publish = (*this)["publish"].boolean(); m_chatter = static_cast((*this)["chatter"].integer()); // If needed later, query output filename now if (read_ahead()) { - (*this)["outcube"].query(); - } - - // Cache sky directions in m_counts - m_dirs.resize(m_counts.npix()); - for (int pix = 0; pix < m_counts.npix(); ++pix) { - m_dirs[pix] = m_counts.inx2dir(pix); - m_dirs[pix].radec(m_dirs[pix].ra(), m_dirs[pix].dec()); + (*this)["outobs"].query(); } // Set number of OpenMP threads @@ -524,6 +699,95 @@ void ctbin::get_parameters(void) /***********************************************************************//** + * @brief Initialise sky direction cache for cube stack + * + * Initialises a cache that holds the sky directions of each counts cube + * pixel for a stacked analysis. This will speed up the computations later. + * + * If no stacking is requested, the method does nothing. + ***************************************************************************/ +void ctbin::init_sky_dir_cache(void) +{ + // Continue only is stacking is requested + if (m_stack) { + + // Create a sky map from task parameters + GSkyMap map = this->ctool::create_map(m_obs); + + // Initialise sky direction + m_dirs.resize(map.npix()); + + // Compute sky direction for each pixel + for (int pixel = 0; pixel < map.npix(); ++pixel) { + m_dirs[pixel] = map.inx2dir(pixel); + } + + } // endif: stacking was requested + + // Return + return; +} + + +/***********************************************************************//** + * @brief Create counts cube sky map for current observation + * + * @param[in] obs CTA observation. + * + * Creates a counts cube sky map for current observations. + ***************************************************************************/ +GSkyMap ctbin::create_cube(const GCTAObservation* obs) +{ + // Read coordinate system and projection + std::string coordsys = (*this)["coordsys"].string(); + std::string proj = (*this)["proj"].string(); + + // Read sky map reference + double xref = 0.0; + double yref = 0.0; + bool usepnt = (*this)["usepnt"].boolean(); + if (!usepnt) { + xref = (*this)["xref"].real(); + yref = (*this)["yref"].real(); + } + + // Read sky map pixel size and number of bins + double binsz = (*this)["binsz"].real(); + int nxpix = (*this)["nxpix"].integer(); + int nypix = (*this)["nypix"].integer(); + + // Read energy boundaries + GEbounds ebounds = create_ebounds(); + + // If requested, get pointing from observations + if (usepnt) { + + // Get pointing direction for observation + GSkyDir pnt = obs->pointing().dir(); + + // Set xref/yref based on the coordinate system + if (gammalib::toupper(coordsys) == "GAL") { + xref = pnt.l_deg(); + yref = pnt.b_deg(); + } + else { + xref = pnt.ra_deg(); + yref = pnt.dec_deg(); + } + + } // endif: got pointing from observation + + // Initialise counts cube + GSkyMap cube = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz, + nxpix, nypix, ebounds.size()); + + + // Return cube + return cube; +} + + +/***********************************************************************//** * @brief Fill events into counts cube * * @param[in] obs CTA observation. @@ -533,10 +797,10 @@ void ctbin::get_parameters(void) * * Fills the events from an event list into the counts cube. ***************************************************************************/ -void ctbin::fill_cube(GCTAObservation* obs) +GSkyMap ctbin::fill_cube(const GCTAObservation* obs) { - // Make sure that the observation holds a CTA event list. If this - // is not the case then throw an exception. + // Make sure that the observation holds a CTA event list. If this is + // not the case then throw an exception. const GCTAEventList* events = dynamic_cast (obs->events()); if (events == NULL) { @@ -547,7 +811,7 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Get the RoI - const GCTARoi& roi = events->roi(); + const GCTARoi& roi = obs->roi(); // Check for RoI sanity if (!roi.is_valid()) { @@ -559,7 +823,7 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Get counts cube usage flags - std::vector usage = cube_layer_usage(m_ebounds, events->ebounds()); + std::vector usage = cube_layer_usage(m_ebounds, obs->ebounds()); // Initialise binning statistics int num_outside_roi = 0; @@ -568,6 +832,13 @@ void ctbin::fill_cube(GCTAObservation* obs) int num_outside_ebds = 0; int num_in_map = 0; + // Create counts cube + GSkyMap counts = create_cube(obs); + + // Extract counts cube boundaries + double pixel_x_max = counts.nx() - 0.5; + double pixel_y_max = counts.ny() - 0.5; + // Fill counts sky map for (int i = 0; i < events->size(); ++i) { @@ -587,7 +858,7 @@ void ctbin::fill_cube(GCTAObservation* obs) // Determine sky pixel GSkyPixel pixel; try { - pixel = m_counts.dir2pix(dir); + pixel = counts.dir2pix(dir); } catch (std::exception &e) { num_invalid_wcs++; @@ -596,8 +867,8 @@ void ctbin::fill_cube(GCTAObservation* obs) // Skip event if corresponding counts cube pixel is outside the // counts cube map range - if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx()-0.5) || - pixel.y() < -0.5 || pixel.y() > (m_counts.ny()-0.5)) { + if (pixel.x() < -0.5 || pixel.x() > pixel_x_max || + pixel.y() < -0.5 || pixel.y() > pixel_y_max) { num_outside_map++; continue; } @@ -614,31 +885,13 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Fill event in skymap - #pragma omp critical(ctbin_fill_cube) - m_counts(pixel, iebin) += 1.0; + counts(pixel, iebin) += 1.0; // Increment number of maps num_in_map++; } // endfor: looped over all events - // Update time information - #pragma omp critical(ctbin_fill_cube) - { - // If GTI is empty then set its time reference from the observation. - // From then on we keep that time reference - if (m_gti.is_empty()) { - m_gti.reference(events->gti().reference()); - } - - // Append GTIs - m_gti.extend(events->gti()); - - // Update ontime and livetime - m_ontime += obs->ontime(); - m_livetime += obs->livetime(); - } - // Log filling results #pragma omp critical(ctbin_fill_cube) { @@ -652,7 +905,7 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Return - return; + return counts; } @@ -667,21 +920,10 @@ void ctbin::fill_cube(GCTAObservation* obs) * Sets the counts cube weights for all bins that are considered for the * specific observation to unity. ***************************************************************************/ -void ctbin::set_weights(GCTAObservation* obs) +GSkyMap ctbin::set_weights(const GCTAObservation* obs) { - // Make sure that the observation holds a CTA event list. If this - // is not the case then throw an exception. - const GCTAEventList* events = dynamic_cast - (obs->events()); - if (events == NULL) { - std::string msg = "CTA Observation does not contain an event " - "list. An event list is needed to fill the " - "counts cube."; - throw GException::invalid_value(G_SET_WEIGHTS, msg); - } - // Get the RoI - const GCTARoi& roi = events->roi(); + const GCTARoi& roi = obs->roi(); // Check for RoI sanity if (!roi.is_valid()) { @@ -692,15 +934,32 @@ void ctbin::set_weights(GCTAObservation* obs) throw GException::invalid_value(G_SET_WEIGHTS, msg); } + // Get the RoI centre and radius in radians + GSkyDir roi_centre = roi.centre().dir(); + double roi_radius = roi.radius() * gammalib::deg2rad; + // Get counts cube usage flags - std::vector usage = cube_layer_usage(m_ebounds, events->ebounds()); + std::vector usage = cube_layer_usage(m_ebounds, obs->ebounds()); + + // Create weights cube + GSkyMap weights = create_cube(obs); + + // Set sky direction cache availability flag + bool cache = (m_stack && m_dirs.size() == weights.npix()); // Loop over all pixels in counts cube - for (int pixel = 0; pixel < m_counts.npix(); ++pixel) { + for (int pixel = 0; pixel < weights.npix(); ++pixel) { // Skip pixel if it is outside the RoI - if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) { - continue; + if (cache) { + if (m_dirs[pixel].dist(roi_centre) > roi_radius) { + continue; + } + } + else { + if (roi_centre.dist(weights.inx2dir(pixel)) > roi_radius) { + continue; + } } // Loop over all energy layers of counts cube @@ -712,15 +971,14 @@ void ctbin::set_weights(GCTAObservation* obs) } // Signal that bin was filled - #pragma omp critical(ctbin_set_weights) - m_weights(pixel, iebin) = 1.0; + weights(pixel, iebin) = 1.0; } // endfor: looped over energy layers of counts cube } // endfor: looped over pixels of counts cube // Return - return; + return weights; } @@ -735,8 +993,11 @@ void ctbin::set_weights(GCTAObservation* obs) * information in case that a single CTA observation is provided to support * binned analysis. ***************************************************************************/ -void ctbin::obs_cube(void) +void ctbin::obs_cube_stacked(void) { + // Set event cube in first slot + m_cubes[0] = GCTAEventCube(m_counts[0], m_weights[0], m_ebounds, m_gti); + // If we have only a single CTA observation in the container, then // keep that observation and just attach the event cube to it. Reset // the filename, otherwise we still will have the old event filename @@ -751,7 +1012,7 @@ void ctbin::obs_cube(void) if (obs->eventtype() == "EventList") { // Assign cube to the observation - obs->events(this->cube()); + obs->events(m_cubes[0]); } @@ -761,7 +1022,8 @@ void ctbin::obs_cube(void) else { // Create empty counts cube - GCTAEventCube cube(m_counts, m_weights, m_ebounds, obs->gti()); + GCTAEventCube cube(m_counts[0], m_weights[0], m_ebounds, + obs->gti()); // Assign empty cube to observation obs->events(cube); @@ -787,7 +1049,7 @@ void ctbin::obs_cube(void) GCTAObservation obs; // Attach event cube to CTA observation - obs.events(this->cube()); + obs.events(m_cubes[0]); // Compute average pointing direction for all CTA event lists double ra = 0.0; diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index 21cd9ab..d3bf9b6 100755 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -39,7 +39,7 @@ * * @brief Event binning tool * - * This class bins event list(s) into a single counts cube. The class can + * This class bins event list(s) into counts cubes. The class can * operate on predefined observation containers, on individual event list * FITS files, and on observation definition XML files. * @@ -47,6 +47,13 @@ * XML definition file, the class will merge these events into a single * counts cube. * + * If the hidden parameter stack=no is used, one counts cube is generated for + * each individual observation. These cubes are saved to a path that can be + * specified by the hidden parameter prefix, followed by the corresponding + * observation id. An observation definition XML file containing the paths + * to the newly generated counts cubes is written to path given in the + * parameter outcube. + * * Results are stored in an observation container that can be written to disk * in form of a single FITS file. On output, the observation container will * have merged the input event lists into a single observation. @@ -74,51 +81,55 @@ public: // Methods void clear(void); void run(void); + int cubes(void) const; + const GCTAEventCube& cube(const int& index = 0) const; void save(void); void publish(const std::string& name = ""); - const GCTAEventCube& cube(void) const; protected: // Protected methods - void init_members(void); - void copy_members(const ctbin& app); - void free_members(void); - void get_parameters(void); - void fill_cube(GCTAObservation* obs); - void set_weights(GCTAObservation* obs); - void obs_cube(void); + void init_members(void); + void copy_members(const ctbin& app); + void free_members(void); + void get_parameters(void); + void init_sky_dir_cache(void); + GSkyMap create_cube(const GCTAObservation* obs); + GSkyMap fill_cube(const GCTAObservation* obs); + GSkyMap set_weights(const GCTAObservation* obs); + void obs_cube_stacked(void); // User parameters - GFilename m_outcube; //!< Output counts map file name - bool m_usepnt; //!< Use pointing instead of xref/yref parameters - bool m_publish; //!< Publish counts cube? - GChatter m_chatter; //!< Chattiness + bool m_usepnt; //!< Use pointing instead of xref/yref parameters + bool m_stack; //!< Output one stacked cube or multiple cubes + std::string m_prefix; //!< Prefix for multiple cubes + bool m_publish; //!< Publish counts cube? + GChatter m_chatter; //!< Chattiness // Protected members - GSkyMap m_counts; //!< Event cube counts - GSkyMap m_weights; //!< Event cube weights - GEbounds m_ebounds; //!< Energy boundaries - GGti m_gti; //!< Good time intervals - GCTAEventCube m_cube; //!< Events cube (for cube() method) - double m_ontime; //!< Total ontime - double m_livetime; //!< Total livetime + std::vector m_cubes; //!< Event cubes + std::vector m_counts; //!< List of event cube counts + std::vector m_weights; //!< List of event cube weights + GEbounds m_ebounds; //!< Energy boundaries + GGti m_gti; //!< Stacked Good time intervals + double m_ontime; //!< Total ontime + double m_livetime; //!< Total livetime // Cache members - std::vector m_dirs; //!< Cached GSkyDir for each pixel in m_counts + std::vector m_dirs; //!< Cached GSkyDir for all pixels }; /***********************************************************************//** - * @brief Return event cube + * @brief Return number of event cubes * - * @return Reference to event cube + * @return Number of event cubes. * - * Returns a reference to the event cube. + * Returns number of event cubes. ***************************************************************************/ inline -const GCTAEventCube& ctbin::cube(void) const +int ctbin::cubes(void) const { - return m_cube; + return m_cubes.size(); } #endif /* CTBIN_HPP */ diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 2a6c521..b46b544 100755 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -16,8 +16,10 @@ # # General parameters #=================== -inobs, f, a, events.fits,,, "Input event list or observation definition XML file" -outcube, f, a, cntcube.fits,,, "Output counts cube file" +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" +outobs, f, a, cntcube.fits,,, "Output counts cube file or observation definition XML file" +stack, b, h, yes,,, "Stack events into one cube?" +prefix, s, h, cntcube_,,, "Prefix for counts cubes in observation definition XML file" # # Energy binning parameters diff --git a/test/data/obs_unbinned_two.xml b/test/data/obs_unbinned_two.xml new file mode 100644 index 0000000..8701042 --- /dev/null +++ b/test/data/obs_unbinned_two.xml @@ -0,0 +1,17 @@ + + + + + + + + + + + + + + + + + diff --git a/test/test_ctbin.py b/test/test_ctbin.py index 125ef94..e7e64bf 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -43,6 +43,9 @@ class Test(test): # Call base class constructor test.__init__(self) + # Set members + self._inobs_two = self._datadir + '/obs_unbinned_two.xml' + # Return return @@ -71,7 +74,7 @@ class Test(test): # Setup ctbin command cmd = ctbin+' inobs="'+self._events+'"'+ \ - ' outcube="cntmap_cmd1.fits"'+\ + ' outobs="cntmap_cmd1.fits"'+\ ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \ ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \ ' xref=83.63 yref=22.01 proj="CAR"'+ \ @@ -87,7 +90,7 @@ class Test(test): # Setup ctbin command cmd = ctbin+' inobs="events_that_do_not_exist.fits"'+ \ - ' outcube="cntmap_cmd2.fits"'+\ + ' outobs="cntmap_cmd2.fits"'+\ ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \ ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \ ' xref=83.63 yref=22.01 proj="CAR"'+ \ @@ -99,7 +102,7 @@ class Test(test): # Setup ctbin command with different ebounds than input cmd = ctbin+' inobs="'+self._events+'"'+ \ - ' outcube="cntmap_cmd3.fits"'+\ + ' outobs="cntmap_cmd3.fits"'+\ ' emin=2.0 emax=20.0 enumbins=10 ebinalg="LOG"'+ \ ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \ ' xref=83.63 yref=22.01 proj="CAR"'+ \ @@ -113,6 +116,40 @@ class Test(test): evt = gammalib.GCTAEventCube('cntmap_cmd3.fits') self._check_cube(evt, 115) + # Check joint (unstacked) binning + + # Setup ctbin command + cmd = ctbin+' inobs="'+self._inobs_two+'"'+ \ + ' outobs="obs_unbinned_two_binned.xml"'+\ + ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \ + ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \ + ' xref=83.63 yref=22.01 proj="CAR"'+ \ + ' logfile="ctbin_cmd4.log" chatter=1 stack=no'+ \ + ' prefix="cntmap_cmd4_"' + + # Check if execution was successful + self.test_value(self._execute(cmd), 0, + 'Check successful execution from command line with stack=no') + + # Load counts cubes and check content. + evt = gammalib.GCTAEventCube('cntmap_cmd4_00001.fits') + self._check_cube(evt, 245) + + # Load counts cubes and check content. + evt = gammalib.GCTAEventCube('cntmap_cmd4_00002.fits') + self._check_cube(evt, 245) + + # Check observation definition XML output file + obs = gammalib.GObservations('obs_unbinned_two_binned.xml') + self.test_value(obs.size(), 2, 'Check for 2 observations in XML file') + + # Check for content of observations file + self.test_value(obs[0].eventfile().file(), 'cntmap_cmd4_cta_00001.fits', + 'Check first counts cube file name') + self.test_value(obs[1].eventfile().file(), 'cntmap_cmd4_cta_00002.fits', + 'Check second counts cube file name') + + # Check ctbin --help self._check_help(ctbin) @@ -131,12 +168,12 @@ class Test(test): # cube self.test_value(binning.obs().size(), 0, 'Check that empty ctbin has an empty observation container') - self.test_value(binning.cube().size(), 0, + self.test_value(binning.cubes(), 0, 'Check that empty ctbin has an empty counts cube') # Check that saving does not nothing binning['logfile'] = 'ctbin_py0.log' - binning['outcube'] = 'ctbin_py0.fits' + binning['outobs'] = 'ctbin_py0.fits' binning.logFileOpen() binning.save() self.test_assert(not os.path.isfile('ctbin_py0.fits'), @@ -150,7 +187,7 @@ class Test(test): # Now set ctbin parameters binning['inobs'] = self._events - binning['outcube'] = 'ctbin_py1.fits' + binning['outobs'] = 'ctbin_py1.fits' binning['ebinalg'] = 'LOG' binning['emin'] = 1.0 binning['emax'] = 100.0 @@ -186,11 +223,13 @@ class Test(test): cpy_bin.logFileOpen() cpy_bin.run() - # Check content of observation container and counts cube. Now an empty - # event cube is expected since on input the observation is binned, and - # any binned observation will be skipped. - self._check_observation(cpy_bin, 0) - self._check_cube(cpy_bin.cube(), 0) + # Check content of observation container and number of counts cubes. + # There should be a single binned observation in the observation + # container, which is the one that was produced in the run before. + # Since ctbin finds no unbinned observation in the container, the + # number of cubes should be zero. + self._check_observation(cpy_bin, 245) + self.test_value(cpy_bin.cubes(), 0, 'Check that there are no counts cubes') # Save counts cube binning.save() @@ -206,7 +245,7 @@ class Test(test): # counts cube self.test_value(binning.obs().size(), 0, 'Check that empty ctbin has an empty observation container') - self.test_value(binning.cube().size(), 0, + self.test_value(binning.cubes(), 0, 'Check that empty ctbin has an empty counts cube') # Prepare observation container for stacking of events into a @@ -215,7 +254,7 @@ class Test(test): # Set-up ctbin using an observation container binning = ctools.ctbin(obs) - binning['outcube'] = 'ctbin_py3.fits' + binning['outobs'] = 'ctbin_py3.fits' binning['ebinalg'] = 'LOG' binning['emin'] = 1.0 binning['emax'] = 100.0 @@ -249,7 +288,7 @@ class Test(test): pars = {'inobs': self._events, 'ebinalg': 'LIN', 'emin': 1.0, 'emax': 100.0, 'enumbins': 10, 'nxpix': 40, 'nypix': 40, 'binsz': 0.1, 'coordsys': 'CEL', 'proj': 'CAR', - 'xref': 83.63, 'yref': 22.01, 'outcube': 'ctbin_py4.fits', + 'xref': 83.63, 'yref': 22.01, 'outobs': 'ctbin_py4.fits', 'logfile': 'ctbin_py4.log', 'chatter': 2} binning = ctools.ctbin() binning.pardict(pars) @@ -260,6 +299,42 @@ class Test(test): evt = gammalib.GCTAEventCube('ctbin_py4.fits') self._check_cube(evt, 245) + # Test unstacked version + pars = {'inobs': self._datadir+'/obs_unbinned_two.xml', 'ebinalg': 'LIN', 'emin': 1.0, + 'emax': 100.0, 'enumbins': 10, 'nxpix': 40, 'nypix': 40, + 'binsz': 0.1, 'coordsys': 'CEL', 'proj': 'CAR', + 'xref': 83.63, 'yref': 22.01, 'outobs': 'ctbin_py5.xml', + 'logfile': 'ctbin_py5.log', 'chatter': 2, + 'prefix': 'cntcube_py5_', 'stack' : False} + binning = ctools.ctbin() + binning.pardict(pars) + + binning.logFileOpen() + binning.run() + + # Check individual cubes + self._check_cube(binning.cube(0), 245) + self._check_cube(binning.cube(1), 245) + + # Save counts cube + binning.save() + + # Load observations + obs = gammalib.GObservations('ctbin_py5.xml') + self.test_value(obs.size(), 2, 'Check number of output observations') + + # Check counts cubes + evt = obs[0].events() + self._check_cube(evt, 245) + evt = obs[1].events() + self._check_cube(evt, 245) + + # Load counts cubes and check content + evt = gammalib.GCTAEventCube('cntcube_py5_00001.fits') + self._check_cube(evt, 245) + evt = gammalib.GCTAEventCube('cntcube_py5_00002.fits') + self._check_cube(evt, 245) + # Return return