Dmitry Malyshev / gammalib

Commit c616af2744cf523eff3fb9d24f90740aeb6a30dc

Authored by Jürgen Knödlseder
1 parent 8aa2978d

Skip observations with zero livetime in stacked cube building for CTA (#1920)

Inline Side-by-side
Showing 6 changed files with 223 additions and 151 deletions
ChangeLog
  View file @c616af2
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.2.0 released
4 4 =========================
5 5  
  6 + Skip zero livetime observations in stacked cubes for CTA (#1920)
6 7 Correct handling of time references in GGti::extend() method (#1922)
7 8 Check that MC energy is comprised in validity range for CTA 3D background template (#1918)
8 9 Add GModelTemporalPhaseCurve class (#1916)
... ...
... ... @@ -297,6 +297,9 @@ The GCTABackground3D::mc() method now verifies that the specified energy is with
297 297 the validity range of the background template. Otherwise an exception is thrown
298 298 (#1918).
299 299  
  300 +Observations with zero livetime are now skipped in building the cubes for stacked
  301 +analysis for CTA (#1920).
  302 +
300 303  
301 304 21. COMPTEL interface
302 305 ---------------------
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @c616af2
1 1 /***************************************************************************
2 2 * GCTACubeBackground.cpp - CTA cube background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015-2016 by Michael Mayer *
  4 + * copyright (C) 2015-2017 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -373,6 +373,22 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log)
373 373 continue;
374 374 }
375 375  
  376 + // Get observation livetime
  377 + double livetime = cta->livetime();
  378 +
  379 + // Skip observation if livetime is zero
  380 + if (livetime == 0.0) {
  381 + if (log != NULL) {
  382 + *log << "Skipping unbinned ";
  383 + *log << cta->instrument();
  384 + *log << " observation ";
  385 + *log << "\"" << cta->name() << "\"";
  386 + *log << " (id=" << cta->id() << ") due to zero livetime";
  387 + *log << std::endl;
  388 + }
  389 + continue;
  390 + }
  391 +
376 392 // Announce observation usage
377 393 if (log != NULL) {
378 394 *log << "Including ";
... ... @@ -399,9 +415,6 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
399 415 // Set GTI of actual observations as the GTI of the event cube
400 416 eventcube.gti(cta->gti());
401 417  
402   - // Get observation livetime
403   - double livetime = cta->livetime();
404   -
405 418 // Loop over all bins in background cube
406 419 for (int i = 0; i < eventcube.size(); ++i) {
407 420  
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @c616af2
1 1 /***************************************************************************
2 2 * GCTACubeEdisp.cpp - CTA cube analysis energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Michael Mayer *
  4 + * copyright (C) 2016-2017 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -880,82 +880,101 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
880 880 throw GException::invalid_value(G_FILL_CUBE, msg);
881 881 }
882 882  
883   - // Announce observation usage
884   - if (log != NULL) {
885   - *log << "Including ";
886   - *log << obs.instrument();
887   - *log << " observation \"" << obs.name();
888   - *log << "\" (id=" << obs.id() << ")";
889   - *log << " in energy dispersion cube computation." << std::endl;
  883 + // Skip observation if livetime is zero
  884 + if (obs.livetime() == 0.0) {
  885 + if (log != NULL) {
  886 + *log << "Skipping unbinned ";
  887 + *log << obs.instrument();
  888 + *log << " observation ";
  889 + *log << "\"" << obs.name() << "\"";
  890 + *log << " (id=" << obs.id() << ") due to zero livetime";
  891 + *log << std::endl;
  892 + }
890 893 }
891 894  
892   - // Loop over all pixels in sky map
893   - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
894   -
895   - // Get pixel sky direction
896   - GSkyDir dir = m_cube.inx2dir(pixel);
  895 + // ... otherwise continue
  896 + else {
897 897  
898   - // Skip pixel if it is outside the RoI
899   - if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
900   - continue;
  898 + // Announce observation usage
  899 + if (log != NULL) {
  900 + *log << "Including ";
  901 + *log << obs.instrument();
  902 + *log << " observation \"" << obs.name();
  903 + *log << "\" (id=" << obs.id() << ")";
  904 + *log << " in energy dispersion cube computation." << std::endl;
901 905 }
902   -
903   - // Compute theta angle with respect to pointing direction in radians
904   - double theta = pnt.dist(dir);
905   -
906   - // Loop over all energy bins
907   - for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
908   -
909   - // Skip energy dispersion cube energy if the energy is
910   - // outside the observation energy. The energy dispersion cube
911   - // energies are true energies while the observation energy
912   - // boundaries are reconstructed energies, hence this is only
913   - // an approximation, but probably the only we can really do.
914   - if (!obs_ebounds.contains(m_energies[iebin])) {
915   - continue;
916   - }
917 906  
918   - // Get log10 of true energy in TeV
919   - double logEsrc = m_energies[iebin].log10TeV();
  907 + // Loop over all pixels in sky map
  908 + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
920 909  
921   - // Compute exposure weight
922   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
923   - obs.livetime();
  910 + // Get pixel sky direction
  911 + GSkyDir dir = m_cube.inx2dir(pixel);
924 912  
925   - // If available, accumulate weights
926   - if (exposure != NULL) {
927   - (*exposure)(pixel, iebin) += weight;
  913 + // Skip pixel if it is outside the RoI
  914 + if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
  915 + continue;
928 916 }
  917 +
  918 + // Compute theta angle with respect to pointing direction in
  919 + // radians
  920 + double theta = pnt.dist(dir);
  921 +
  922 + // Loop over all energy bins
  923 + for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
  924 +
  925 + // Skip energy dispersion cube energy if the energy is
  926 + // outside the observation energy. The energy dispersion
  927 + // cube energies are true energies while the observation
  928 + // energy boundaries are reconstructed energies, hence
  929 + // this is only an approximation, but probably the only
  930 + // we can really do.
  931 + if (!obs_ebounds.contains(m_energies[iebin])) {
  932 + continue;
  933 + }
929 934  
930   - // Loop over delta values
931   - for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
  935 + // Get log10 of true energy in TeV
  936 + double logEsrc = m_energies[iebin].log10TeV();
932 937  
933   - // Compute energy migration fraction
934   - double migra = m_migras[imigra];
  938 + // Compute exposure weight
  939 + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
  940 + obs.livetime();
935 941  
936   - // Skip migra bin if zero
937   - if (migra <= 0.0) {
938   - continue;
  942 + // If available, accumulate weights
  943 + if (exposure != NULL) {
  944 + (*exposure)(pixel, iebin) += weight;
939 945 }
940 946  
941   - // Compute log10 of reconstructed energy in TeV
942   - double logEobs = std::log10(migra) + logEsrc;
943   - GEnergy Eobs(std::pow(10.0, logEobs), "TeV");
  947 + // Loop over delta values
  948 + for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
944 949  
945   - // Set map index
946   - int imap = offset(imigra, iebin);
  950 + // Compute energy migration fraction
  951 + double migra = m_migras[imigra];
  952 +
  953 + // Skip migra bin if zero
  954 + if (migra <= 0.0) {
  955 + continue;
  956 + }
  957 +
  958 + // Compute log10 of reconstructed energy in TeV
  959 + double logEobs = std::log10(migra) + logEsrc;
  960 + GEnergy Eobs(std::pow(10.0, logEobs), "TeV");
947 961  
948   - // Add energy dispersion cube value
949   - m_cube(pixel, imap) += rsp->edisp(Eobs,
950   - theta, 0.0,
951   - 0.0, 0.0,
952   - logEsrc) * weight;
  962 + // Set map index
  963 + int imap = offset(imigra, iebin);
953 964  
954   - } // endfor: looped over migration bins
  965 + // Add energy dispersion cube value
  966 + m_cube(pixel, imap) += rsp->edisp(Eobs,
  967 + theta, 0.0,
  968 + 0.0, 0.0,
  969 + logEsrc) * weight;
  970 +
  971 + } // endfor: looped over migration bins
955 972  
956   - } // endfor: looped over energy bins
  973 + } // endfor: looped over energy bins
957 974  
958   - } // endfor: looped over all pixels
  975 + } // endfor: looped over all pixels
  976 +
  977 + } // endelse: livetime was not zero
959 978  
960 979 } // endif: observation contained an event list
961 980  
... ...
inst/cta/src/GCTACubeExposure.cpp
  View file @c616af2
1 1 /***************************************************************************
2 2 * GCTACubeExposure.cpp - CTA cube analysis exposure class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2017 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -656,56 +656,74 @@ void GCTACubeExposure::fill_cube(const GCTAObservation&amp; obs, GLog* log)
656 656 throw GException::invalid_value(G_FILL_CUBE, msg);
657 657 }
658 658  
659   - // Announce observation usage
660   - if (log != NULL) {
661   - *log << "Including ";
662   - *log << obs.instrument();
663   - *log << " observation \"" << obs.name();
664   - *log << "\" (id=" << obs.id() << ")";
665   - *log << " in exposure cube computation." << std::endl;
  659 + // Skip observation if livetime is zero
  660 + if (obs.livetime() == 0.0) {
  661 + if (log != NULL) {
  662 + *log << "Skipping unbinned ";
  663 + *log << obs.instrument();
  664 + *log << " observation ";
  665 + *log << "\"" << obs.name() << "\"";
  666 + *log << " (id=" << obs.id() << ") due to zero livetime";
  667 + *log << std::endl;
  668 + }
666 669 }
667 670  
668   - // Loop over all pixels in sky map
669   - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
670   -
671   - // Get pixel sky direction
672   - GSkyDir dir = m_cube.inx2dir(pixel);
  671 + // ... otherwise continue
  672 + else {
673 673  
674   - // Skip pixel if it is outside the RoI
675   - if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
676   - continue;
  674 + // Announce observation usage
  675 + if (log != NULL) {
  676 + *log << "Including ";
  677 + *log << obs.instrument();
  678 + *log << " observation \"" << obs.name();
  679 + *log << "\" (id=" << obs.id() << ")";
  680 + *log << " in exposure cube computation." << std::endl;
677 681 }
678   -
679   - // Compute theta angle with respect to pointing direction in radians
680   - double theta = pnt.dist(dir);
681   -
682   - // Loop over all exposure cube energies
683   - for (int iebin = 0; iebin < m_energies.size(); ++iebin){
684   -
685   - // Skip exposure cube energy if the energy is outside the
686   - // observation energy. The exposure cube energies are true
687   - // energies while the observation energy boundaries are
688   - // reconstructed energies, hence this is only an approximation,
689   - // but probably the only we can really do.
690   - if (!obs_ebounds.contains(m_energies[iebin])) {
  682 +
  683 + // Loop over all pixels in sky map
  684 + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
  685 +
  686 + // Get pixel sky direction
  687 + GSkyDir dir = m_cube.inx2dir(pixel);
  688 +
  689 + // Skip pixel if it is outside the RoI
  690 + if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
691 691 continue;
692 692 }
  693 +
  694 + // Compute theta angle with respect to pointing direction in
  695 + // radians
  696 + double theta = pnt.dist(dir);
693 697  
694   - // Get logE/TeV
695   - double logE = m_energies[iebin].log10TeV();
  698 + // Loop over all exposure cube energies
  699 + for (int iebin = 0; iebin < m_energies.size(); ++iebin){
696 700  
697   - // Add to exposure cube (effective area * livetime)
698   - m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
699   - obs.livetime();
  701 + // Skip exposure cube energy if the energy is outside the
  702 + // observation energy. The exposure cube energies are true
  703 + // energies while the observation energy boundaries are
  704 + // reconstructed energies, hence this is only an
  705 + // approximation, but probably the only we can really do.
  706 + if (!obs_ebounds.contains(m_energies[iebin])) {
  707 + continue;
  708 + }
700 709  
701   - } // endfor: looped over energy bins
  710 + // Get logE/TeV
  711 + double logE = m_energies[iebin].log10TeV();
702 712  
703   - } // endfor: looped over all pixels
  713 + // Add to exposure cube (effective area * livetime)
  714 + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  715 + obs.livetime();
  716 +
  717 + } // endfor: looped over energy bins
  718 +
  719 + } // endfor: looped over all pixels
  720 +
  721 + // Append GTIs and increment livetime
  722 + m_gti.extend(obs.gti());
  723 + m_livetime += obs.livetime();
  724 +
  725 + } // endelse: livetime was not zero
704 726  
705   - // Append GTIs and increment livetime
706   - m_gti.extend(obs.gti());
707   - m_livetime += obs.livetime();
708   -
709 727 } // endif: observation contained an event list
710 728  
711 729 // Return
... ...
inst/cta/src/GCTACubePsf.cpp
  View file @c616af2
1 1 /***************************************************************************
2 2 * GCTACubePsf.cpp - CTA cube analysis point spread function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2017 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -787,71 +787,89 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
787 787 throw GException::invalid_value(G_FILL_CUBE, msg);
788 788 }
789 789  
790   - // Announce observation usage
791   - if (log != NULL) {
792   - *log << "Including ";
793   - *log << obs.instrument();
794   - *log << " observation \"" << obs.name();
795   - *log << "\" (id=" << obs.id() << ")";
796   - *log << " in point spread function cube computation." << std::endl;
  790 + // Skip observation if livetime is zero
  791 + if (obs.livetime() == 0.0) {
  792 + if (log != NULL) {
  793 + *log << "Skipping unbinned ";
  794 + *log << obs.instrument();
  795 + *log << " observation ";
  796 + *log << "\"" << obs.name() << "\"";
  797 + *log << " (id=" << obs.id() << ") due to zero livetime";
  798 + *log << std::endl;
  799 + }
797 800 }
798 801  
799   - // Loop over all pixels in sky map
800   - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
  802 + // ... otherwise continue
  803 + else {
801 804  
802   - // Get pixel sky direction
803   - GSkyDir dir = m_cube.inx2dir(pixel);
804   -
805   - // Skip pixel if it is outside the RoI
806   - if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
807   - continue;
  805 + // Announce observation usage
  806 + if (log != NULL) {
  807 + *log << "Including ";
  808 + *log << obs.instrument();
  809 + *log << " observation \"" << obs.name();
  810 + *log << "\" (id=" << obs.id() << ")";
  811 + *log << " in point spread function cube computation." << std::endl;
808 812 }
809 813  
810   - // Compute theta angle with respect to pointing direction in radians
811   - double theta = pnt.dist(dir);
812   -
813   - // Loop over all energy bins
814   - for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
  814 + // Loop over all pixels in sky map
  815 + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
815 816  
816   - // Skip PSF cube energy if the energy is outside the observation
817   - // energy. The PSF cube energies are true energies while the
818   - // observation energy boundaries are reconstructed energies,
819   - // hence this is only an approximation, but probably the only
820   - // we can really do.
821   - if (!obs_ebounds.contains(m_energies[iebin])) {
  817 + // Get pixel sky direction
  818 + GSkyDir dir = m_cube.inx2dir(pixel);
  819 +
  820 + // Skip pixel if it is outside the RoI
  821 + if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
822 822 continue;
823 823 }
824 824  
825   - // Get logE/TeV
826   - double logE = m_energies[iebin].log10TeV();
  825 + // Compute theta angle with respect to pointing direction in
  826 + // radians
  827 + double theta = pnt.dist(dir);
827 828  
828   - // Compute exposure weight
829   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
830   - obs.livetime();
  829 + // Loop over all energy bins
  830 + for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
831 831  
832   - // If available, accumulate weights
833   - if (exposure != NULL) {
834   - (*exposure)(pixel, iebin) += weight;
835   - }
  832 + // Skip PSF cube energy if the energy is outside the
  833 + // observation energy. The PSF cube energies are true
  834 + // energies while the observation energy boundaries are
  835 + // reconstructed energies, hence this is only an
  836 + // approximation, but probably the only we can really do.
  837 + if (!obs_ebounds.contains(m_energies[iebin])) {
  838 + continue;
  839 + }
836 840  
837   - // Loop over delta values
838   - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
  841 + // Get logE/TeV
  842 + double logE = m_energies[iebin].log10TeV();
839 843  
840   - // Compute delta in radians
841   - double delta = m_deltas[idelta] * gammalib::deg2rad;
  844 + // Compute exposure weight
  845 + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  846 + obs.livetime();
842 847  
843   - // Set map index
844   - int imap = offset(idelta, iebin);
  848 + // If available, accumulate weights
  849 + if (exposure != NULL) {
  850 + (*exposure)(pixel, iebin) += weight;
  851 + }
845 852  
846   - // Add on PSF cube
847   - m_cube(pixel, imap) +=
848   - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
  853 + // Loop over delta values
  854 + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
849 855  
850   - } // endfor: looped over delta bins
  856 + // Compute delta in radians
  857 + double delta = m_deltas[idelta] * gammalib::deg2rad;
851 858  
852   - } // endfor: looped over energy bins
  859 + // Set map index
  860 + int imap = offset(idelta, iebin);
853 861  
854   - } // endfor: looped over all pixels
  862 + // Add on PSF cube
  863 + m_cube(pixel, imap) +=
  864 + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
  865 +
  866 + } // endfor: looped over delta bins
  867 +
  868 + } // endfor: looped over energy bins
  869 +
  870 + } // endfor: looped over all pixels
  871 +
  872 + } // endelse: livetime was not zero
855 873  
856 874 } // endif: observation contained an event list
857 875  
... ...