35 #define G_GET_PARAMETERS "ctskymap::get_parameter()"
36 #define G_FILL_MAPS_COUNTS "ctskymap::fill_maps_counts(GCTAObservation*)"
37 #define G_FILL_MAPS_ACCEPTANCE "ctskymap::fill_maps_acceptance("\
39 #define G_RING_BOUNDING_BOX "ctskymap::ring_bounding_box(int&, int&, int&, "\
41 #define G_RING_KERNEL "ctskymap::ring_kernel(double&, double&)"
194 this->GApplication::clear();
255 log_header1(TERSE,
"Save sky map");
258 if ((*
this)[
"outmap"].is_valid() && !
m_skymap.is_empty()) {
261 m_outmap = (*this)[
"outmap"].filename();
271 fname.append(
" (map is empty, no file created)");
273 log_value(NORMAL,
"Sky map file", fname);
336 log_header1(TERSE,
"Publish sky map");
339 std::string user_name(name);
340 if (user_name.empty()) {
345 log_value(NORMAL,
"Sky map name", user_name);
473 if ((*
this)[
"inmap"].is_valid()) {
476 GFilename inmap = (*this)[
"inmap"].filename();
482 if (fits.contains(
"COUNTS") && fits.contains(
"ACCEPTANCE")) {
485 GFitsHDU &counts = *fits[
"COUNTS"];
486 GFitsHDU &acceptance = *fits[
"ACCEPTANCE"];
496 m_emin = (counts.has_card(
"E_MIN")) ? counts.real(
"E_MIN") : 0.0;
497 m_emax = (counts.has_card(
"E_MAX")) ? counts.real(
"E_MAX") : 0.0;
503 std::string msg =
"Input sky map \""+inmap.url()+
"\" does not "
504 "contain \"COUNTS\" and \"ACCEPTANCE\" "
523 m_emin = (*this)[
"emin"].real();
524 m_emax = (*this)[
"emax"].real();
543 (*this)[
"inexclusion"].query();
544 m_usefft = (*this)[
"usefft"].boolean();
548 std::string msg(
"\"roiradius\" must be smaller than \"inradius\"");
552 std::string msg(
"\"inradius\" must be smaller than \"outradius\"");
565 m_publish = (*this)[
"publish"].boolean();
566 m_chatter =
static_cast<GChatter
>((*this)[
"chatter"].integer());
570 (*this)[
"outmap"].query();
607 for (
int i = 0; i <
m_counts.npix(); ++i) {
642 if ((*
this)[
"inexclusion"].is_valid()) {
694 GSkyRegions regions(filename);
697 for (
int i = 0; i <
m_exclmap.npix(); ++i) {
703 if (regions.contains(dir)) {
730 for (
int i = 0; i < dest_map.npix(); ++i) {
733 GSkyDir dir = dest_map.pix2dir(i);
765 log_header1(TERSE, gammalib::number(
"Find unbinned observation",
769 std::vector<GCTAObservation*> obs_list(0);
774 obs_list.push_back(
obs);
777 std::string msg =
" Including unbinned "+
obs->instrument()+
779 log_string(NORMAL, msg);
783 log_header1(TERSE, gammalib::number(
"Fill map from observation",
787 #pragma omp parallel for
788 for (
int i = 0; i < obs_list.size(); ++i) {
791 GCTAObservation*
obs = obs_list[i];
802 obs->dispose_events();
825 GCTAEventList* events =
dynamic_cast<GCTAEventList*
>
826 (
const_cast<GEvents*
>(obs->events()));
830 if (events == NULL) {
831 std::string msg =
"Specified observation does not contain a CTA event "
832 "list. Please specify a CTA observation containing "
833 "a CTA event list as argument.";
845 int num_outside_roi = 0;
846 int num_outside_map = 0;
847 int num_outside_erange = 0;
851 GCTARoi roi = obs->roi();
854 for (
int i = 0; i < events->size(); ++i) {
857 GCTAEventAtom*
event = (*events)[i];
860 if (event->energy() < emin ||
event->energy() > emax) {
861 num_outside_erange++;
866 GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
867 GSkyDir dir = inst->dir();
868 GSkyPixel pixel =
m_counts.dir2pix(dir);
871 if (pixel.x() < -0.5 || pixel.x() > (
m_counts.nx() - 0.5) ||
872 pixel.y() < -0.5 || pixel.y() > (
m_counts.ny() - 0.5)) {
878 if (roi.is_valid() && !roi.contains(*inst)) {
884 #pragma omp critical(ctskymap_fill_maps_counts_1)
891 #pragma omp critical(ctskymap_fill_maps_counts_2)
894 log_value(NORMAL,
"Events in list", obs->events()->size());
895 log_value(NORMAL,
"Events in map", num_in_map);
896 log_value(NORMAL,
"Events outside RoI", num_outside_roi);
897 log_value(NORMAL,
"Events outside map area", num_outside_map);
898 log_value(NORMAL,
"Events outside energies", num_outside_erange);
901 log_header1(EXPLICIT,
"Sky map");
933 const GCTABackground* bkg = NULL;
936 const GCTAResponseIrf* rsp =
dynamic_cast<const GCTAResponseIrf*
>
941 std::string msg =
"No response information available for "+
943 "Please specify response information or use "
944 "another background subtraction method.";
952 bkg = rsp->background();
956 std::string msg =
"No IRF background template found in instrument "
957 "response function for "+
959 "instrument response function containing a "
960 "background template.";
967 GEnergy emin(
m_emin,
"TeV");
968 GEnergy emax(
m_emax,
"TeV");
975 GSkyDir roi_centre = obs->roi().centre().dir();
976 double cos_roi_radius = std::cos(obs->roi().radius() * gammalib::deg2rad);
979 double exposure = obs->livetime();
988 GSkyDir& skydir =
m_dirs[i];
991 if (roi_centre.cos_dist(skydir) < cos_roi_radius) {
996 GCTAInstDir instdir = obs->pointing().instdir(skydir);
999 double value = (bkg != NULL) ? bkg->rate_ebin(instdir, emin, emax) : 1.0;
1005 #pragma omp critical(ctskymap_fill_maps_acceptance_1)
1014 #pragma omp critical(ctskymap_fill_maps_acceptance_2)
1017 std::string msg = (bkg != NULL) ?
"Use IRF background template"
1018 :
"Use constant background rate";
1021 log_value(NORMAL,
"Background estimate", msg);
1022 log_value(NORMAL,
"Events in background",
int(total+0.5));
1073 log_header1(TERSE,
"Compute maps");
1107 for (
int i = 0; i <
m_exclmap.npix(); ++i) {
1140 log_header3(NORMAL,
"Computing ring background map (FFT method)");
1151 int num_bad_alpha = 0;
1157 for (
int i = 0; i <
m_counts.npix(); ++i) {
1159 excl_counts(i) = 0.0;
1160 excl_acceptance(i) = 0.0;
1171 for (
int i = 0; i <
m_sigmap.npix(); ++i) {
1174 double n_on = on_counts(i);
1175 double n_off = off_counts(i);
1177 if (off_alpha(i) == 0.0) {
1183 alpha = on_alpha(i) / off_alpha(i);
1203 log_value(NORMAL,
"Total number of pixels",
m_sigmap.npix());
1204 log_value(NORMAL,
"Pixels with alpha=0", num_bad_alpha);
1226 log_header3(NORMAL,
"Computing ring background map (direct method)");
1227 log_value(NORMAL,
"Total pixels to process",
m_counts.npix());
1238 int num_bad_alpha = 0;
1243 int n_logpix =
m_counts.npix() / 10;
1244 if (n_logpix > 100000) {
1249 #pragma omp parallel for
1250 for (
int i = 0; i <
m_sigmap.npix(); ++i) {
1258 if (i % n_logpix == 0 && n_logpix > 10000) {
1259 #pragma omp critical(ctskymap_map_significance_ring)
1260 log_value(NORMAL,
"Pixels remaining",
m_sigmap.npix()-i);
1273 #pragma omp critical(ctskymap_compute_maps_ring_direct)
1282 #pragma omp critical(ctskymap_compute_maps_ring_direct)
1290 log_value(NORMAL,
"Pixels with alpha=0", num_bad_alpha);
1324 const GSkyMap& counts,
1325 const GSkyMap& background,
1336 double alpha_on = 0.0;
1337 double alpha_off = 0.0;
1340 GSkyDir& position =
m_dirs[ipixel];
1350 int nx = counts.nx();
1353 for (
int ix_nowrap = ix_start; ix_nowrap < ix_stop; ++ix_nowrap) {
1361 else if (ix >= nx) {
1366 int i = ix + iy_start * nx;
1369 for (
int iy = iy_start; iy < iy_stop; ++iy, i += nx) {
1372 GSkyDir& skydir =
m_dirs[i];
1386 alpha_off += background(i);
1398 alpha_on += background(i);
1410 if (alpha_off == 0.0) {
1416 alpha = alpha_on / alpha_off;
1445 int& iy1,
int& iy2)
const
1452 int ix0 = ipixel % nx;
1453 int iy0 = ipixel / nx;
1456 const GWcs* wcs =
dynamic_cast<const GWcs*
>(
m_counts.projection());
1458 std::string msg =
"Sky map is not a WCS projection. Method is only "
1459 "valid for WCS projections.";
1464 double dx_binsz = wcs->cdelt(0);
1465 double dy_binsz = wcs->cdelt(1);
1469 int ix = (ix0 > 0) ? ix0 - 1 : ix0 + 1;
1473 if (dx < dx_binsz) {
1479 int iy = (iy0 > 0) ? iy0 - 1 : iy0 + 1;
1483 if (dy < dy_binsz) {
1495 if (2*nx_bbox < nx) {
1496 ix1 = ix0 - nx_bbox;
1497 ix2 = ix0 + nx_bbox;
1505 if (2*ny_bbox < ny) {
1506 iy1 = iy0 - ny_bbox;
1507 iy2 = iy0 + ny_bbox;
1536 const double& rmax)
const
1539 GSkyMap conv_map = map;
1545 GNdarray array(conv_map.nx(), conv_map.ny());
1546 double *dst = array.data();
1547 for (
int i = 0; i < conv_map.npix(); ++i) {
1548 *dst++ = conv_map(i);
1552 GFft fft_array = GFft(array);
1555 GFft fft_smooth = fft_array * fft_kernel;
1558 GNdarray smooth = fft_smooth.backward();
1561 double *src = smooth.data();
1562 for (
int i = 0; i < conv_map.npix(); ++i) {
1563 conv_map(i) = *src++;
1589 const GWcs* wcs =
dynamic_cast<const GWcs*
>(
m_counts.projection());
1591 std::string msg =
"Sky map is not a WCS projection. Method is only "
1592 "valid for WCS projections.";
1597 double dx = wcs->cdelt(0);
1598 double dy = wcs->cdelt(1);
1604 for (
int ix1 = 0, ix2 =
m_counts.nx(); ix1 <
m_counts.nx(); ++ix1, --ix2) {
1605 double x = ix1 * dx;
1607 for (
int iy1 = 0, iy2 =
m_counts.ny(); iy1 <
m_counts.ny(); ++iy1, --iy2) {
1608 double y = iy1 * dy;
1609 double r = std::sqrt(xqs + y*y);
1610 if ((r >= rmin) && (r <= rmax)) {
1611 kern(ix1,iy1) += 1.0;
1613 kern(ix2,iy1) += 1.0;
1616 kern(ix1,iy2) += 1.0;
1619 kern(ix2,iy2) += 1.0;
1641 const double& n_off,
1642 const double& alpha)
const
1649 sigma = -2.0 * n_off * std::log(1.0+alpha);
1653 else if (n_off == 0.0) {
1654 sigma = 2.0 * n_on * std::log((1.0+alpha)/alpha);
1659 sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
1660 (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
1661 n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
1681 GFitsHDU* hdu = map.write(fits);
1685 hdu->extname(extname);
1710 hdu->card(
"BKGSUB",
m_bkgsubtract,
"Background substraction method");
1711 hdu->card(
"E_MIN",
m_emin,
"[TeV] Lower energy boundary");
1712 hdu->card(
"E_MAX",
m_emax,
"[TeV] Upper energy boundary");
1713 hdu->card(
"EUNIT",
"TeV",
"Units for E_MIN and E_MAX");
1717 hdu->card(
"USEFFT",
m_usefft,
"Use FFT for RING background");
1718 hdu->card(
"ROIRAD",
m_roiradius,
"[deg] Source region radius");
1719 hdu->card(
"INRAD",
m_inradius,
"[deg] Inner background ring radius");
1720 hdu->card(
"OUTRAD",
m_outradius,
"[deg] Outer background ring radius");
1721 hdu->card(
"ITER",
m_iterations,
"Exclusion map iterations");
1722 hdu->card(
"THRES",
m_threshold,
"Exclusion map threshold");
1723 hdu->card(
"EXCLMAP",
m_inexclusion.url(),
"Exclusion map name");
double m_roiradius
Region of interest radius for RING bkg.
void setup_exclusion_map(void)
Generates map of pixel exclusions.
void setup_maps(void)
Setup maps.
void clear(void)
Clear sky mapping tool.
GNdarray ring_kernel(const double &rmin, const double &rmax) const
Return FFT kernel for background ring.
int m_iterations
Number of iterations for RING background.
#define G_FILL_MAPS_ACCEPTANCE
GFilename m_outmap
Output file name.
GFits m_fits
Output GFits object.
void setup_exclusion_map_region(const GFilename &filename)
Fills exclusions map from DS9 region file.
const GFits & fits(void) const
Return fits container.
bool m_publish
Publish sky map?
GFilename m_inexclusion
Exclusion map file name.
double m_threshold
Threshold for RING background.
double m_inradius
Inner ring radius for RING background.
bool m_usefft
Use FFT for RING background.
void ring_bounding_box(const int &ipixel, int &ix1, int &ix2, int &iy1, int &iy2) const
Computes bounding box for RING background computation.
void free_members(void)
Delete class members.
const GObservations & obs(void) const
Return observation container.
GSkyMap m_bkgmap
Background map.
void fill_maps_acceptance(GCTAObservation *obs)
Compute background acceptance sky map based on IRF template.
void setup_exclusion_map_fits(const GFilename &filename)
Fills exclusions map from FITS image.
GCTAObservation * next_unbinned_observation(void)
Return next unbinned CTA observation.
GSkyMap m_counts
Counts map.
virtual ~ctskymap(void)
Destructor.
GSkyMap ring_convolve(const GSkyMap &map, const double &rmin, const double &rmax) const
Return FFT kernel for background ring.
std::string m_bkgsubtract
Background subtraction method.
void fill_maps(void)
Fill maps from observation container.
void copy_members(const ctskymap &app)
Copy class members.
double m_cos_outradius
Cosine of outer ring radius.
void write_map(GFits &fits, const GSkyMap &map, const std::string &extname) const
Write sky map into FITS file.
GChatter m_chatter
Chattiness.
Base class for observation tools.
void read_ogip_keywords(GFitsHDU *hdu) const
Read OGIP keywords from FITS HDU.
void compute_ring_values(const int &ipixel, const GSkyMap &counts, const GSkyMap &background, double &non, double &noff, double &alpha)
Computes Non, Noff and alpha for a counts map and sensitivity map.
void write_hdu_keywords(GFitsHDU *hdu) const
Write keywords in FITS HDU.
double m_cos_roiradius
Cosine of RoI radius.
void free_members(void)
Delete class members.
#define G_FILL_MAPS_COUNTS
std::vector< GSkyDir > m_dirs
Cached pixel directions.
void get_parameters(void)
Get application parameters.
void save(void)
Save sky map.
double sigma_li_ma(const double &n_on, const double &n_off, const double &alpha) const
Compute significance following Li & Ma.
void init_members(void)
Initialise class members.
void write_ogip_keywords(GFitsHDU *hdu) const
Write OGIP keywords in FITS HDU.
void process(void)
Process the sky mapping tool.
void fill_maps_counts(GCTAObservation *obs)
Fill events into counts map.
double m_outradius
Outer ring radius for RING background.
void init_members(void)
Initialise class members.
std::vector< double > m_solidangle
Cached pixel solid angles.
GCTAObservation * first_unbinned_observation(void)
Return first unbinned CTA observation.
GSkyMap m_acceptance
Acceptance map.
void compute_maps(void)
Compute sky map, background map and significance map.
ctobservation & operator=(const ctobservation &app)
Assignment operator.
double m_emax
Maximum energy (TeV)
ctskymap(void)
Void constructor.
double m_emin
Minimum energy (TeV)
bool m_has_inmap
Has valid input map.
void compute_maps_ring_fft(void)
Compute the maps for RING background using a FFT.
double m_cos_inradius
Cosine of inner ring radius.
void publish(const std::string &name="")
Publish sky map.
GObservations m_obs
Observation container.
void adjust_exclusion_map(void)
Verifys that exclusion map points in fov of counts.
void compute_maps_ring_direct(void)
Compute the pixel significance for RING background.
GSkyMap m_exclmap
Exclusion map for RING background.
ctskymap & operator=(const ctskymap &app)
Assignment operator.
GSkyMap m_sigmap
Significance map.
Sky mapping tool definition.
void construct_fits(void)
Construct GFits object consisting of all maps.
#define G_RING_BOUNDING_BOX