diff --git a/ChangeLog b/ChangeLog index d50d41b..201ba00 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-11-03 +2018-11-07 * Version 1.6.0 released ======================== + Enable ctbin to produce counts cubes for a joint binned analysis (#2714) Add code generator (#2709) Add csbkgmodel script (#2694) Compute DETX and DETY in ctmodel if they are not available (#2694) diff --git a/NEWS b/NEWS index 44d5d2a..9a99096 100755 --- a/NEWS +++ b/NEWS @@ -55,6 +55,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/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 77dbf6d..2c8a5fa 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -41,6 +41,8 @@ Bug fixes Improvements ------------ +* [`2714 `_] - + Enable :ref:`ctbin` to produce counts cubes for a joint binned analysis * [`2685 `_] - Add ``ctbkgcube::cntcube()`` method Account for varying energy thresholds in :ref:`ctbin` weight computation diff --git a/pyext/ctbin.i b/pyext/ctbin.i index b3cab33..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,10 +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; - const GCTAEventCube& cube(const size_t& index) const; }; diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index d41128c..20e13b1 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"); @@ -255,33 +259,33 @@ void ctbin::run(void) log_string(NORMAL, msg); } - // Initialise cubes and ids if unstacked - if (!m_stack) { - m_cubes = std::vector(obs_list.size()); - m_ids = std::vector(obs_list.size()); + // 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); + m_ids = std::vector(nobs); + } + else { + m_counts.clear(); + m_weights.clear(); + m_cubes.clear(); + m_ids.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]; - // 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_FILL_CUBE, msg); - } - // Fill the cube GSkyMap counts = fill_cube(obs); @@ -291,88 +295,118 @@ void ctbin::run(void) // Dispose events to free memory obs->dispose_events(); + // Post process cubes #pragma omp critical(ctbin_fill_cube) - if (m_stack) { - - - m_counts += counts; - - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - for (int pixel = 0; pixel < m_counts.npix(); ++pixel) { - if (weights(pixel, iebin) == 1.0) { - m_weights(pixel, iebin) = 1.0; - } - } - } - - - // If GTI is empty then set its time reference from the observation. - // From then on we keep that time reference + { + // 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, weights and observation ID + else { + m_counts[i] = counts; + m_weights[i] = weights; + m_ids[i] = obs->id(); + } + + } // end: omp critical section + + } // endfor: looped over observations + + // 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]; + + // 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(events->gti().reference()); + m_gti.reference(gti.reference()); } - + // Append GTIs - m_gti.extend(events->gti()); - + m_gti.extend(gti); + // Update ontime and livetime m_ontime += obs->ontime(); m_livetime += obs->livetime(); - - } else { // don't stack cubes - // Construct a cube from this observation - GCTAEventCube cube(counts, weights, m_ebounds, events->gti()); - m_cubes[i] = cube; - obs->events(cube); + } // 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 (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; + } + } + } + } - // Store id of the observation - m_ids[i] = obs->id(); - - } + // Set a single cube in the observation container + obs_cube_stacked(); + } // endif: stacked cube post processing - } // endfor: looped over observations + // Post process data for non-stacked cubes + else { - if (m_stack) { - // Compute total livetime in all observations - double livetime = 0.0; - for (int i = 0; i < obs_list.size(); ++i) { - livetime += obs_list[i]->livetime(); - } - - // 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(); - } - } - 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; - } - } - } - } + // Loop over all relevant observations + for (int i = 0; i < nobs; ++i) { - // 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); + // Get pointer to observation + GCTAObservation* obs = obs_list[i]; - // Set a single cube in the observation container - obs_cube(); - } + // 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]); + + // Store ID of the observation for saving of XML file + m_ids[i] = obs->id(); + + } // endfor: looped over relevant observations + + } // endelse: non-stacked cube post processing // Write resulting observation container into logger log_observations(NORMAL, m_obs, "Binned observation"); @@ -388,72 +422,108 @@ 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"); - if (!m_stack) { - - // Loop over all observations - for (size_t i = 0; i < m_obs.size(); ++i) { - - // Loop over all ids of observations that were binned - for (size_t j = 0; j < m_ids.size(); ++j) { - - // Check if ids are equal - if (m_ids[j] == m_obs[i]->id()) { - - // Construct filename - std::string filename(m_ids[j]); - filename.insert(0, m_prefix); - filename.append(".fits"); - - // Set filename for counts cube - GCTAObservation* obs = dynamic_cast(m_obs[i]); - obs->eventfile(filename); - - // Save countscube - obs->save(filename, clobber()); - } - } - } - if ((*this)["outcube"].is_valid()) { - // Get counts cube filename - m_outcube = (*this)["outcube"].filename(); - m_obs.save(m_outcube); - } - - } else { // else: stack - - // Save only if filename is valid and if there is at least one - // observation - if ((*this)["outcube"].is_valid() && m_obs.size() > 0) { - - // Get counts cube filename - m_outcube = (*this)["outcube"].filename(); - - // 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", m_outcube.url()); - - // Save cube - obs->save(m_outcube, clobber()); - - } // endif: observation was valid - - } // endif: outcube file was valid - - } // endelse: stack + // If counts cube was stacked then save first observation + if (m_stack) { + + // Save only if filename is valid and if there is at least one + // observation + if ((*this)["outcube"].is_valid() && m_obs.size() > 0) { + + // Get counts cube filename + m_outcube = (*this)["outcube"].filename(); + + // 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", m_outcube.url()); + + // Save cube + obs->save(m_outcube, clobber()); + + } // endif: observation was valid + + } // endif: outcube 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) { + + // Loop over all ids of observations that were binned + for (int j = 0; j < m_ids.size(); ++j) { + + // Skip if ID differes + if (m_ids[j] != m_obs[i]->id()) { + continue; + } + + // Construct filename + std::string filename(m_ids[j]); + filename.insert(0, m_prefix); + filename.append(".fits"); + + // Set filename for counts cube + GCTAObservation* obs = dynamic_cast(m_obs[i]); + obs->eventfile(filename); + + // Save countscube + obs->save(filename, clobber()); + + } // endfor: looped over all IDs + + } // endfor: looped over all observations + + // Write XML file + if ((*this)["outcube"].is_valid()) { + m_outcube = (*this)["outcube"].filename(); + m_obs.save(m_outcube); + } + + } // endelse: counts cube were not stacked + // Return return; } @@ -478,8 +548,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; @@ -505,15 +577,14 @@ void ctbin::init_members(void) m_chatter = static_cast(2); // Initialise protected members + m_cubes.clear(); m_counts.clear(); m_weights.clear(); + m_ids.clear(); m_ebounds.clear(); m_gti.clear(); - m_cube.clear(); m_ontime = 0.0; m_livetime = 0.0; - m_cubes.clear(); - m_ids.clear(); // Initialise cache members m_dirs.clear(); @@ -538,15 +609,14 @@ void ctbin::copy_members(const ctbin& app) m_chatter = app.m_chatter; // Copy protected members + m_cubes = app.m_cubes; m_counts = app.m_counts; m_weights = app.m_weights; + m_ids = app.m_ids; 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; - m_cubes = app.m_cubes; - m_ids = app.m_ids; // Copy cache members m_dirs = app.m_dirs; @@ -581,14 +651,8 @@ 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(); @@ -608,13 +672,6 @@ void ctbin::get_parameters(void) (*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()); - } - // Set number of OpenMP threads #ifdef _OPENMP int nthreads = (*this)["nthreads"].integer(); @@ -632,6 +689,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. @@ -641,10 +787,10 @@ void ctbin::get_parameters(void) * * Fills the events from an event list into the counts cube. ***************************************************************************/ -GSkyMap 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) { @@ -655,7 +801,7 @@ GSkyMap 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()) { @@ -667,7 +813,7 @@ GSkyMap 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; @@ -676,9 +822,12 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs) int num_outside_ebds = 0; int num_in_map = 0; - // Initialise new empty counts map - GSkyMap counts(m_counts); - counts *= 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) { @@ -699,7 +848,7 @@ GSkyMap 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++; @@ -708,8 +857,8 @@ GSkyMap 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; } @@ -726,7 +875,6 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs) } // Fill event in skymap - //#pragma omp critical(ctbin_fill_cube) counts(pixel, iebin) += 1.0; // Increment number of maps @@ -734,23 +882,6 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs) } // 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) { @@ -779,21 +910,10 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs) * Sets the counts cube weights for all bins that are considered for the * specific observation to unity. ***************************************************************************/ -GSkyMap 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()) { @@ -804,19 +924,32 @@ GSkyMap 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); - // Initialise new empty weights map - GSkyMap weights(m_weights); - weights *= 0.; + // 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 @@ -828,7 +961,6 @@ GSkyMap ctbin::set_weights(GCTAObservation* obs) } // Signal that bin was filled - //#pragma omp critical(ctbin_set_weights) weights(pixel, iebin) = 1.0; } // endfor: looped over energy layers of counts cube @@ -851,8 +983,11 @@ GSkyMap 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 @@ -865,9 +1000,9 @@ void ctbin::obs_cube(void) // Change the event type if we had an unbinned observation if (obs->eventtype() == "EventList") { - + // Assign cube to the observation - obs->events(this->cube()); + obs->events(m_cubes[0]); } @@ -877,7 +1012,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); @@ -899,49 +1035,49 @@ void ctbin::obs_cube(void) // Allocate observation container GObservations container; - // Allocate CTA observation. - GCTAObservation obs; - - // Attach event cube to CTA observation - obs.events(this->cube()); - - // Compute average pointing direction for all CTA event lists - double ra = 0.0; - double dec = 0.0; - double number = 0.0; - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* cta = dynamic_cast(m_obs[i]); - if ((cta != NULL) && (cta->eventtype() == "EventList")) { - ra += cta->pointing().dir().ra(); - dec += cta->pointing().dir().dec(); - number += 1.0; - } - } - if (number > 0.0) { - ra /= number; - dec /= number; - } - GSkyDir dir; - dir.radec(ra, dec); - GCTAPointing pointing(dir); - - // Compute deadtime correction - double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; - - // Set CTA observation attributes - obs.pointing(pointing); - obs.ra_obj(dir.ra_deg()); //!< Dummy - obs.dec_obj(dir.dec_deg()); //!< Dummy - obs.ontime(m_ontime); - obs.livetime(m_livetime); - obs.deadc(deadc); - - // Append CTA observation - container.append(obs); + // Allocate CTA observation. + GCTAObservation obs; + + // Attach event cube to CTA observation + obs.events(m_cubes[0]); + + // Compute average pointing direction for all CTA event lists + double ra = 0.0; + double dec = 0.0; + double number = 0.0; + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* cta = dynamic_cast(m_obs[i]); + if ((cta != NULL) && (cta->eventtype() == "EventList")) { + ra += cta->pointing().dir().ra(); + dec += cta->pointing().dir().dec(); + number += 1.0; + } + } + if (number > 0.0) { + ra /= number; + dec /= number; + } + GSkyDir dir; + dir.radec(ra, dec); + GCTAPointing pointing(dir); + + // Compute deadtime correction + double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; + + // Set CTA observation attributes + obs.pointing(pointing); + obs.ra_obj(dir.ra_deg()); //!< Dummy + obs.dec_obj(dir.dec_deg()); //!< Dummy + obs.ontime(m_ontime); + obs.livetime(m_livetime); + obs.deadc(deadc); // Set models in observation container container.models(m_obs.models()); + // Append CTA observation + container.append(obs); + // Copy over all remaining non-CTA observations for (int i = 0; i < m_obs.size(); ++i) { GCTAObservation* obs = dynamic_cast(m_obs[i]); diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index fa4034d..cc872a5 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 * @@ -81,10 +81,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; - const GCTAEventCube& cube(const size_t& index) const; protected: // Protected methods @@ -92,60 +92,46 @@ protected: void copy_members(const ctbin& app); void free_members(void); void get_parameters(void); - GSkyMap fill_cube(GCTAObservation* obs); - GSkyMap set_weights(GCTAObservation* obs); - void obs_cube(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_stack; //!< Output one stacked cube or multiple cubes - std::string m_prefix; //!< Prefix for output path of multiple cubes - bool m_publish; //!< Publish counts cube? - GChatter m_chatter; //!< Chattiness + GFilename m_outcube; //!< Output counts map file name + 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 output path of 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 (for cube() method) - std::vector m_ids; //!< Observation ids + std::vector m_cubes; //!< Event cubes + std::vector m_counts; //!< List of event cube counts + std::vector m_weights; //!< List of event cube weights + std::vector m_ids; //!< List of observation ids + 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; -} - - -/***********************************************************************//** - * @brief Return event cube at index - * - * @param[in] index Index of cube - * @return Reference to event cube - * - * Returns a reference to the event cube at the given index. - ***************************************************************************/ -inline -const GCTAEventCube& ctbin::cube(const size_t& index) const -{ - return m_cubes[index]; + return m_cubes.size(); } #endif /* CTBIN_HPP */ diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 6af6516..351e0a6 100755 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -18,7 +18,7 @@ #=================== inobs, f, a, events.fits,,, "Input event list or observation definition XML file" outcube, f, a, cntcube.fits,,, "Output counts cube file" -stack, b, h, yes,,, "Stack events to one cube?" +stack, b, h, yes,,, "Stack events into one cube?" prefix, s, h, binned_,,, "Prefix for counts cubes in observation definition XML file" # diff --git a/test/test_ctbin.py b/test/test_ctbin.py index e1dd911..7ad96e0 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -165,7 +165,7 @@ 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 @@ -220,11 +220,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() @@ -240,7 +242,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