diff --git a/ChangeLog b/ChangeLog index 74265b9..3fdc03d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.2.0 released ========================= + Skip zero livetime observations in stacked cubes for CTA (#1920) Correct handling of time references in GGti::extend() method (#1922) Check that MC energy is comprised in validity range for CTA 3D background template (#1918) Add GModelTemporalPhaseCurve class (#1916) diff --git a/NEWS b/NEWS index e5b8770..7bb1123 100644 --- a/NEWS +++ b/NEWS @@ -297,6 +297,9 @@ The GCTABackground3D::mc() method now verifies that the specified energy is with the validity range of the background template. Otherwise an exception is thrown (#1918). +Observations with zero livetime are now skipped in building the cubes for stacked +analysis for CTA (#1920). + 21. COMPTEL interface --------------------- diff --git a/inst/cta/src/GCTACubeBackground.cpp b/inst/cta/src/GCTACubeBackground.cpp index 3556439..9c8a204 100644 --- a/inst/cta/src/GCTACubeBackground.cpp +++ b/inst/cta/src/GCTACubeBackground.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTACubeBackground.cpp - CTA cube background class * * ----------------------------------------------------------------------- * - * copyright (C) 2015-2016 by Michael Mayer * + * copyright (C) 2015-2017 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -373,6 +373,22 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) continue; } + // Get observation livetime + double livetime = cta->livetime(); + + // Skip observation if livetime is zero + if (livetime == 0.0) { + if (log != NULL) { + *log << "Skipping unbinned "; + *log << cta->instrument(); + *log << " observation "; + *log << "\"" << cta->name() << "\""; + *log << " (id=" << cta->id() << ") due to zero livetime"; + *log << std::endl; + } + continue; + } + // Announce observation usage if (log != NULL) { *log << "Including "; @@ -399,9 +415,6 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) // Set GTI of actual observations as the GTI of the event cube eventcube.gti(cta->gti()); - // Get observation livetime - double livetime = cta->livetime(); - // Loop over all bins in background cube for (int i = 0; i < eventcube.size(); ++i) { diff --git a/inst/cta/src/GCTACubeEdisp.cpp b/inst/cta/src/GCTACubeEdisp.cpp index 1900a00..8ddba8d 100644 --- a/inst/cta/src/GCTACubeEdisp.cpp +++ b/inst/cta/src/GCTACubeEdisp.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTACubeEdisp.cpp - CTA cube analysis energy dispersion class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Michael Mayer * + * copyright (C) 2016-2017 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -880,82 +880,101 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation& obs, throw GException::invalid_value(G_FILL_CUBE, msg); } - // Announce observation usage - if (log != NULL) { - *log << "Including "; - *log << obs.instrument(); - *log << " observation \"" << obs.name(); - *log << "\" (id=" << obs.id() << ")"; - *log << " in energy dispersion cube computation." << std::endl; + // Skip observation if livetime is zero + if (obs.livetime() == 0.0) { + if (log != NULL) { + *log << "Skipping unbinned "; + *log << obs.instrument(); + *log << " observation "; + *log << "\"" << obs.name() << "\""; + *log << " (id=" << obs.id() << ") due to zero livetime"; + *log << std::endl; + } } - // Loop over all pixels in sky map - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { - - // Get pixel sky direction - GSkyDir dir = m_cube.inx2dir(pixel); + // ... otherwise continue + else { - // Skip pixel if it is outside the RoI - if (roi.centre().dir().dist_deg(dir) > roi.radius()) { - continue; + // Announce observation usage + if (log != NULL) { + *log << "Including "; + *log << obs.instrument(); + *log << " observation \"" << obs.name(); + *log << "\" (id=" << obs.id() << ")"; + *log << " in energy dispersion cube computation." << std::endl; } - - // 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) { - - // Skip energy dispersion cube energy if the energy is - // outside the observation energy. The energy dispersion cube - // energies are true energies while the observation energy - // boundaries are reconstructed energies, hence this is only - // an approximation, but probably the only we can really do. - if (!obs_ebounds.contains(m_energies[iebin])) { - continue; - } - // Get log10 of true energy in TeV - double logEsrc = m_energies[iebin].log10TeV(); + // Loop over all pixels in sky map + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { - // Compute exposure weight - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * - obs.livetime(); + // Get pixel sky direction + GSkyDir dir = m_cube.inx2dir(pixel); - // If available, accumulate weights - if (exposure != NULL) { - (*exposure)(pixel, iebin) += weight; + // Skip pixel if it is outside the RoI + if (roi.centre().dir().dist_deg(dir) > roi.radius()) { + continue; } + + // 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) { + + // Skip energy dispersion cube energy if the energy is + // outside the observation energy. The energy dispersion + // cube energies are true energies while the observation + // energy boundaries are reconstructed energies, hence + // this is only an approximation, but probably the only + // we can really do. + if (!obs_ebounds.contains(m_energies[iebin])) { + continue; + } - // Loop over delta values - for (int imigra = 0; imigra < m_migras.size(); ++imigra) { + // Get log10 of true energy in TeV + double logEsrc = m_energies[iebin].log10TeV(); - // Compute energy migration fraction - double migra = m_migras[imigra]; + // Compute exposure weight + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * + obs.livetime(); - // Skip migra bin if zero - if (migra <= 0.0) { - continue; + // If available, accumulate weights + if (exposure != NULL) { + (*exposure)(pixel, iebin) += weight; } - // Compute log10 of reconstructed energy in TeV - double logEobs = std::log10(migra) + logEsrc; - GEnergy Eobs(std::pow(10.0, logEobs), "TeV"); + // Loop over delta values + for (int imigra = 0; imigra < m_migras.size(); ++imigra) { - // Set map index - int imap = offset(imigra, iebin); + // Compute energy migration fraction + double migra = m_migras[imigra]; + + // Skip migra bin if zero + if (migra <= 0.0) { + continue; + } + + // Compute log10 of reconstructed energy in TeV + double logEobs = std::log10(migra) + logEsrc; + GEnergy Eobs(std::pow(10.0, logEobs), "TeV"); - // Add energy dispersion cube value - m_cube(pixel, imap) += rsp->edisp(Eobs, - theta, 0.0, - 0.0, 0.0, - logEsrc) * weight; + // Set map index + int imap = offset(imigra, iebin); - } // endfor: looped over migration bins + // Add energy dispersion cube value + m_cube(pixel, imap) += rsp->edisp(Eobs, + theta, 0.0, + 0.0, 0.0, + logEsrc) * weight; + + } // endfor: looped over migration bins - } // endfor: looped over energy bins + } // endfor: looped over energy bins - } // endfor: looped over all pixels + } // endfor: looped over all pixels + + } // endelse: livetime was not zero } // endif: observation contained an event list diff --git a/inst/cta/src/GCTACubeExposure.cpp b/inst/cta/src/GCTACubeExposure.cpp index c8d35f5..2a34b8f 100644 --- a/inst/cta/src/GCTACubeExposure.cpp +++ b/inst/cta/src/GCTACubeExposure.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTACubeExposure.cpp - CTA cube analysis exposure class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Chia-Chun Lu * + * copyright (C) 2014-2017 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -656,56 +656,74 @@ void GCTACubeExposure::fill_cube(const GCTAObservation& obs, GLog* log) throw GException::invalid_value(G_FILL_CUBE, msg); } - // Announce observation usage - if (log != NULL) { - *log << "Including "; - *log << obs.instrument(); - *log << " observation \"" << obs.name(); - *log << "\" (id=" << obs.id() << ")"; - *log << " in exposure cube computation." << std::endl; + // Skip observation if livetime is zero + if (obs.livetime() == 0.0) { + if (log != NULL) { + *log << "Skipping unbinned "; + *log << obs.instrument(); + *log << " observation "; + *log << "\"" << obs.name() << "\""; + *log << " (id=" << obs.id() << ") due to zero livetime"; + *log << std::endl; + } } - // Loop over all pixels in sky map - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { - - // Get pixel sky direction - GSkyDir dir = m_cube.inx2dir(pixel); + // ... otherwise continue + else { - // Skip pixel if it is outside the RoI - if (roi.centre().dir().dist_deg(dir) > roi.radius()) { - continue; + // Announce observation usage + if (log != NULL) { + *log << "Including "; + *log << obs.instrument(); + *log << " observation \"" << obs.name(); + *log << "\" (id=" << obs.id() << ")"; + *log << " in exposure cube computation." << std::endl; } - - // 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){ - - // Skip exposure cube energy if the energy is outside the - // observation energy. The exposure cube energies are true - // energies while the observation energy boundaries are - // reconstructed energies, hence this is only an approximation, - // but probably the only we can really do. - if (!obs_ebounds.contains(m_energies[iebin])) { + + // Loop over all pixels in sky map + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { + + // Get pixel sky direction + GSkyDir dir = m_cube.inx2dir(pixel); + + // Skip pixel if it is outside the RoI + if (roi.centre().dir().dist_deg(dir) > roi.radius()) { 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 exposure cube energies + 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) * - obs.livetime(); + // Skip exposure cube energy if the energy is outside the + // observation energy. The exposure cube energies are true + // energies while the observation energy boundaries are + // reconstructed energies, hence this is only an + // approximation, but probably the only we can really do. + if (!obs_ebounds.contains(m_energies[iebin])) { + continue; + } - } // endfor: looped over energy bins + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - } // endfor: looped over all pixels + // Add to exposure cube (effective area * livetime) + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + obs.livetime(); + + } // endfor: looped over energy bins + + } // endfor: looped over all pixels + + // Append GTIs and increment livetime + m_gti.extend(obs.gti()); + m_livetime += obs.livetime(); + + } // endelse: livetime was not zero - // Append GTIs and increment livetime - m_gti.extend(obs.gti()); - m_livetime += obs.livetime(); - } // endif: observation contained an event list // Return diff --git a/inst/cta/src/GCTACubePsf.cpp b/inst/cta/src/GCTACubePsf.cpp index 91b38bf..b37b78a 100644 --- a/inst/cta/src/GCTACubePsf.cpp +++ b/inst/cta/src/GCTACubePsf.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTACubePsf.cpp - CTA cube analysis point spread function class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Chia-Chun Lu * + * copyright (C) 2014-2017 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -787,71 +787,89 @@ void GCTACubePsf::fill_cube(const GCTAObservation& obs, throw GException::invalid_value(G_FILL_CUBE, msg); } - // Announce observation usage - if (log != NULL) { - *log << "Including "; - *log << obs.instrument(); - *log << " observation \"" << obs.name(); - *log << "\" (id=" << obs.id() << ")"; - *log << " in point spread function cube computation." << std::endl; + // Skip observation if livetime is zero + if (obs.livetime() == 0.0) { + if (log != NULL) { + *log << "Skipping unbinned "; + *log << obs.instrument(); + *log << " observation "; + *log << "\"" << obs.name() << "\""; + *log << " (id=" << obs.id() << ") due to zero livetime"; + *log << std::endl; + } } - // Loop over all pixels in sky map - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { + // ... otherwise continue + else { - // Get pixel sky direction - GSkyDir dir = m_cube.inx2dir(pixel); - - // Skip pixel if it is outside the RoI - if (roi.centre().dir().dist_deg(dir) > roi.radius()) { - continue; + // Announce observation usage + if (log != NULL) { + *log << "Including "; + *log << obs.instrument(); + *log << " observation \"" << obs.name(); + *log << "\" (id=" << obs.id() << ")"; + *log << " in point spread function cube computation." << std::endl; } - // 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) { + // Loop over all pixels in sky map + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) { - // Skip PSF cube energy if the energy is outside the observation - // energy. The PSF cube energies are true energies while the - // observation energy boundaries are reconstructed energies, - // hence this is only an approximation, but probably the only - // we can really do. - if (!obs_ebounds.contains(m_energies[iebin])) { + // Get pixel sky direction + GSkyDir dir = m_cube.inx2dir(pixel); + + // Skip pixel if it is outside the RoI + if (roi.centre().dir().dist_deg(dir) > roi.radius()) { continue; } - // Get logE/TeV - double logE = m_energies[iebin].log10TeV(); + // 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, logE) * - obs.livetime(); + // Loop over all energy bins + for (int iebin = 0; iebin < m_energies.size(); ++iebin) { - // If available, accumulate weights - if (exposure != NULL) { - (*exposure)(pixel, iebin) += weight; - } + // Skip PSF cube energy if the energy is outside the + // observation energy. The PSF cube energies are true + // energies while the observation energy boundaries are + // reconstructed energies, hence this is only an + // approximation, but probably the only we can really do. + if (!obs_ebounds.contains(m_energies[iebin])) { + continue; + } - // Loop over delta values - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { + // Get logE/TeV + double logE = m_energies[iebin].log10TeV(); - // Compute delta in radians - double delta = m_deltas[idelta] * gammalib::deg2rad; + // Compute exposure weight + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + obs.livetime(); - // Set map index - int imap = offset(idelta, iebin); + // If available, accumulate weights + if (exposure != NULL) { + (*exposure)(pixel, iebin) += weight; + } - // Add on PSF cube - m_cube(pixel, imap) += - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; + // Loop over delta values + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { - } // endfor: looped over delta bins + // Compute delta in radians + double delta = m_deltas[idelta] * gammalib::deg2rad; - } // endfor: looped over energy bins + // Set map index + int imap = offset(idelta, iebin); - } // endfor: looped over all pixels + // Add on PSF cube + m_cube(pixel, imap) += + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; + + } // endfor: looped over delta bins + + } // endfor: looped over energy bins + + } // endfor: looped over all pixels + + } // endelse: livetime was not zero } // endif: observation contained an event list