Iftach Sadeh / ctools

Commit 4e9d5a1a7f9f2a8999ed3d7352f6d2ca11c9d562

Authored by Jürgen Knödlseder
1 parent 7534afae

Compute mean pointing direction in ctbin

Inline Side-by-side
Showing 3 changed files with 59 additions and 34 deletions
cscripts/obsutils.py
  View file @4e9d5a1
... ... @@ -232,7 +232,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.98, \
232 232 irf='South_50h', caldb='prod2', obsid='000000', mjdref=51544.5):
233 233 """
234 234 Set a single CTA observation
235   -
  235 +
236 236 The function sets a single CTA observation containing an empty CTA
237 237 event list. By looping over this function CTA observations can be
238 238 added to the observation container.
... ... @@ -409,7 +409,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.98, \
409 409 def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
410 410 """
411 411 Sets a number of standard patterns
412   -
  412 +
413 413 Parameters
414 414 ----------
415 415 pattern : str
... ... @@ -466,7 +466,7 @@ def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb,
466 466 deadc=0.98, pattern='single', offset=1.5):
467 467 """
468 468 Set an observation container filled with CTA observations
469   -
  469 +
470 470 Parameters
471 471 ----------
472 472 ra : float
... ... @@ -577,8 +577,17 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
577 577 else:
578 578 usepnt = False
579 579  
  580 + # Set energy limits, with larger etrue limits for the case that energy
  581 + # dispersion is requested
  582 + if edisp:
  583 + _emin = 0.5 * emin
  584 + _emax = 1.5 * emax
  585 + else:
  586 + _emin = emin
  587 + _emax = emax
  588 +
580 589 # Set number of energy bins to at least 30 per energy decade
581   - _enumbins = int((math.log10(emax) - math.log10(emin)) * 30.0)
  590 + _enumbins = int((math.log10(_emax) - math.log10(_emin)) * 30.0)
582 591 if enumbins > _enumbins:
583 592 _enumbins = enumbins
584 593  
... ... @@ -598,8 +607,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
598 607 expcube['nxpix'] = nxpix
599 608 expcube['nypix'] = nypix
600 609 expcube['enumbins'] = _enumbins
601   - expcube['emin'] = emin
602   - expcube['emax'] = emax
  610 + expcube['emin'] = _emin
  611 + expcube['emax'] = _emax
603 612 expcube['coordsys'] = coordsys
604 613 expcube['proj'] = proj
605 614 expcube['addbounds'] = addbounds
... ... @@ -621,8 +630,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
621 630 psfcube['nxpix'] = psf_nxpix
622 631 psfcube['nypix'] = psf_nypix
623 632 psfcube['enumbins'] = _enumbins
624   - psfcube['emin'] = emin
625   - psfcube['emax'] = emax
  633 + psfcube['emin'] = _emin
  634 + psfcube['emax'] = _emax
626 635 psfcube['coordsys'] = coordsys
627 636 psfcube['proj'] = proj
628 637 psfcube['addbounds'] = addbounds
... ... @@ -673,8 +682,8 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
673 682 edispcube['nxpix'] = psf_nxpix
674 683 edispcube['nypix'] = psf_nypix
675 684 edispcube['enumbins'] = _enumbins
676   - edispcube['emin'] = emin
677   - edispcube['emax'] = emax
  685 + edispcube['emin'] = _emin
  686 + edispcube['emax'] = _emax
678 687 edispcube['coordsys'] = coordsys
679 688 edispcube['proj'] = proj
680 689 edispcube['addbounds'] = addbounds
... ... @@ -723,7 +732,7 @@ def get_stacked_obs(cls, obs):
723 732 # Write header
724 733 if cls._logExplicit():
725 734 cls._log.header3('Binning events')
726   -
  735 +
727 736 # Bin events
728 737 cntcube = ctools.ctbin(obs)
729 738 cntcube['usepnt'] = False
... ...
src/ctbin/ctbin.cpp
  View file @4e9d5a1
... ... @@ -247,12 +247,21 @@ void ctbin::run(void)
247 247  
248 248 // Find all unbinned CTA observations in m_obs
249 249 std::vector<GCTAObservation*> obs_list(0);
  250 + double mean_ra = 0.0;
  251 + double mean_dec = 0.0;
250 252 for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
251 253 obs = next_unbinned_observation()) {
252 254  
253 255 // Push observation into list
254 256 obs_list.push_back(obs);
255 257  
  258 + // Get pointing
  259 + const GCTAPointing& pnt = obs->pointing();
  260 +
  261 + // Add to coordinates
  262 + mean_ra += pnt.dir().ra_deg();
  263 + mean_dec += pnt.dir().dec_deg();
  264 +
256 265 // Write message
257 266 std::string msg = " Including unbinned "+obs->instrument()+
258 267 " observation";
... ... @@ -263,6 +272,13 @@ void ctbin::run(void)
263 272 // Set number of relevant observations
264 273 int nobs = obs_list.size();
265 274  
  275 + // Compute mean pointing
  276 + if (nobs > 0) {
  277 + mean_ra /= double(nobs);
  278 + mean_dec /= double(nobs);
  279 + m_mean_pnt.radec_deg(mean_ra, mean_dec);
  280 + }
  281 +
266 282 // Initialise arrays
267 283 if (nobs > 0) {
268 284 m_counts = std::vector<GSkyMap>(nobs);
... ... @@ -275,6 +291,13 @@ void ctbin::run(void)
275 291 m_cubes.clear();
276 292 }
277 293  
  294 + // Log observation selection results
  295 + std::string use = (*this)["usepnt"].boolean() ? "yes" : "no";
  296 + log_header3(NORMAL, "Summary");
  297 + log_value(NORMAL, "Number of observations", nobs);
  298 + log_value(NORMAL, "Mean pointing", m_mean_pnt.print());
  299 + log_value(NORMAL, "Use mean pointing", use);
  300 +
278 301 // Write header into logger
279 302 log_header1(TERSE, gammalib::number("Bin observation", nobs));
280 303  
... ... @@ -584,7 +607,7 @@ void ctbin::publish(const std::string&amp; name)
584 607 void ctbin::init_members(void)
585 608 {
586 609 // Initialise members
587   - m_usepnt = false;
  610 + //m_usepnt = false;
588 611 m_stack = true;
589 612 m_prefix.clear();
590 613 m_publish = false;
... ... @@ -594,6 +617,7 @@ void ctbin::init_members(void)
594 617 m_cubes.clear();
595 618 m_counts.clear();
596 619 m_weights.clear();
  620 + m_mean_pnt.clear();
597 621 m_ebounds.clear();
598 622 m_gti.clear();
599 623 m_ontime = 0.0;
... ... @@ -615,7 +639,7 @@ void ctbin::init_members(void)
615 639 void ctbin::copy_members(const ctbin& app)
616 640 {
617 641 // Copy attributes
618   - m_usepnt = app.m_usepnt;
  642 + //m_usepnt = app.m_usepnt;
619 643 m_stack = app.m_stack;
620 644 m_prefix = app.m_prefix;
621 645 m_publish = app.m_publish;
... ... @@ -625,6 +649,7 @@ void ctbin::copy_members(const ctbin&amp; app)
625 649 m_cubes = app.m_cubes;
626 650 m_counts = app.m_counts;
627 651 m_weights = app.m_weights;
  652 + m_mean_pnt = app.m_mean_pnt;
628 653 m_ebounds = app.m_ebounds;
629 654 m_gti = app.m_gti;
630 655 m_ontime = app.m_ontime;
... ... @@ -746,7 +771,17 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs)
746 771 double xref = 0.0;
747 772 double yref = 0.0;
748 773 bool usepnt = (*this)["usepnt"].boolean();
749   - if (!usepnt) {
  774 + if (usepnt) {
  775 + if (gammalib::toupper(coordsys) == "GAL") {
  776 + xref = m_mean_pnt.l_deg();
  777 + yref = m_mean_pnt.b_deg();
  778 + }
  779 + else {
  780 + xref = m_mean_pnt.ra_deg();
  781 + yref = m_mean_pnt.dec_deg();
  782 + }
  783 + }
  784 + else {
750 785 xref = (*this)["xref"].real();
751 786 yref = (*this)["yref"].real();
752 787 }
... ... @@ -759,29 +794,10 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs)
759 794 // Read energy boundaries
760 795 GEbounds ebounds = create_ebounds();
761 796  
762   - // If requested, get pointing from observations
763   - if (usepnt) {
764   -
765   - // Get pointing direction for observation
766   - GSkyDir pnt = obs->pointing().dir();
767   -
768   - // Set xref/yref based on the coordinate system
769   - if (gammalib::toupper(coordsys) == "GAL") {
770   - xref = pnt.l_deg();
771   - yref = pnt.b_deg();
772   - }
773   - else {
774   - xref = pnt.ra_deg();
775   - yref = pnt.dec_deg();
776   - }
777   -
778   - } // endif: got pointing from observation
779   -
780 797 // Initialise counts cube
781 798 GSkyMap cube = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz,
782 799 nxpix, nypix, ebounds.size());
783 800  
784   -
785 801 // Return cube
786 802 return cube;
787 803 }
... ...
src/ctbin/ctbin.hpp
  View file @4e9d5a1
... ... @@ -99,7 +99,6 @@ protected:
99 99 void obs_cube_stacked(void);
100 100  
101 101 // User parameters
102   - bool m_usepnt; //!< Use pointing instead of xref/yref parameters
103 102 bool m_stack; //!< Output one stacked cube or multiple cubes
104 103 std::string m_prefix; //!< Prefix for multiple cubes
105 104 bool m_publish; //!< Publish counts cube?
... ... @@ -109,6 +108,7 @@ protected:
109 108 std::vector<GCTAEventCube> m_cubes; //!< Event cubes
110 109 std::vector<GSkyMap> m_counts; //!< List of event cube counts
111 110 std::vector<GSkyMap> m_weights; //!< List of event cube weights
  111 + GSkyDir m_mean_pnt; //!< Mean pointing
112 112 GEbounds m_ebounds; //!< Energy boundaries
113 113 GGti m_gti; //!< Stacked Good time intervals
114 114 double m_ontime; //!< Total ontime
... ...