Olga Sergijenko / ctools

Commit 217ed01922444c6dfc06a57bc2e6e7ae0c864da1

Authored by Jürgen Knödlseder
1 parent bc051b10

Speed up IRF computation by using GSkyDir::cos_dist() to select relevant pixels.… 

… Also implement OpenMP support for multiple observations.
Inline Side-by-side
Showing 4 changed files with 343 additions and 365 deletions
doc/source/users/reference_manual/ctskymap.rst
  View file @217ed01
... ... @@ -82,13 +82,13 @@ Background subtraction configuration parameters
82 82 Background subtraction method.
83 83  
84 84 ``roiradius [real]``
85   - Source region radius for estimating on-counts (in degrees).
  85 + Source region radius for ``RING`` subtraction (in degrees).
86 86  
87 87 ``inradius [real]``
88   - Inner ring radius for RING subtraction (in degrees).
  88 + Inner background ring radius for ``RING`` subtraction (in degrees).
89 89  
90 90 ``outradius [real]``
91   - Outer ring radius for RING subtraction (in degrees).
  91 + Outer background ring radius for ``RING`` subtraction (in degrees).
92 92  
93 93 ``(inexclusion = NONE) [file]``
94 94 Exclusion region file as either a FITS map or DS9 region file.
... ...
src/ctskymap/ctskymap.cpp
  View file @217ed01
... ... @@ -35,8 +35,6 @@
35 35 #define G_GET_PARAMETERS "ctskymap::get_parameter()"
36 36 #define G_MAP_EVENTS "ctskymap::map_events(GCTAObservation*)"
37 37 #define G_MAP_BACKGROUND_IRF "ctskymap::map_background_irf(GCTAObservation*)"
38   -#define G_MAP_BACKGROUND_RING "ctskymap::map_background_ring("\
39   - "GCTAObservation*)"
40 38  
41 39 /* __ Debug definitions __________________________________________________ */
42 40  
... ... @@ -221,20 +219,44 @@ void ctskymap::run(void)
221 219 // Get parameters
222 220 get_parameters();
223 221  
  222 + // Setup maps
  223 + setup_maps();
  224 +
224 225 // Write input observation container into logger
225 226 log_observations(NORMAL, m_obs, "Input observation");
226 227  
227 228 // Write header into logger
  229 + log_header1(TERSE, gammalib::number("Find unbinned observation",
  230 + m_obs.size()));
  231 +
  232 + // Find all unbinned CTA observations in m_obs
  233 + std::vector<GCTAObservation*> obs_list(0);
  234 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  235 + obs = next_unbinned_observation()) {
  236 +
  237 + // Push observation into list
  238 + obs_list.push_back(obs);
  239 +
  240 + // Write message
  241 + std::string msg = " Including unbinned "+obs->instrument()+
  242 + " observation";
  243 + log_string(NORMAL, msg);
  244 + }
  245 +
  246 + // Write header into logger
228 247 log_header1(TERSE, gammalib::number("Map observation", m_obs.size()));
229 248  
230 249 // Loop over all unbinned CTA observations in the container
231   - for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
232   - obs = next_unbinned_observation()) {
  250 + #pragma omp parallel for
  251 + for (int i = 0; i < obs_list.size(); ++i) {
  252 +
  253 + // Get pointer to observation
  254 + GCTAObservation* obs = obs_list[i];
233 255  
234 256 // Map events into sky map
235 257 map_events(obs);
236 258  
237   - // Compute background sky map
  259 + // Compute background sky map (only if m_bkgsubtract != "NONE")
238 260 map_background(obs);
239 261  
240 262 // Dispose events to free memory
... ... @@ -250,7 +272,8 @@ void ctskymap::run(void)
250 272  
251 273 // Subtract background map from counts map
252 274 m_skymap -= m_bkgmap;
253   - }
  275 +
  276 + } // endif: background subtraction is selected
