diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst index 077332c..aed4d0f 100644 --- a/doc/source/users/reference_manual/ctskymap.rst +++ b/doc/source/users/reference_manual/ctskymap.rst @@ -74,7 +74,7 @@ General parameters ``yref [real]`` Declination / Galactic latitude of image centre (J2000, in degrees). -``bkgsubtract [string]`` +``bkgsubtract [string]`` Background subtraction method. diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 97883e3..201fa99 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -390,6 +390,14 @@ void ctskymap::copy_members(const ctskymap& app) m_bkgsubtract = app.m_bkgsubtract; m_publish = app.m_publish; m_chatter = app.m_chatter; + m_solidangle = app.m_solidangle; + + // Variables for RING bkgd subtraction + m_exclmap = app.m_exclmap; + m_roiradius = app.m_roiradius; + m_inradius = app.m_inradius; + m_outradius = app.m_outradius; + m_cnts = app.m_cnts; // Return return; @@ -442,9 +450,22 @@ void ctskymap::get_parameters(void) // Read ahead parameters if (read_ahead()) { m_outmap = (*this)["outmap"].filename(); - m_roiradius = (*this)["roiradius"].real(); - m_inradius = (*this)["inradius"].real(); - m_outradius = (*this)["outradius"].real(); + + if (m_bkgsubtract == "RING") { + // Get the integration region radii + m_roiradius = (*this)["roiradius"].real(); + m_inradius = (*this)["inradius"].real(); + m_outradius = (*this)["outradius"].real(); + + // Make sure the values are reasonable + if (m_roiradius > m_inradius) { + std::string err("'roiradius' must be smaller than inner ring radius"); + throw GException::invalid_value(err); + } else if (m_inradius > m_outradius) { + std::string err("Inner ring radius must be smaller than outer ring radius"); + throw GException::invalid_value(err); + } // endif: Ring radius checks + } } // Create background map and significance map if background subtraction @@ -749,11 +770,8 @@ void ctskymap::map_background_ring(GCTAObservation* obs) double exposure = obs->livetime() ; // Cached background IRF sensitivities - std::vector sens(m_cnts.size(), 0.0); - - // Keep track of the minimum and maximum pixels that need to be looped over - int imin(m_bkgmap.npix()); - int imax(0); + std::vector sens(m_cnts.size(), 0.0); + std::vector obsoverlap(m_cnts.size(), false); // Cache the IRF sensitivities // Loop over all map pixels @@ -771,15 +789,14 @@ void ctskymap::map_background_ring(GCTAObservation* obs) continue; } + // Notify that this bin overlaps the observation + obsoverlap[i] = true; + // If bin is excluded in the exclusion, skip it if (m_exclmap(i) == 0.0) { continue; } - // Update pixel ranges - if (i > imax) imax = i; - if (i < imin) imin = i; - // Setup integration function ctskymap::irf_kern integrand(bkg, &instdir); GIntegral integral(&integrand); @@ -802,18 +819,14 @@ void ctskymap::map_background_ring(GCTAObservation* obs) GSkyRegionCircle inner_reg(0.0, 0.0, m_inradius); GSkyRegionCircle outer_reg(0.0, 0.0, m_outradius); +log_string(NORMAL, "LOOPING ON BKG IRF DONE"); +log_string(NORMAL, "LOOPING ON EVENTS!"); // Loop over all map pixels - for (int i = imin; i <= imax; ++i) { - - // Get sky direction of pixel - GSkyDir skydir = m_bkgmap.inx2dir(i); + for (int i = 0; i < m_bkgmap.npix(); ++i) { - // Convert sky direction in instrument direction - GCTAInstDir instdir = obs->pointing().instdir(skydir); + log_value(NORMAL, "PIXEL #", i); - // If RoI is valid and instrument direction is not within RoI then - // skip pixel - if (roi.is_valid() && !roi.contains(instdir)) { + if (!obsoverlap[i]) { continue; } @@ -825,20 +838,16 @@ void ctskymap::map_background_ring(GCTAObservation* obs) // Loop again over every pixel in the observation to compute Non/Noff int n_on(0), n_off(0); double alpha_on(0.0), alpha_off(0.0); - for (int j = imin; j < imax; ++j) { - - // Get sky direction of pixel - skydir = m_bkgmap.inx2dir(i); + for (int j = 0; j < m_bkgmap.npix(); ++j) { - // Convert sky direction in instrument direction - instdir = obs->pointing().instdir(skydir); - - // If RoI is valid and instrument direction is not within RoI then - // skip pixel - if (roi.is_valid() && !roi.contains(instdir)) { + // Check if the observation overlaps this pixel + if (!obsoverlap[j]) { continue; } + // Get sky direction of pixel + skydir = m_bkgmap.inx2dir(i); + // Check if pixel is inside source region if (roi_reg.contains(skydir)) { // Update n_on @@ -870,9 +879,6 @@ void ctskymap::map_background_ring(GCTAObservation* obs) } // endfor: looped over background map pixels - // Update the likelihood values in the - - // Zero out the alpha map to prepare it for the next observation std::vector::iterator iter; for (iter=m_cnts.begin(); iter!=m_cnts.end(); ++iter) { diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index 381a733..c446cfb 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -43,9 +43,9 @@ yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC # Background subtraction parameters #================================== bkgsubtract, s, a, NONE, NONE|IRF|RING,, "Background subtraction method" -roiradius, r, h, 0.0,,, "Region of interest radius" -inradius, r, h, 0.6,,, "Inner ring radius" -outradius, r, h, 0.8,,, "Outer ring radius" +roiradius, r, a, 0.0,,, "Region of interest radius" +inradius, r, a, 0.6,,, "Inner ring radius" +outradius, r, a, 0.8,,, "Outer ring radius" regfile, f, h, NONE,,, "Exclusion regions definition file" #