Olga Sergijenko / gammalib

Commit 579b020818f9c001062d96e62198743d97599529

Authored by Andreas Specovius
Committed by Jürgen Knödlseder
1 parent 9f76aaf7

Implemented the map production for rotated GSkyRegionMap::set_region_rect() maki… 

…ng now use of the new PA.
Inline Side-by-side
Showing 1 changed file with 43 additions and 18 deletions
src/sky/GSkyRegionMap.cpp
  View file @579b020
... ... @@ -637,17 +637,49 @@ void GSkyRegionMap::set_region_circle(const GSkyRegionCircle* circle)
637 637 * Sets a region map from a region rectangle.
638 638 * A skymap is allocated with a resolution goal of 0.01 degrees (same for both
639 639 * axes).
640   - * The particular resolution of either axes depends on the w/h of the rectangle.
641   - * The sky map will be in celestial coordinates and TAN projection.
  640 + * The particular resolution of either axes depends on the w/h of the RoI of the
  641 + * rotated rectangle.
  642 + * The map extent is chosen dynamically to match the RoI of the rotated rectangle.
  643 + * Three sparse pixels are added at either side of the map.
  644 + * The sky map is created in celestial coordinates and TAN projection.
642 645 ***************************************************************************/
643 646 void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect)
644 647 {
645 648 // Define desired resolution (in degrees)
646 649 const double resolution_goal = 0.01;
647 650  
  651 + // Try to estimate extent of the (rotated) rectangle:
  652 + // Loop over the rectangle corners to get the extension
  653 + double dxmax=0;
  654 + double dymax=0;
  655 + for(int icorner=0; icorner<4; ++icorner) {
  656 +
  657 + // Get skydir of current corner
  658 + GSkyDir corner = rect->get_corner(icorner);
  659 +
  660 + // Get extension along declination
  661 + double dy = std::abs(corner.dec() - rect->centre().dec());
  662 + if (dy > dymax) {
  663 + dymax = dy;
  664 + }
  665 +
  666 + // Get extension along right ascension, correct for spherical coordsys
  667 + double dx = std::abs(corner.ra() - rect->centre().ra()) * std::cos(corner.dec());
  668 + if (dx > dxmax) {
  669 + dxmax = dx;
  670 + }
  671 + }
  672 +
  673 + // Compute width of map RoI in degrees
  674 + double map_width = 2*dxmax *gammalib::rad2deg;
  675 + double map_height = 2*dymax *gammalib::rad2deg;
  676 + // std::cout<<"Generate map for rectangle:\n"<<(*rect)<<std::endl;
  677 + // std::cout<<"Found map extension: (w,h)=("<<map_width<<","<<map_height<<")\n";
  678 +
648 679 // Compute num of bins in x direction, round
649   - int nbins_x = int(rect->width() / resolution_goal +0.5);
650   - int nbins_y = int(rect->height() / resolution_goal +0.5);
  680 + int nbins_x = int(map_width / resolution_goal +0.5);
  681 + int nbins_y = int(map_height / resolution_goal +0.5);
  682 + // std::cout<<" nbins: x="<<nbins_x<<", y="<<nbins_y<<std::endl;
651 683  
652 684 // Ensure a minimum num of bins
653 685 if (nbins_x < 10) {
... ... @@ -658,33 +690,26 @@ void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect)
658 690 }
659 691  
660 692 // Compute actual resolution using nbins
661   - double resolution_x = rect->width() / nbins_x;
662   - double resolution_y = rect->height() / nbins_y;
  693 + double resolution_x = map_width / nbins_x;
  694 + double resolution_y = map_height / nbins_y;
  695 + // std::cout<<" resolution for roi: x="<<resolution_x<<", y="<<resolution_y<<std::endl;
663 696  
664 697 // Add leakage pixels
665 698 nbins_x += 6;
666 699 nbins_y += 6;
667 700  
668   - // Create coarse 9-pixel map to later avoid "bending" of fine box-map
669   - GSkyMap map_coarse("TAN", "CEL", rect->ra(), rect->dec(),
670   - rect->width(), rect->height(),
671   - 3, 3
672   - );
673   - map_coarse(4) = 1.0;
674   - GSkyRegionMap regionmap_coarse = GSkyRegionMap(map_coarse);
675   -
676 701 // Create fine sky map
677 702 GSkyMap map_fine("TAN", "CEL", rect->ra(), rect->dec(),
678   - resolution_x, resolution_y,
  703 + -resolution_x, resolution_y,
679 704 nbins_x, nbins_y
680 705 );
681 706  
682 707 // Set sky map pixels
683   - for (int i = 0; i < map_fine.npix(); ++i) {
  708 + for (int ipix = 0; ipix < map_fine.npix(); ++ipix) {
684 709  
685 710 // If pixel is contained in rectangle then set map pixel to 1
686   - if (regionmap_coarse.contains(map_fine.inx2dir(i))) {
687   - map_fine(i) = 1.0;
  711 + if (rect->contains(map_fine.inx2dir(ipix))) {
  712 + map_fine(ipix) = 1.0;
688 713 }
689 714  
690 715 } // endfor: looped over sky map pixels
... ...