254 277  
255 278 // Optionally publish sky map
256 279 if (m_publish) {
... ... @@ -292,7 +315,7 @@ void ctskymap::save(void)
292 315 // If background subtraction is requested then write background map
293 316 // and significance map to FITS file
294 317 if (m_bkgsubtract != "NONE") {
295   -
  318 +
296 319 // Write background map into FITS file
297 320 hdu = m_bkgmap.write(fits);
298 321  
... ... @@ -375,12 +398,8 @@ void ctskymap::publish(const std::string&amp; name)
375 398 void ctskymap::init_members(void)
376 399 {
377 400 // Initialise members
378   - m_skymap.clear();
379   - m_bkgmap.clear();
380   - m_sigmap.clear();
381   - m_exclmap.clear();
382   - m_alphamap.clear();
383   - m_onmap.clear();
  401 + m_outmap.clear();
  402 + m_inexclusion.clear();
384 403 m_emin = 0.0;
385 404 m_emax = 0.0;
386 405 m_bkgsubtract = "NONE";
... ... @@ -390,6 +409,12 @@ void ctskymap::init_members(void)
390 409 m_publish = false;
391 410 m_chatter = static_cast<GChatter>(2);
392 411  
  412 + // Initialise maps
  413 + m_skymap.clear();
  414 + m_bkgmap.clear();
  415 + m_sigmap.clear();
  416 + m_exclmap.clear();
  417 +
393 418 // Initialise cache
394 419 m_solidangle.clear();
395 420 m_dirs.clear();
... ... @@ -406,13 +431,9 @@ void ctskymap::init_members(void)
406 431 ***************************************************************************/
407 432 void ctskymap::copy_members(const ctskymap& app)
408 433 {
409   - // Copy attributes
410   - m_skymap = app.m_skymap;
411   - m_bkgmap = app.m_bkgmap;
412   - m_sigmap = app.m_sigmap;
413   - m_exclmap = app.m_exclmap;
414   - m_alphamap = app.m_alphamap;
415   - m_onmap = app.m_onmap;
  434 + // Copy members
  435 + m_outmap = app.m_outmap;
  436 + m_inexclusion = app.m_inexclusion;
416 437 m_emin = app.m_emin;
417 438 m_emax = app.m_emax;
418 439 m_bkgsubtract = app.m_bkgsubtract;
... ... @@ -422,9 +443,15 @@ void ctskymap::copy_members(const ctskymap&amp; app)
422 443 m_publish = app.m_publish;
423 444 m_chatter = app.m_chatter;
424 445  
  446 + // Copy maps
  447 + m_skymap = app.m_skymap;
  448 + m_bkgmap = app.m_bkgmap;
  449 + m_sigmap = app.m_sigmap;
  450 + m_exclmap = app.m_exclmap;
  451 +
425 452 // Copy cache
426   - m_solidangle = app.m_solidangle;
427   - m_dirs = app.m_dirs;
  453 + m_solidangle = app.m_solidangle;
  454 + m_dirs = app.m_dirs;
428 455  
429 456 // Return
430 457 return;
... ... @@ -456,7 +483,7 @@ void ctskymap::get_parameters(void)
456 483 // information and do not accept counts cubes.
457 484 setup_observations(m_obs, false, true, false);
458 485  
459   - // Create sky map based on task parameters
  486 + // Create sky map. Note that the method reads the relevant User parameters.
460 487 m_skymap = create_map(m_obs);
461 488  
462 489 // Get further parameters
... ... @@ -468,9 +495,10 @@ void ctskymap::get_parameters(void)
468 495 if (m_bkgsubtract == "RING") {
469 496  
470 497 // Get parameters
471   - m_roiradius = (*this)["roiradius"].real();
472   - m_inradius = (*this)["inradius"].real();
473   - m_outradius = (*this)["outradius"].real();
  498 + m_roiradius = (*this)["roiradius"].real();
  499 + m_inradius = (*this)["inradius"].real();
  500 + m_outradius = (*this)["outradius"].real();
  501 + m_inexclusion = (*this)["inexclusion"].filename();
474 502  
475 503 // Make sure that (roiradius < inradius < outradius)
476 504 if (m_roiradius > m_inradius) {
... ... @@ -484,8 +512,8 @@ void ctskymap::get_parameters(void)
484 512  
485 513 } // endif: read RING background parameters
486 514  
487   - // If IRF background subtraction is requested then make sure that the
488   - // CTA observations in the observation container have response information
  515 + // If background subtraction is requested then make sure that the CTA
  516 + // observations in the observation container have response information
489 517 if (m_bkgsubtract != "NONE") {
490 518 set_response(m_obs);
491 519 }
... ... @@ -496,19 +524,38 @@ void ctskymap::get_parameters(void)
496 524  
497 525 // Read ahead parameters
498 526 if (read_ahead()) {
499   - m_outmap = (*this)["outmap"].filename();
  527 + m_outmap = (*this)["outmap"].filename();
500 528 }
501 529  
  530 + // Write parameters into logger
  531 + log_parameters(TERSE);
  532 +
  533 + // Return
  534 + return;
  535 +}
  536 +
  537 +
  538 +/***********************************************************************//**
  539 + * @brief Setup maps
  540 + *
  541 + * Setup maps for sky map generation.
  542 + ***************************************************************************/
  543 +void ctskymap::setup_maps(void)
  544 +{
  545 + // Clear cache
  546 + m_solidangle.clear();
  547 + m_dirs.clear();
  548 +
502 549 // Create background map and significance map if background subtraction
503 550 // is requested
504 551 if (m_bkgsubtract != "NONE") {
505 552  
506   - // Create backgrond and significance maps
  553 + // Create background and significance maps
507 554 m_bkgmap = create_map(m_obs);
508 555 m_sigmap = create_map(m_obs);
509   -
  556 +
510 557 // Setup the exclusions map
511   - map_exclusions((*this)["inexclusion"].filename());
  558 + map_exclusions(m_inexclusion);
512 559  
513 560 // Cache the pixel solid angles and sky directions
514 561 m_solidangle.reserve(m_bkgmap.npix());
... ... @@ -518,17 +565,107 @@ void ctskymap::get_parameters(void)
518 565 m_dirs.push_back(m_bkgmap.inx2dir(i));
519 566 }
520 567  
521   - // If doing a ring background subtraction, generate an alpha map and
522   - // make sure to cache the counts for each observation
523   - if (m_bkgsubtract == "RING") {
524   - m_alphamap = create_map(m_obs);
525   - m_onmap = create_map(m_obs);
  568 + } // endif: background subtraction selected
  569 +
  570 + // Return
  571 + return;
  572 +}
  573 +
  574 +
  575 +/***********************************************************************//**
  576 + * @brief Generates map of pixel exclusions
  577 + *
  578 + * @param[in] filename Exclusion file name.
  579 + *
  580 + * Generates a sky map of the pixels that are to be excluded from the
  581 + * background estimation. Pixels with values different from 0 will be
  582 + * excluded.
  583 + ***************************************************************************/
  584 +void ctskymap::map_exclusions(const GFilename& filename)
  585 +{
  586 + // Create exlusion map
  587 + m_exclmap = create_map(m_obs);
  588 +
  589 + // Set all pixels to 0 (no pixel excluded)
  590 + m_exclmap = 0.0;
  591 +
  592 + // Make sure the exclusions filename is valid
  593 + if (is_valid_filename(filename)) {
  594 +
  595 + // Fill the exclusions based on the regions supplied
  596 + if (filename.is_fits()) {
  597 + map_exclusions_fits(filename);
  598 + }
  599 + else {
  600 + map_exclusions_region(filename);
526 601 }
527 602  
528   - } // endif: background subtraction selected
  603 + } // endif: filename was valid
  604 +
  605 + // Return
  606 + return;
  607 +}
529 608  
530   - // Write parameters into logger
531   - log_parameters(TERSE);
  609 +
  610 +/***********************************************************************//**
  611 + * @brief Fills exclusions map from FITS image
  612 + *
  613 + * @param[in] filename FITS image file name.
  614 + *
  615 + * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
  616 + * the exclusion sky map FITS file.
  617 + ***************************************************************************/
  618 +void ctskymap::map_exclusions_fits(const GFilename& filename)
  619 +{
  620 + // Load the fits image
  621 + GSkyMap inmap(filename);
  622 +
  623 + // Loop through the individual pixels in the exclusion map
  624 + for (int i = 0; i < m_exclmap.npix(); ++i) {
  625 +
  626 + // Get the pixel direction
  627 + GSkyDir dir = m_exclmap.pix2dir(i);
  628 +
  629 + // Check this sky position in the fits map
  630 + if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
  631 +
  632 + // Set the pixel to 1
  633 + m_exclmap(i) = 1.0;
  634 +
  635 + } // endif: pixel,region overlap check
  636 +
  637 + } // endfor: looped over exclusion map pixels
  638 +
  639 + // Return
  640 + return;
  641 +}
  642 +
  643 +
  644 +/***********************************************************************//**
  645 + * @brief Fills exclusions map from DS9 region file
  646 + *
  647 + * @param[in] filename DS9 region file name.
  648 + *
  649 + * Sets all exclusion map pixels to 1 that are contained in any of the DS9
  650 + * regions.
  651 + ***************************************************************************/
  652 +void ctskymap::map_exclusions_region(const GFilename& filename)
  653 +{
  654 + // Load the exclusion regions
  655 + GSkyRegions regions(filename);
  656 +
  657 + // Loop through the individual pixels in the exclusion map
  658 + for (int i = 0; i < m_exclmap.npix(); ++i) {
  659 +
  660 + // Get the pixel position
  661 + GSkyDir dir = m_exclmap.pix2dir(i);
  662 +
  663 + // If pixel position overlaps with the regions
  664 + if (regions.contains(dir)) {
  665 + m_exclmap(i) = 1.0;
  666 + }
  667 +
  668 + } // endfor: looped over exclusion map pixels
532 669  
533 670 // Return
534 671 return;
... ... @@ -604,21 +741,26 @@ void ctskymap::map_events(GCTAObservation* obs)
604 741 }
605 742  
606 743 // Fill event in skymap
607   - m_skymap(pixel, 0) += 1.0;
  744 + #pragma omp critical(ctskymap_map_events)
  745 + m_skymap(pixel) += 1.0;
608 746 num_in_map++;
609   -
  747 +
610 748 } // endfor: looped over all events
611 749  
612 750 // Log binning results
613   - log_value(NORMAL, "Events in list", obs->events()->size());
614   - log_value(NORMAL, "Events in map", num_in_map);
615   - log_value(NORMAL, "Events outside RoI", num_outside_roi);
616   - log_value(NORMAL, "Events outside map area", num_outside_map);
617   - log_value(NORMAL, "Events outside energies", num_outside_erange);
618   -
619   - // Write sky map into header
620   - log_header1(EXPLICIT, "Sky map");
621   - log_string(EXPLICIT, m_skymap.print(m_chatter));
  751 + #pragma omp critical(ctskymap_map_events)
  752 + {
  753 + log_header3(TERSE, get_obs_header(obs));
  754 + log_value(NORMAL, "Events in list", obs->events()->size());
  755 + log_value(NORMAL, "Events in map", num_in_map);
  756 + log_value(NORMAL, "Events outside RoI", num_outside_roi);
  757 + log_value(NORMAL, "Events outside map area", num_outside_map);
  758 + log_value(NORMAL, "Events outside energies", num_outside_erange);
  759 +
  760 + // Write sky map into header
  761 + log_header1(EXPLICIT, "Sky map");
  762 + log_string(EXPLICIT, m_skymap.print(m_chatter));
  763 + }
622 764  
623 765 // Return
624 766 return;
... ... @@ -646,7 +788,7 @@ void ctskymap::map_background(GCTAObservation* obs)
646 788 map_background_irf(obs);
647 789 }
648 790 else if (m_bkgsubtract == "RING") {
649   - map_background_ring(obs);
  791 + map_background_irf(obs);
650 792 }
651 793  
652 794 // Return
... ... @@ -670,7 +812,8 @@ void ctskymap::map_background(GCTAObservation* obs)
670 812 void ctskymap::map_background_irf(GCTAObservation* obs)
671 813 {
672 814 // Get IRF response
673   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
  815 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
  816 + (obs->response());
674 817  
675 818 // Throw an exception if observation has no instrument response function
676 819 if (rsp == NULL) {
... ... @@ -697,8 +840,16 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
697 840 GEnergy emin(m_emin, "TeV");
698 841 GEnergy emax(m_emax, "TeV");
699 842  
700   - // Extract region of interest from observation
701   - GCTARoi roi = obs->roi();
  843 + // Extract region of interest center and radius from observation and
  844 + // compute the cosine of the RoI radius. This allows to mimimize as
  845 + // much as possible the trigonometric computations. Note that the cosine
  846 + // of an angle is maximal for an angle of zero. This explains later why
  847 + // we chose ">=" to test whether an angle is smaller than a given radius.
  848 + GSkyDir roi_centre = obs->roi().centre().dir();
  849 + double cos_roi_radius = std::cos(obs->roi().radius() * gammalib::deg2rad);
  850 +
  851 + // Extract exposure
  852 + double exposure = obs->livetime();
702 853  
703 854 // Initialise statistics
704 855 double total = 0.0;
... ... @@ -707,25 +858,25 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
707 858 for (int i = 0; i < m_bkgmap.npix(); ++i) {
708 859  
709 860 // Get sky direction of pixel
710   - GSkyDir skydir = m_bkgmap.inx2dir(i);
711   -
712   - // Convert sky direction in instrument direction
713   - GCTAInstDir instdir = obs->pointing().instdir(skydir);
  861 + GSkyDir& skydir = m_dirs[i];
714 862  
715   - // If RoI is valid and instrument direction is not within RoI then
716   - // skip pixel
717   - if (roi.is_valid() && !roi.contains(instdir)) {
  863 + // Skip pixel if it is not within the RoI
  864 + if (roi_centre.cos_dist(skydir) < cos_roi_radius) {
718 865 continue;
719 866 }
720 867  
  868 + // Convert sky direction to instrument direction
  869 + GCTAInstDir instdir = obs->pointing().instdir(skydir);
  870 +
721 871 // Compute background value
722 872 double value = bkg->rate_ebin(instdir, emin, emax);
723 873  
724 874 // Multiply background rate with livetime and solid angle
725   - value *= obs->livetime() * m_solidangle[i];
  875 + value *= exposure * m_solidangle[i];
726 876  
727 877 // Add number of background events to background map
728   - m_bkgmap(i,0) += value;
  878 + #pragma omp critical(ctskymap_map_background_irf)
  879 + m_bkgmap(i) += value;
729 880  
730 881 // Update total number of background events
731 882 total += value;
... ... @@ -733,7 +884,10 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
733 884 } // endfor: looped over background map pixels
734 885  
735 886 // Log background subtraction results
736   - log_value(NORMAL, "Events in background", int(total+0.5));
  887 + #pragma omp critical(ctskymap_map_events)
  888 + {
  889 + log_value(NORMAL, "Events in background", int(total+0.5));
  890 + }
737 891  
738 892 // Return
739 893 return;
... ... @@ -741,196 +895,122 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
741 895  
742 896  
743 897 /***********************************************************************//**
744   - * @brief Estimates the background in sky map based on the ring background
745   - * method
746   - *
747   - * @param[in] obs CTA observation.
748   - *
749   - * @exception GException::invalid_value
750   - * No response information available for observation.
751   - * No background template available in instrument response function.
  898 + * @brief Computes the bin-by-bin significance if background subtraction
  899 + * requested
752 900 *
753   - * Estimates the background in the sky map by summing the events within a
754   - * ring centered at a given pixel's position. The pixels in the ring are
755   - * weighted also by the background IRF value.
  901 + * Method computes the bin-by bin significance if a background subtraction
  902 + * method is specified. If method is "IRF", Poisson statistics (in the
  903 + * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance
  904 + * is computed for each bin.
756 905 ***************************************************************************/
757   -void ctskymap::map_background_ring(GCTAObservation* obs)
  906 +void ctskymap::map_significance(void)
758 907 {
759   - // Get IRF response (to scale background counts)
760   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
761   -
762   - // Throw an exception if observation has no instrument response function
763   - if (rsp == NULL) {
764   - std::string msg = "No response information available for "+
765   - get_obs_header(obs)+" to compute IRF background. "
766   - "Please specify response information or use "
767   - "another background subtraction method.";
768   - throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
  908 + // Compute significance from "RING" method (Li & Ma eq. 17)
  909 + if (m_bkgsubtract == "RING") {
  910 + map_significance_ring();
769 911 }
770 912  
771   - // Get IRF background template
772   - const GCTABackground* bkg = rsp->background();
773   -
774   - // Throw an exception if observation has no IRF background template
775   - if (bkg == NULL) {
776   - std::string msg = "No IRF background template found in instrument "
777   - "response function for "+
778   - get_obs_header(obs)+". Please specify an instrument "
779   - "response function containing a background template.";
780   - throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
  913 + // Compute significance using Poisson statistics in the Gaussian limit
  914 + else {
  915 + m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap);
781 916 }
782 917  
783   - // Set minimum and maximum energy as GEnergy instances
784   - GEnergy emin(m_emin, "TeV");
785   - GEnergy emax(m_emax, "TeV");
786   -
787   - // Extract region of interest from observation
788   - GCTARoi roi = obs->roi();
789   -
790   - // Extract exposure
791   - double exposure = obs->livetime();
792   -
793   - // Loop over all map pixels
794   - for (int i = 0; i < m_bkgmap.npix(); ++i) {
795   -
796   - // Get sky direction of pixel
797   - GSkyDir& skydir = m_dirs[i];
798   -
799   - // Convert sky direction in instrument direction
800   - GCTAInstDir instdir = obs->pointing().instdir(skydir);
801   -
802   - // If RoI is valid and instrument direction is not within RoI then
803   - // skip pixel
804   - if (roi.is_valid() && !roi.contains(instdir)) {
805   - continue;
806   - }
807   -
808   - // Compute background value
809   - m_alphamap(i,0) += bkg->rate_ebin(instdir, emin, emax) *
810   - m_solidangle[i] * exposure;
811   -
812   - // Store the counts in this bin
813   - m_onmap(i,0) += m_skymap(i,0);
814   -
815   - } // endfor: Loop for caching bkg IRF sensitivity
816   -
817   - // Zero out the counts map to prepare it for the next observation
818   - m_skymap = 0.0;
819   -
820 918 // Return
821 919 return;
822 920 }
823 921  
824 922  
825 923 /***********************************************************************//**
826   - * @brief Write keywords in FITS HDU
827   - *
828   - * @param[in,out] hdu Pointer to FITS HDU.
  924 + * @brief Computes the bin-by-bin significance for RING background
829 925 *
830   - * Writes keywords in FITS HDU.
  926 + * Computes the a Li & Ma significance for each sky map bin and replaces
  927 + * the sky and background maps by the On- and Off-count maps.
831 928 ***************************************************************************/
832   -void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
  929 +void ctskymap::map_significance_ring(void)
833 930 {
834   - // Continue only if pointer is valid
835   - if (hdu != NULL) {
  931 + // Log message about what is being done
  932 + log_header1(NORMAL, "Computing Ring background map");
  933 + log_value(NORMAL, "Total pixels to process", m_skymap.npix());
836 934  
837   - // Set keywords
838   - hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
839   - hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary");
840   - hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary");
841   - hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX");
842   -
843   - } // endif: pointer was valid
  935 + // Initialise On- and Off-count maps
  936 + GSkyMap onmap = m_skymap;
  937 + GSkyMap offmap = m_skymap;
  938 + onmap = 0.0;
  939 + offmap = 0.0;
844 940  
845   - // Return
846   - return;
847   -}
  941 + // Initialise statistics
  942 + int num_bad_alpha = 0;
  943 +
  944 + // Set progress logging frequency. By default the user is updated on the
  945 + // progress when another 10% of pixels is processed, however, if there
  946 + // are many pixels the user is updated each 100000 pixels.
  947 + int n_logpix = m_skymap.npix() / 10;
  948 + if (n_logpix > 100000) {
  949 + n_logpix = 100000;
  950 + }
848 951  
  952 + // Loop over sky map pixels
  953 + #pragma omp parallel for
  954 + for (int i = 0; i < m_skymap.npix(); ++i) {
849 955  
850   -/***********************************************************************//**
851   - * @brief Computes the bin-by-bin significance if background subtraction
852   - * requested
853   - *
854   - * Method computes the bin-by bin significance if a background subtraction
855   - * method is specified. If method is "IRF", Poisson statistics (in the
856   - * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance
857   - * is computed for each bin (as well as a map of the computed alpha values).
858   - ***************************************************************************/
859   -void ctskymap::map_significance(void)
860   -{
861   - // Compute significance from "RING" method (Li & Ma eq. 17)
862   - if (m_bkgsubtract == "RING") {
863   -
864   - // Log message about what is being done
865   - log_header1(NORMAL, "Computing Ring background map");
866   - log_value(NORMAL, "Total pixels to process", m_onmap.npix());
  956 + // Initialise the on/off counts and alpha for this bin
  957 + double n_on = 0.0;
  958 + double n_off = 0.0;
  959 + double alpha = 0.0;
867 960  
868   - // Store the number of bins with an inappropriate alpha parameter
869   - int num_bad_alpha = 0;
  961 + // Update user about progress
  962 + if (i % n_logpix == 0) {
  963 + #pragma omp critical(ctskymap_map_significance_ring)
  964 + log_value(NORMAL, "Pixels remaining", m_skymap.npix()-i);
  965 + }
870 966  
871   - // Loop through each bin in the on-counts map
872   - for (int i = 0; i < m_onmap.npix(); ++i) {
  967 + // Get bin coordinates
  968 + GSkyDir& skydir = m_dirs[i];
873 969  
874   - // Initialise the on/off counts and alpha for this bin
875   - double n_on = 0.0;
876   - double n_off = 0.0;
877   - double alpha = 0.0;
  970 + // Compute the alpha and counts for this bin
  971 + compute_ring_values(m_skymap, m_bkgmap, skydir, n_on, n_off, alpha);
878 972  
879   - // Since this can take a long time, keep the user updated on
880   - // the progress when another 10% of pixels is processed
881   - if (i % (m_onmap.npix()/10) == 0) {
882   - log_value(NORMAL, "Pixels remaining", m_onmap.npix()-i);
883   - }
  973 + // Store the On and alpha-weighted Off-counts
  974 + onmap(i) = n_on;
  975 + offmap(i) = alpha * n_off;
884 976  
885   - // Get bin coordinates
886   - GSkyDir& skydir = m_dirs[i];
  977 + // If alpha is zero then increment the bad alpha counter
  978 + if (alpha == 0.0) {
  979 + num_bad_alpha++;
  980 + }
887 981  
888   - // Compute the alpha and counts for this bin
889   - compute_ring_values(m_onmap, m_alphamap, skydir,
890   - n_on, n_off, alpha);
  982 + // ... otherwise store the results
  983 + else {
891 984  
892   - // If alpha is zero then increment the bad alpha counter
893   - if (alpha == 0.0) {
894   - num_bad_alpha++;
  985 + // Compute and store significance (Li & Ma eq. 17)
  986 + double sigma;
  987 + if (n_on == 0.0) {
  988 + sigma = -2.0 * n_off * std::log(1.0+alpha);
895 989 }
896   -
897   - // ... otherwise store the results
898 990 else {
  991 + sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
  992 + (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
  993 + n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
  994 + }
  995 + #pragma omp critical(ctskymap_map_significance_ring)
  996 + m_sigmap(i) = sigma;
899 997  
900   - // Store the on-counts and alpha-weighted off-counts
901   - m_skymap(i,0) = n_on;
902   - m_bkgmap(i,0) = alpha * n_off;
903   -
904   - // Compute and store significance (Li & Ma eq. 17)
905   - if (n_on == 0.0) {
906   - m_sigmap(i,0) = -2.0 * n_off * std::log(1.0+alpha);
907   - }
908   - else {
909   - m_sigmap(i,0) = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
910   - (n_on * std::log((1.0+alpha) *
911   - n_on / (alpha * (n_on+n_off))) +
912   - n_off * std::log((1.0+alpha) *
913   - n_off / (n_on+n_off)));
914   - }
915   -
916   - } // endelse: alpha was non-zero
917   -
918   - } // endfor: looped over all pixels
  998 + } // endelse: alpha was non-zero
919 999  
920   - // Log the number of bad-alpha bins
921   - log_value(NORMAL, "Bins with alpha=0", num_bad_alpha);
  1000 + } // endfor: looped over all pixels
922 1001  
923   - // Now take the square root. Since some bins can be negative, first
924   - // take the absolute value of the map, square root that, then multiply
925   - // each bin by its sign to preserve the +/- significance.
926   - m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
  1002 + // Log the number of bad-alpha bins
  1003 + log_value(NORMAL, "Bins with alpha=0", num_bad_alpha);
927 1004  
928   - } // endif: ring background selected
  1005 + // Now take the square root. Since some bins can be negative, first
  1006 + // take the absolute value of the map, square root that, then multiply
  1007 + // each bin by its sign to preserve the +/- significance.
  1008 + m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
929 1009  
930   - // Compute significance using Poisson statistics in the Gaussian limit
931   - else {
932   - m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap);
933   - }
  1010 + // Terminate by storing the on-counts map as sky map and the off-counts
  1011 + // map as background map
  1012 + m_skymap = onmap;
  1013 + m_bkgmap = offmap;
934 1014  
935 1015 // Return
936 1016 return;
... ... @@ -940,18 +1020,18 @@ void ctskymap::map_significance(void)
940 1020 /***********************************************************************//**
941 1021 * @brief Computes Non, Noff and alpha for a counts map and sensitivity map
942 1022 *
943   - * @param[in] counts Counts map
944   - * @param[in] sensitivity Sensitivity map
945   - * @param[in] position Position of the ring & roi centers
946   - * @param[out] non Returned estimate of ON counts
947   - * @param[out] noff Returned estimate of OFF counts
948   - * @param[out] alpha Returned estimate of alpha
  1023 + * @param[in] counts Counts map
  1024 + * @param[in] background Background map
  1025 + * @param[in] position Position of the ring & roi centers
  1026 + * @param[out] non Returned estimate of ON counts
  1027 + * @param[out] noff Returned estimate of OFF counts
  1028 + * @param[out] alpha Returned estimate of alpha
949 1029 *
950 1030 * This method computes the Non counts Noff values at a given position. It
951 1031 * also computes the alpha parameter from the passed sensitivity map.
952 1032 ***************************************************************************/
953 1033 void ctskymap::compute_ring_values(const GSkyMap& counts,
954   - const GSkyMap& sensitivity,
  1034 + const GSkyMap& background,
955 1035 const GSkyDir& position,
956 1036 double& non,
957 1037 double& noff,
... ... @@ -984,13 +1064,13 @@ void ctskymap::compute_ring_values(const GSkyMap&amp; counts,
984 1064 if (position.cos_dist(skydir) >= cos_outradius) {
985 1065  
986 1066 // Check if pixel is inside the background region
987   - if ((m_exclmap(j,0) == 0.0) && position.cos_dist(skydir) < cos_inradius) {
  1067 + if ((m_exclmap(j) == 0.0) && position.cos_dist(skydir) < cos_inradius) {
988 1068  
989 1069 // Update n_off
990   - noff += counts(j,0);
  1070 + noff += counts(j);
991 1071  
992 1072 // Update alpha_off
993   - alpha_off += sensitivity(j,0);
  1073 + alpha_off += background(j);
994 1074  
995 1075 }
996 1076  
... ... @@ -998,10 +1078,10 @@ void ctskymap::compute_ring_values(const GSkyMap&amp; counts,
998 1078 else if (position.cos_dist(skydir) >= cos_roiradius) {
999 1079  
1000 1080 // Update n_on for significance computation
1001   - non += counts(j,0);
  1081 + non += counts(j);
1002 1082  
1003 1083 // Update alpha_on
1004   - alpha_on += sensitivity(j,0);
  1084 + alpha_on += background(j);
1005 1085  
1006 1086 } // endif: source and background region check
1007 1087  
... ... @@ -1026,121 +1106,32 @@ void ctskymap::compute_ring_values(const GSkyMap&amp; counts,
1026 1106  
1027 1107  
1028 1108 /***********************************************************************//**
1029   - * @brief Generates map of pixel exclusions
1030   - *
1031   - * @param[in] filename Exclusion file name.
1032   - *
1033   - * Generates a sky map of the pixels that are to be excluded from the
1034   - * background estimation. Pixels with values different from 0 will be
1035   - * excluded.
1036   - ***************************************************************************/
1037   -void ctskymap::map_exclusions(const GFilename& filename)
1038   -{
1039   - // Create exlusion map
1040   - m_exclmap = create_map(m_obs);
1041   -
1042   - // Set all pixels to 0 (no pixel excluded)
1043   - m_exclmap = 0.0;
1044   -
1045   - // Make sure the exclusions filename is valid
1046   - if (is_valid_filename(filename)) {
1047   -
1048   - // Fill the exclusions based on the regions supplied
1049   - if (filename.is_fits()) {
1050   - map_exclusions_fits(filename);
1051   - }
1052   - else {
1053   - map_exclusions_reg(filename);
1054   - }
1055   -
1056   - } // endif: filename was valid
1057   -
1058   - // Return
1059   - return;
1060   -}
1061   -
1062   -
1063   -/***********************************************************************//**
1064   - * @brief Fills exclusions map from FITS image
1065   - *
1066   - * @param[in] filename FITS image file name.
1067   - *
1068   - * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
1069   - * the exclusion sky map FITS file.
1070   - ***************************************************************************/
1071   -void ctskymap::map_exclusions_fits(const GFilename& filename)
1072   -{
1073   - // Load the fits image
1074   - GSkyMap inmap(filename);
1075   -
1076   - // Loop through the individual pixels in the exclusion map
1077   - for (int i = 0; i < m_exclmap.npix(); ++i) {
1078   -
1079   - // Get the pixel direction
1080   - GSkyDir dir = m_exclmap.pix2dir(i);
1081   -
1082   - // Check this sky position in the fits map
1083   - if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
1084   -
1085   - // Set the pixel to 1
1086   - m_exclmap(i) = 1.0;
1087   -
1088   - } // endif: pixel,region overlap check
1089   -
1090   - } // endfor: looped over exclusion map pixels
1091   -
1092   - // Return
1093   - return;
1094   -}
1095   -
1096   -
1097   -/***********************************************************************//**
1098   - * @brief Fills exclusions map from DS9 region file
  1109 + * @brief Write keywords in FITS HDU
1099 1110 *
1100   - * @param[in] filename DS9 region file name.
  1111 + * @param[in,out] hdu Pointer to FITS HDU.
1101 1112 *
1102   - * Sets all exclusion map pixels to 1 that are contained in any of the DS9
1103   - * regions.
  1113 + * Writes keywords in FITS HDU.
1104 1114 ***************************************************************************/
1105   -void ctskymap::map_exclusions_reg(const GFilename& filename)
  1115 +void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
1106 1116 {
1107   - // Load the exclusion regions
1108   - GSkyRegions regions(filename);
1109   -
1110   - // Loop through the individual pixels in the exclusion map
1111   - for (int i = 0; i < m_exclmap.npix(); ++i) {
  1117 + // Continue only if pointer is valid
  1118 + if (hdu != NULL) {
1112 1119  
1113   - // Get the pixel position
1114   - GSkyDir dir = m_exclmap.pix2dir(i);
  1120 + // Set keywords
  1121 + hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
  1122 + hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary");
  1123 + hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary");
  1124 + hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX");
1115 1125  
1116   - // If pixel position overlaps with the regions
1117   - if (regions.contains(dir)) {
1118   - m_exclmap(i) = 1.0;
  1126 + // Write RING background keywords
  1127 + if (m_bkgsubtract == "RING") {
  1128 + hdu->card("ROIRAD", m_roiradius, "[deg] Source region radius");
  1129 + hdu->card("INRAD", m_inradius, "[deg] Inner background ring radius");
  1130 + hdu->card("OUTRAD", m_outradius, "[deg] Outer background ring radius");
1119 1131 }
1120 1132  
1121   - } // endfor: looped over exclusion map pixels
  1133 + } // endif: pointer was valid
1122 1134  
1123 1135 // Return
1124 1136 return;
1125 1137 }
1126   -
1127   -
1128   -/***********************************************************************//**
1129   - * @brief Estimates the background in sky map based on IRF template
1130   - *
1131   - * @param[in] lnE Natural logarithm of energy in MeV.
1132   - ***************************************************************************/
1133   -double ctskymap::irf_kern::eval(const double& lnE)
1134   -{
1135   - // Get log10 of energy in TeV
1136   - double logE = lnE * gammalib::inv_ln10 - 6;
1137   -
1138   - // Get function value
1139   - double value = (*m_bgd)(logE, m_dir->detx(), m_dir->dety());
1140   -
1141   - // Correct for variable substitution
1142   - value *= std::exp(lnE);
1143   -
1144   - // Return value
1145   - return value;
1146   -}
... ...
src/ctskymap/ctskymap.hpp
  View file @217ed01
... ... @@ -68,35 +68,24 @@ protected:
68 68 void copy_members(const ctskymap& app);
69 69 void free_members(void);
70 70 void get_parameters(void);
  71 + void setup_maps(void);
  72 + void map_exclusions(const GFilename& filename);
  73 + void map_exclusions_fits(const GFilename& filename);
  74 + void map_exclusions_region(const GFilename& filename);
71 75 void map_events(GCTAObservation* obs);
72 76 void map_background(GCTAObservation* obs);
73 77 void map_background_irf(GCTAObservation* obs);
74   - void map_background_ring(GCTAObservation* obs);
75 78 void map_significance(void);
76   - void map_exclusions(const GFilename& filename);
77   - void map_exclusions_fits(const GFilename& filename);
78   - void map_exclusions_reg(const GFilename& filename);
79   - void compute_ring_values(const GSkyMap& counts,
80   - const GSkyMap& sensitivity,
  79 + void map_significance_ring(void);
  80 + void compute_ring_values(const GSkyMap& counts,
  81 + const GSkyMap& background,
81 82 const GSkyDir& position,
82 83 double& non, double& noff, double& alpha);
83 84 void write_hdu_keywords(GFitsHDU* hdu) const;
84 85  
85   - // Background integration kernel
86   - class irf_kern : public GFunction {
87   - public:
88   - irf_kern(const GCTABackground* bgd,
89   - const GCTAInstDir* dir) :
90   - m_bgd(bgd),
91   - m_dir(dir) { }
92   - double eval(const double& lnE);
93   - protected:
94   - const GCTABackground* m_bgd; //!< Pointer to background
95   - const GCTAInstDir* m_dir; //!< Pointer to instrument direction
96   - };
97   -
98 86 // User parameters
99 87 GFilename m_outmap; //!< Output file name
  88 + GFilename m_inexclusion; //!< Exclusion map file name
100 89 double m_emin; //!< Minimum energy (TeV)
101 90 double m_emax; //!< Maximum energy (TeV)
102 91 std::string m_bkgsubtract; //!< Background subtraction method
... ... @@ -111,8 +100,6 @@ protected:
111 100 GSkyMap m_bkgmap; //!< Background map
112 101 GSkyMap m_sigmap; //!< Significance map
113 102 GSkyMap m_exclmap; //!< Exclusion map for RING background
114   - GSkyMap m_alphamap; //!< Alpha map values for RING background
115   - GSkyMap m_onmap; //!< On counts map for RING background
116 103  
117 104 // Caching variables to prevent multiple computations of the same thing
118 105 std::vector<double> m_solidangle; //!< Cached pixel solid angles
... ...
src/ctskymap/ctskymap.par
  View file @217ed01
... ... @@ -44,8 +44,8 @@ yref, r, a, 22.01,-90,90, &quot;Second coordinate of image center in degrees (DEC
44 44 #==================================
45 45 bkgsubtract, s, a, NONE, NONE|IRF|RING,, "Background subtraction method"
46 46 roiradius, r, a, 0.2,,, "Source region radius for estimating on-counts (degrees)"
47   -inradius, r, a, 0.6,,, "Inner ring radius (degrees)"
48   -outradius, r, a, 0.8,,, "Outer ring radius (degrees)"
  47 +inradius, r, a, 0.6,,, "Inner background ring radius (degrees)"
  48 +outradius, r, a, 0.8,,, "Outer background ring radius (degrees)"
49 49 inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name"
50 50  
51 51 #
... ...