diff --git a/ChangeLog b/ChangeLog index cf7e8b0..ed975f6 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-06-20 Juergen Knoedlseder +2016-06-21 Juergen Knoedlseder * Version 1.1.0 released ========================= + CTA cube response no longer selects on RoI and energies (#1789) Add weights to CTA event cube and event bin (#1789) Remove GCTACubeExposure, GCTACubePsf and GCTACubeBackground elogmeans() methods GCTACubeExposure, GCTACubePsf and GCTACubeBackground use energy values (#1788) diff --git a/NEWS b/NEWS index 2b6e1fe..54ec728 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0 Author(s) : Juergen Knoedlseder -20 June 2016 +21 June 2016 1. Introduction @@ -408,6 +408,10 @@ multiplied to the bin size in the GCTAEventBin::size() method so that model valu weigthed according to the bin weight. This enables the usage of only partially covered bins, arising for example from different energy thresholds or ROIs (#1789). +The cube response classes GCTACubeExposure, GCTACubePsf, GCTACubeEdisp and +GCTACubeBackground do no longer select on RoI and energy boundaries as this is now +handled by the weight array (#1789). + 21. COMPTEL interface --------------------- diff --git a/README.md b/README.md index 74ee090..20aafc5 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ GammaLib information ==================== -* Version: 1.1.0.dev1 (20 June 2016) +* Version: 1.1.0.dev1 (21 June 2016) * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/) diff --git a/inst/cta/src/GCTACubeBackground.cpp b/inst/cta/src/GCTACubeBackground.cpp index 6741e24..67050b5 100644 --- a/inst/cta/src/GCTACubeBackground.cpp +++ b/inst/cta/src/GCTACubeBackground.cpp @@ -339,20 +339,6 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) *log << " in background cube computation." << std::endl; } - // Extract region of interest from CTA observation - GCTARoi roi = cta->roi(); - - // Extract energy boundaries from CTA observation - GEbounds obs_ebounds = cta->ebounds(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " - "\""+cta->name()+"\". Run ctselect to specify " - "an RoI for this observation"; - throw GException::invalid_value(G_FILL, msg); - } - // Set GTI of actual observations as the GTI of the event cube eventcube.gti(cta->gti()); @@ -365,31 +351,20 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) // Get event bin GCTAEventBin* bin = eventcube[i]; - // Skip if energy is not contained within the energy boundaries of - // the observation - if (!obs_ebounds.contains(bin->energy())) { - continue; - } - - // Continue only if binned in contained in ROI - if (roi.contains(*bin)) { - - // Compute model value for event bin. The model value is - // given in counts/MeV/s/sr. - double model = obs.models().eval(*bin, *cta); - - // Multiply by livetime to get the correct weighting for - // each observation. We divide by the total livetime later - // to get the background model in units of counts/MeV/s/sr. - model *= livetime; + // Compute model value for event bin. The model value is + // given in counts/MeV/s/sr. + double model = obs.models().eval(*bin, *cta); - // Add existing number of counts - model += bin->counts(); + // Multiply by livetime to get the correct weighting for + // each observation. We divide by the total livetime later + // to get the background model in units of counts/MeV/s/sr. + model *= livetime; - // Store cumulated value (units: counts/MeV/sr) - bin->counts(model); + // Add existing number of counts + model += bin->counts(); - } // endif: bin was contained in RoI + // Store cumulated value (units: counts/MeV/sr) + bin->counts(model); } // endfor: looped over all bins diff --git a/inst/cta/src/GCTACubeEdisp.cpp b/inst/cta/src/GCTACubeEdisp.cpp index dcb544d..d9af607 100644 --- a/inst/cta/src/GCTACubeEdisp.cpp +++ b/inst/cta/src/GCTACubeEdisp.cpp @@ -428,16 +428,6 @@ void GCTACubeEdisp::set(const GCTAObservation& obs) "valid event list."; throw GException::invalid_value(G_SET, msg); } - const GCTARoi& roi = list->roi(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No valid Region of Interest found in CTA " - "observation \""+obs.name()+"\". Run ctselect " - "to set the Region of Interest for this " - "observation."; - throw GException::invalid_value(G_SET, msg); - } // Get CTA response const GCTAResponseIrf* rsp = dynamic_cast(obs.response()); @@ -469,46 +459,41 @@ void GCTACubeEdisp::set(const GCTAObservation& obs) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { + // Compute theta angle with respect to pointing direction + // in radians + double theta = pnt.dist(dir); - // Compute theta angle with respect to pointing direction - // in radians - double theta = pnt.dist(dir); + // Loop over all energy dispersion cube true energies + for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin){ - // Loop over all energy dispersion cube true energies - for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin){ + // Get log10 of true energy in TeV + double logEsrc = m_elogmeans[iebin]; - // Get log10 of true energy in TeV - double logEsrc = m_elogmeans[iebin]; - - // Loop over migration bins - for (int imigra = 0; imigra < m_migras.size(); ++imigra) { - - // Get migration for this migration bin - double migra = m_migras[imigra]; + // Loop over migration bins + for (int imigra = 0; imigra < m_migras.size(); ++imigra) { - // Skip migration bin if it's not positive - if (migra <= 0.0) { - continue; - } + // Get migration for this migration bin + double migra = m_migras[imigra]; - // Compute log10 of reconstructed energy in TeV - double logEobs = std::log10(migra) - logEsrc; + // Skip migration bin if it's not positive + if (migra <= 0.0) { + continue; + } - // Set map index - int imap = offset(imigra, iebin); + // Compute log10 of reconstructed energy in TeV + double logEobs = std::log10(migra) - logEsrc; - // Set energy dispersion cube value - m_cube(pixel, imap) = (*edisp)(logEobs, logEsrc, - theta, 0.0, - 0.0, 0.0); + // Set map index + int imap = offset(imigra, iebin); - } // endfor: looped over migration bins + // Set energy dispersion cube value + m_cube(pixel, imap) = (*edisp)(logEobs, logEsrc, + theta, 0.0, + 0.0, 0.0); - } // endfor: looped over energy bins + } // endfor: looped over migration bins - } // endif: pixel was within RoI + } // endfor: looped over energy bins } // endfor: looped over all pixels @@ -576,23 +561,6 @@ void GCTACubeEdisp::fill(const GObservations& obs, GLog* log) continue; } - // Extract region of interest from CTA observation - GCTARoi roi = cta->roi(); - - // Extract energy boundaries from CTA observation - GEbounds obs_ebounds = cta->ebounds(); - GEnergy obs_emin = obs_ebounds.emin() - g_energy_margin; - GEnergy obs_emax = obs_ebounds.emax() + g_energy_margin; - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No valid Region of Interest found in CTA " - "observation \""+cta->name()+"\". Run ctselect " - "to set the Region of Interest for this " - "observation."; - throw GException::invalid_value(G_FILL, msg); - } - // Get CTA response const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); if (rsp == NULL) { @@ -632,60 +600,47 @@ void GCTACubeEdisp::fill(const GObservations& obs, GLog* log) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { - - // Compute theta angle with respect to pointing - // direction in radians - double theta = pnt.dist(dir); - - // Loop over all energy bins - for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin) { - - // Skip if energy is not within observation energy - // boundaries - if (m_energies[iebin] < obs_emin || - m_energies[iebin] > obs_emax) { - continue; - } - - // Get log10 of true energy in TeV - double logEsrc = m_elogmeans[iebin]; + // Compute theta angle with respect to pointing direction in radians + double theta = pnt.dist(dir); - // Compute exposure weight - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * - cta->livetime(); + // Loop over all energy bins + for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin) { - // Accumulate weights - exposure(pixel, iebin) += weight; + // Get log10 of true energy in TeV + double logEsrc = m_elogmeans[iebin]; - // Loop over delta values - for (int imigra = 0; imigra < m_migras.size(); ++imigra) { + // Compute exposure weight + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * + cta->livetime(); - // Compute delta in radians - double migra = m_migras[imigra]; + // Accumulate weights + exposure(pixel, iebin) += weight; - // Skip migra bin if zero - if (migra <= 0.0) { - continue; - } + // Loop over delta values + for (int imigra = 0; imigra < m_migras.size(); ++imigra) { - // Compute log10 of reconstructed energy in TeV - double logEobs = std::log10(migra) - logEsrc; + // Compute delta in radians + double migra = m_migras[imigra]; - // Set map index - int imap = offset(imigra, iebin); + // Skip migra bin if zero + if (migra <= 0.0) { + continue; + } - // Add energy dispersion cube value - m_cube(pixel, imap) += (*edisp)(logEobs, logEsrc, - theta, 0.0, - 0.0, 0.0) * weight; + // Compute log10 of reconstructed energy in TeV + double logEobs = std::log10(migra) - logEsrc; - } // endfor: looped over migration bins + // Set map index + int imap = offset(imigra, iebin); - } // endfor: looped over energy bins + // Add energy dispersion cube value + m_cube(pixel, imap) += (*edisp)(logEobs, logEsrc, + theta, 0.0, + 0.0, 0.0) * weight; - } // endif: pixel was within RoI + } // endfor: looped over migration bins + + } // endfor: looped over energy bins } // endfor: looped over all pixels diff --git a/inst/cta/src/GCTACubeExposure.cpp b/inst/cta/src/GCTACubeExposure.cpp index b69a434..ea2dbb6 100644 --- a/inst/cta/src/GCTACubeExposure.cpp +++ b/inst/cta/src/GCTACubeExposure.cpp @@ -300,18 +300,17 @@ GCTACubeExposure* GCTACubeExposure::clone(void) const /***********************************************************************//** - * @brief Set exposure cube from one CTA observation + * @brief Set exposure cube for one CTA observation * * @param[in] obs CTA observation. * - * Set the exposure cube from a single CTA observations. The cube pixel - * values are computed as product of the effective area and the livetime. + * Set the exposure cube for a single CTA observation. The cube pixel values + * are computed as product of the effective area and the livetime. * * @todo: Throw an exception if response is not valid ***************************************************************************/ void GCTACubeExposure::set(const GCTAObservation& obs) { - // Only continue if we have an unbinned observation if (obs.eventtype() == "EventList") { @@ -320,17 +319,6 @@ void GCTACubeExposure::set(const GCTAObservation& obs) m_livetime = 0.0; m_cube = 0.0; - // Extract region of interest from CTA observation - GCTARoi roi = obs.roi(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " - "\""+obs.name()+"\". Run ctselect to specify " - "an RoI for this observation"; - throw GException::invalid_value(G_SET, msg); - } - // Get references on CTA response and pointing direction const GCTAResponseIrf* rsp = dynamic_cast(obs.response()); const GSkyDir& pnt = obs.pointing().dir(); @@ -344,26 +332,21 @@ void GCTACubeExposure::set(const GCTAObservation& obs) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { - - // Compute theta angle with respect to pointing direction - // in radians - double theta = pnt.dist(dir); - - // Loop over all exposure cube energies - for (int iebin = 0; iebin < m_energies.size(); ++iebin){ + // Compute theta angle with respect to pointing direction + // in radians + double theta = pnt.dist(dir); - // Get logE/TeV - double logE = m_energies[iebin].log10TeV(); + // Loop over all exposure cube energies + for (int iebin = 0; iebin < m_energies.size(); ++iebin){ - // Set exposure cube (effective area * livetime) - m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * - obs.livetime(); + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - } // endfor: looped over energy bins + // Set exposure cube (effective area * livetime) + m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + obs.livetime(); - } // endif: pixel was within RoI + } // endfor: looped over energy bins } // endfor: looped over all pixels @@ -430,20 +413,6 @@ void GCTACubeExposure::fill(const GObservations& obs, GLog* log) continue; } - // Extract region of interest from CTA observation - GCTARoi roi = cta->roi(); - - // Extract energy boundaries from CTA observation - GEbounds obs_ebounds = cta->ebounds(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " - "\""+cta->name()+"\". Run ctselect to specify " - "an RoI for this observation"; - throw GException::invalid_value(G_FILL, msg); - } - // Get references on CTA response and pointing direction const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); const GSkyDir& pnt = cta->pointing().dir(); @@ -476,36 +445,20 @@ void GCTACubeExposure::fill(const GObservations& obs, GLog* log) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { - - // Compute theta angle with respect to pointing - // direction in radians - double theta = pnt.dist(dir); + // Compute theta angle with respect to pointing direction in radians + double theta = pnt.dist(dir); - // Loop over all exposure cube energy bins - for (int iebin = 0; iebin < m_energies.size(); ++iebin){ - - // Only add exposure if energy is inside observation energy - // boundaries - // TODO: The exposure cube energies are true energies, while - // the observation energy boundaries are reconstructed - // energies. In principle one would need to use the - // energy dispersion information to do this correctly. - if (obs_ebounds.contains(m_energies[iebin])) { - - // Get logE/TeV - double logE = m_energies[iebin].log10TeV(); + // Loop over all exposure cube energy bins + for (int iebin = 0; iebin < m_energies.size(); ++iebin){ - // Add to exposure cube (effective area * livetime) - m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * - cta->livetime(); + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - } // endif: energy bin was inside observation energy boundaries - - } // endfor: looped over energy bins + // Add to exposure cube (effective area * livetime) + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + cta->livetime(); - } // endif: pixel within RoI + } // endfor: looped over energy bins } // endfor: looped over all pixels diff --git a/inst/cta/src/GCTACubePsf.cpp b/inst/cta/src/GCTACubePsf.cpp index d7a5c8d..80e5aa9 100644 --- a/inst/cta/src/GCTACubePsf.cpp +++ b/inst/cta/src/GCTACubePsf.cpp @@ -51,7 +51,6 @@ /* __ Debug definitions __________________________________________________ */ /* __ Constants __________________________________________________________ */ -const GEnergy g_energy_margin(1.0e-12, "TeV"); /*========================================================================== @@ -378,25 +377,6 @@ void GCTACubePsf::set(const GCTAObservation& obs) // Only continue if we have an unbinned observation if (obs.eventtype() == "EventList") { - // Extract region of interest from CTA observation - const GCTAEventList* list = dynamic_cast(obs.events()); - if (list == NULL) { - std::string msg = "CTA Observation does not contain an event " - "list. Event list information is needed to " - "retrieve the Region of Interest for each " - "CTA observation."; - throw GException::invalid_value(G_SET, msg); - } - const GCTARoi& roi = list->roi(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " - "\""+obs.name()+"\". Run ctselect to specify " - "an RoI for this observation"; - throw GException::invalid_value(G_SET, msg); - } - // Get references on CTA response and pointing direction const GCTAResponseIrf* rsp = dynamic_cast(obs.response()); const GSkyDir& pnt = obs.pointing().dir(); @@ -410,36 +390,31 @@ void GCTACubePsf::set(const GCTAObservation& obs) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { - - // Compute theta angle with respect to pointing direction - // in radians - double theta = pnt.dist(dir); - - // Loop over all exposure cube energy bins - for (int iebin = 0; iebin < m_energies.size(); ++iebin){ + // Compute theta angle with respect to pointing direction + // in radians + double theta = pnt.dist(dir); - // Get logE/TeV - double logE = m_energies[iebin].log10TeV(); + // Loop over all exposure cube energy bins + for (int iebin = 0; iebin < m_energies.size(); ++iebin){ - // Loop over delta values - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { - - // Compute delta in radians - double delta = m_deltas[idelta] * gammalib::deg2rad; + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - // Set map index - int imap = offset(idelta, iebin); + // Loop over delta values + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { - // Set PSF cube - m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE); + // Compute delta in radians + double delta = m_deltas[idelta] * gammalib::deg2rad; - } // endfor: looped over delta bins + // Set map index + int imap = offset(idelta, iebin); - } // endfor: looped over energy bins + // Set PSF cube + m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE); - } // endif: pixel was within RoI + } // endfor: looped over delta bins + + } // endfor: looped over energy bins } // endfor: looped over all pixels @@ -501,20 +476,6 @@ void GCTACubePsf::fill(const GObservations& obs, GLog* log) continue; } - // Extract region of interest from CTA observation - GCTARoi roi = cta->roi(); - - // Extract energy boundaries from CTA observation - GEbounds obs_ebounds = cta->ebounds(); - - // Check for RoI sanity - if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " - "\""+cta->name()+"\". Run ctselect to specify " - "an RoI for this observation"; - throw GException::invalid_value(G_FILL, msg); - } - // Get references on CTA response and pointing direction const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); const GSkyDir& pnt = cta->pointing().dir(); @@ -547,54 +508,38 @@ void GCTACubePsf::fill(const GObservations& obs, GLog* log) // Get pixel sky direction GSkyDir dir = m_cube.inx2dir(pixel); - // Continue only if pixel is within RoI - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) { - - // Compute theta angle with respect to pointing - // direction in radians - double theta = pnt.dist(dir); - - // Loop over all energy bins - for (int iebin = 0; iebin < m_energies.size(); ++iebin) { - - // Only add exposure if energy is inside observation energy - // boundaries - // TODO: The exposure cube energies are true energies, while - // the observation energy boundaries are reconstructed - // energies. In principle one would need to use the - // energy dispersion information to do this correctly. - if (!obs_ebounds.contains(m_energies[iebin])) { - continue; - } + // Compute theta angle with respect to pointing direction in radians + double theta = pnt.dist(dir); - // Get logE/TeV - double logE = m_energies[iebin].log10TeV(); + // Loop over all energy bins + for (int iebin = 0; iebin < m_energies.size(); ++iebin) { - // Compute exposure weight - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * - cta->livetime(); + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - // Accumulate weights - exposure(pixel, iebin) += weight; + // Compute exposure weight + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + cta->livetime(); - // Loop over delta values - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { + // Accumulate weights + exposure(pixel, iebin) += weight; - // Compute delta in radians - double delta = m_deltas[idelta] * gammalib::deg2rad; + // Loop over delta values + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { - // Set map index - int imap = offset(idelta, iebin); + // Compute delta in radians + double delta = m_deltas[idelta] * gammalib::deg2rad; - // Add on PSF cube - m_cube(pixel, imap) += - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; + // Set map index + int imap = offset(idelta, iebin); - } // endfor: looped over delta bins + // Add on PSF cube + m_cube(pixel, imap) += + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; - } // endfor: looped over energy bins + } // endfor: looped over delta bins - } // endif: pixel was within RoI + } // endfor: looped over energy bins } // endfor: looped over all pixels