diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index dbba584..270d8f5 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -232,7 +232,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.98, \ irf='South_50h', caldb='prod2', obsid='000000', mjdref=51544.5): """ Set a single CTA observation - + The function sets a single CTA observation containing an empty CTA event list. By looping over this function CTA observations can be added to the observation container. @@ -409,7 +409,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.98, \ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): """ Sets a number of standard patterns - + Parameters ---------- pattern : str @@ -466,7 +466,7 @@ def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb, deadc=0.98, pattern='single', offset=1.5): """ Set an observation container filled with CTA observations - + Parameters ---------- ra : float @@ -577,8 +577,17 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, else: usepnt = False + # Set energy limits, with larger etrue limits for the case that energy + # dispersion is requested + if edisp: + _emin = 0.5 * emin + _emax = 1.5 * emax + else: + _emin = emin + _emax = emax + # Set number of energy bins to at least 30 per energy decade - _enumbins = int((math.log10(emax) - math.log10(emin)) * 30.0) + _enumbins = int((math.log10(_emax) - math.log10(_emin)) * 30.0) if enumbins > _enumbins: _enumbins = enumbins @@ -598,8 +607,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, expcube['nxpix'] = nxpix expcube['nypix'] = nypix expcube['enumbins'] = _enumbins - expcube['emin'] = emin - expcube['emax'] = emax + expcube['emin'] = _emin + expcube['emax'] = _emax expcube['coordsys'] = coordsys expcube['proj'] = proj expcube['addbounds'] = addbounds @@ -621,8 +630,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, psfcube['nxpix'] = psf_nxpix psfcube['nypix'] = psf_nypix psfcube['enumbins'] = _enumbins - psfcube['emin'] = emin - psfcube['emax'] = emax + psfcube['emin'] = _emin + psfcube['emax'] = _emax psfcube['coordsys'] = coordsys psfcube['proj'] = proj psfcube['addbounds'] = addbounds @@ -673,8 +682,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, edispcube['nxpix'] = psf_nxpix edispcube['nypix'] = psf_nypix edispcube['enumbins'] = _enumbins - edispcube['emin'] = emin - edispcube['emax'] = emax + edispcube['emin'] = _emin + edispcube['emax'] = _emax edispcube['coordsys'] = coordsys edispcube['proj'] = proj edispcube['addbounds'] = addbounds @@ -723,7 +732,7 @@ def get_stacked_obs(cls, obs): # Write header if cls._logExplicit(): cls._log.header3('Binning events') - + # Bin events cntcube = ctools.ctbin(obs) cntcube['usepnt'] = False diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 27c20a0..30a8cfc 100755 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -247,12 +247,21 @@ void ctbin::run(void) // Find all unbinned CTA observations in m_obs std::vector obs_list(0); + double mean_ra = 0.0; + double mean_dec = 0.0; for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; obs = next_unbinned_observation()) { // Push observation into list obs_list.push_back(obs); + // Get pointing + const GCTAPointing& pnt = obs->pointing(); + + // Add to coordinates + mean_ra += pnt.dir().ra_deg(); + mean_dec += pnt.dir().dec_deg(); + // Write message std::string msg = " Including unbinned "+obs->instrument()+ " observation"; @@ -263,6 +272,13 @@ void ctbin::run(void) // Set number of relevant observations int nobs = obs_list.size(); + // Compute mean pointing + if (nobs > 0) { + mean_ra /= double(nobs); + mean_dec /= double(nobs); + m_mean_pnt.radec_deg(mean_ra, mean_dec); + } + // Initialise arrays if (nobs > 0) { m_counts = std::vector(nobs); @@ -275,6 +291,13 @@ void ctbin::run(void) m_cubes.clear(); } + // Log observation selection results + std::string use = (*this)["usepnt"].boolean() ? "yes" : "no"; + log_header3(NORMAL, "Summary"); + log_value(NORMAL, "Number of observations", nobs); + log_value(NORMAL, "Mean pointing", m_mean_pnt.print()); + log_value(NORMAL, "Use mean pointing", use); + // Write header into logger log_header1(TERSE, gammalib::number("Bin observation", nobs)); @@ -584,7 +607,7 @@ void ctbin::publish(const std::string& name) void ctbin::init_members(void) { // Initialise members - m_usepnt = false; + //m_usepnt = false; m_stack = true; m_prefix.clear(); m_publish = false; @@ -594,6 +617,7 @@ void ctbin::init_members(void) m_cubes.clear(); m_counts.clear(); m_weights.clear(); + m_mean_pnt.clear(); m_ebounds.clear(); m_gti.clear(); m_ontime = 0.0; @@ -615,7 +639,7 @@ void ctbin::init_members(void) void ctbin::copy_members(const ctbin& app) { // Copy attributes - m_usepnt = app.m_usepnt; + //m_usepnt = app.m_usepnt; m_stack = app.m_stack; m_prefix = app.m_prefix; m_publish = app.m_publish; @@ -625,6 +649,7 @@ void ctbin::copy_members(const ctbin& app) m_cubes = app.m_cubes; m_counts = app.m_counts; m_weights = app.m_weights; + m_mean_pnt = app.m_mean_pnt; m_ebounds = app.m_ebounds; m_gti = app.m_gti; m_ontime = app.m_ontime; @@ -746,7 +771,17 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs) double xref = 0.0; double yref = 0.0; bool usepnt = (*this)["usepnt"].boolean(); - if (!usepnt) { + if (usepnt) { + if (gammalib::toupper(coordsys) == "GAL") { + xref = m_mean_pnt.l_deg(); + yref = m_mean_pnt.b_deg(); + } + else { + xref = m_mean_pnt.ra_deg(); + yref = m_mean_pnt.dec_deg(); + } + } + else { xref = (*this)["xref"].real(); yref = (*this)["yref"].real(); } @@ -759,29 +794,10 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs) // Read energy boundaries GEbounds ebounds = create_ebounds(); - // If requested, get pointing from observations - if (usepnt) { - - // Get pointing direction for observation - GSkyDir pnt = obs->pointing().dir(); - - // Set xref/yref based on the coordinate system - if (gammalib::toupper(coordsys) == "GAL") { - xref = pnt.l_deg(); - yref = pnt.b_deg(); - } - else { - xref = pnt.ra_deg(); - yref = pnt.dec_deg(); - } - - } // endif: got pointing from observation - // Initialise counts cube GSkyMap cube = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, ebounds.size()); - // Return cube return cube; } diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index d3bf9b6..e60f01a 100755 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -99,7 +99,6 @@ protected: void obs_cube_stacked(void); // User parameters - bool m_usepnt; //!< Use pointing instead of xref/yref parameters bool m_stack; //!< Output one stacked cube or multiple cubes std::string m_prefix; //!< Prefix for multiple cubes bool m_publish; //!< Publish counts cube? @@ -109,6 +108,7 @@ protected: std::vector m_cubes; //!< Event cubes std::vector m_counts; //!< List of event cube counts std::vector m_weights; //!< List of event cube weights + GSkyDir m_mean_pnt; //!< Mean pointing GEbounds m_ebounds; //!< Energy boundaries GGti m_gti; //!< Stacked Good time intervals double m_ontime; //!< Total ontime