Michael Mayer / gammalib

Commit 0b61c2e92d06128862d3c4a0db871ef4da9525de

Authored by Michael Mayer
1 parent 5dae729e

Added possibility to compute rotation angles between instrument and sky system

Inline Side-by-side
Showing 11 changed files with 294 additions and 44 deletions
inst/cta/include/GCTACubeBackground.hpp
  View file @0b61c2e
... ... @@ -75,7 +75,8 @@ public:
75 75 void clear(void);
76 76 GCTACubeBackground* clone(void) const;
77 77 std::string classname(void) const;
78   - void fill(const GObservations& obs, GLog* log = NULL);
  78 + void fill(const GObservations& obs, GLog* log = NULL,
  79 + const double& delta_rot_angle = 5.0);
79 80 double integral(const double& logE) const;
80 81 void read(const GFits& fits);
81 82 void write(GFits& file) const;
... ...
inst/cta/include/GCTAEventList.hpp
  View file @0b61c2e
... ... @@ -90,6 +90,7 @@ public:
90 90 void irf_cache(const std::string& name, const int& index,
91 91 const double& irf) const;
92 92 const std::string& gtiname(void) const;
  93 + bool has_detxy(void) const;
93 94  
94 95 protected:
95 96 // Protected methods
... ... @@ -196,4 +197,16 @@ const std::string& GCTAEventList::gtiname(void) const
196 197 return (m_gti_extname);
197 198 }
198 199  
  200 +
  201 +/***********************************************************************//**
  202 + * @brief Returns whether event list has DETX/Y information
  203 + *
  204 + * @return Boolean whether event list has DETX/Y information
  205 + ***************************************************************************/
  206 +inline
  207 +bool GCTAEventList::has_detxy(void) const
  208 +{
  209 + return (m_has_detxy);
  210 +}
  211 +
199 212 #endif /* GCTAEVENTLIST_HPP */
... ...
inst/cta/include/GCTAPointing.hpp
  View file @0b61c2e
... ... @@ -41,6 +41,7 @@ class GHorizDir;
41 41 class GFitsTable;
42 42 class GXmlElement;
43 43 class GCTAInstDir;
  44 +class GCTAEventList;
44 45  
45 46  
46 47 /***********************************************************************//**
... ... @@ -75,8 +76,8 @@ public:
75 76 std::string classname(void) const;
76 77 const GSkyDir& dir(void) const;
77 78 void dir(const GSkyDir& dir);
78   - GCTAInstDir instdir(const GSkyDir& skydir) const;
79   - GSkyDir skydir(const GCTAInstDir& instdir) const;
  79 + GCTAInstDir instdir(const GSkyDir& skydir, const GTime& time = GTime()) const;
  80 + GSkyDir skydir(const GCTAInstDir& instdir, const GTime& time = GTime()) const;
80 81 const GMatrix& rot(void) const;
81 82 const double& zenith(void) const;
82 83 const double& azimuth(void) const;
... ... @@ -88,6 +89,9 @@ public:
88 89 void read(const GXmlElement& xml);
89 90 void write(GXmlElement& xml) const;
90 91 std::string print(const GChatter& chatter = NORMAL) const;
  92 + void compute_rotation_angles(const GCTAEventList* events);
  93 + const double& rotation_range(void) const;
  94 + double rotation_angle(const GTime& time) const;
91 95  
92 96 protected:
93 97 // Protected methods
... ... @@ -106,6 +110,9 @@ protected:
106 110 std::vector<double> m_table_alt; //!< Table of altitudes (rad)
107 111 GTime m_table_tmin; //!< Min time bound in table
108 112 GTime m_table_tmax; //!< Max time bound in table
  113 + GNodeArray m_time_nodes; //!< Event time nodes
  114 + std::vector<double> m_rotangles; //!< Rotation angles between sky and instrument system
  115 + double m_rot_range; //!< Range of rotation angles between sky and instrument system
109 116 GTimeReference m_reference; //!< Time reference
110 117  
111 118 // Cached members
... ... @@ -183,4 +190,16 @@ void GCTAPointing::azimuth(const double&amp; azimuth)
183 190 {
184 191 m_azimuth = azimuth;
185 192 }
  193 +
  194 +
  195 +/***********************************************************************//**
  196 + * @brief Return range of rotation angles (radians)
  197 + *
  198 + * Returns range of rotation angles (radians)
  199 + ***************************************************************************/
  200 +inline
  201 +const double& GCTAPointing::rotation_range(void) const
  202 +{
  203 + return (m_rot_range);
  204 +}
186 205 #endif /* GCTAPOINTING_HPP */
... ...
inst/cta/pyext/GCTACubeBackground.i
  View file @0b61c2e
... ... @@ -61,7 +61,8 @@ public:
61 61 void clear(void);
62 62 GCTACubeBackground* clone(void) const;
63 63 std::string classname(void) const;
64   - void fill(const GObservations& obs, GLog* log = NULL);
  64 + void fill(const GObservations& obs, GLog* log = NULL,
  65 + const double& delta_rotangle = 5.0);
65 66 double integral(const double& logE) const;
66 67 void read(const GFits& fits);
67 68 void write(GFits& file) const;
... ...
inst/cta/pyext/GCTAEventList.i
  View file @0b61c2e
... ... @@ -69,6 +69,7 @@ public:
69 69 void irf_cache(const std::string& name, const int& index,
70 70 const double& irf) const;
71 71 const std::string& gtiname(void) const;
  72 + bool has_detxy(void) const;
72 73 };
73 74  
74 75  
... ...
inst/cta/pyext/GCTAPointing.i
  View file @0b61c2e
... ... @@ -51,8 +51,10 @@ public:
51 51 std::string classname(void) const;
52 52 const GSkyDir& dir(void) const;
53 53 void dir(const GSkyDir& dir);
54   - GCTAInstDir instdir(const GSkyDir& skydir) const;
55   - GSkyDir skydir(const GCTAInstDir& instdir) const;
  54 + GCTAInstDir instdir(const GSkyDir& skydir,
  55 + const GTime& time = GTime()) const;
  56 + GSkyDir skydir(const GCTAInstDir& instdir,
  57 + const GTime& time = GTime()) const;
56 58 const GMatrix& rot(void) const;
57 59 const double& zenith(void) const;
58 60 const double& azimuth(void) const;
... ... @@ -63,6 +65,9 @@ public:
63 65 void read(const GFitsTable& table);
64 66 void read(const GXmlElement& xml);
65 67 void write(GXmlElement& xml) const;
  68 + void compute_rotation_angles(const GCTAEventList* events);
  69 + const double& rotation_range(void) const;
  70 + double rotation_angle(const GTime& time) const;
66 71 };
67 72  
68 73  
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @0b61c2e
... ... @@ -318,6 +318,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
318 318 *
319 319 * @param[in] obs Observation container.
320 320 * @param[in] log Pointer to logger (optional).
  321 + * @param[in] delta_rot_angle Minimum change in instrument rotation angle for time binning
321 322 *
322 323 * @exception GException::invalid_value
323 324 * No event list found in CTA observations.
... ... @@ -325,8 +326,12 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
325 326 * Set the background cube by computing the livetime weighted background rate
326 327 * for all CTA observations in an observation container. The cube pixel
327 328 * values are computed as the sum over the background rates.
  329 + * This method creates an internal time binning to account for the time-dependent
  330 + * rotation of the sky system with respect to the instrument system. The
  331 + * binning can be steered with the delta_rot_angle_parameter
328 332 ***************************************************************************/
329   -void GCTACubeBackground::fill(const GObservations& obs, GLog* log)
  333 +void GCTACubeBackground::fill(const GObservations& obs, GLog* log,
  334 + const double& delta_rot_angle)
330 335 {
331 336 // Clear background cube
332 337 m_cube = 0.0;
... ... @@ -388,6 +393,9 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
388 393 // Extract energy boundaries from CTA observation
389 394 GEbounds obs_ebounds = cta->ebounds();
390 395  
  396 + // Extract GTI information from CTA observation
  397 + GGti obs_gti = cta->gti();
  398 +
391 399 // Check for RoI sanity
392 400 if (!roi.is_valid()) {
393 401 std::string msg = "No RoI information found in input observation "
... ... @@ -396,44 +404,111 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
396 404 throw GException::invalid_value(G_FILL, msg);
397 405 }
398 406  
399   - // Set GTI of actual observations as the GTI of the event cube
400   - eventcube.gti(cta->gti());
  407 + // Retrieve CTA pointing
  408 + const GCTAPointing& pnt = cta->pointing();
401 409  
402   - // Get observation livetime
403   - double livetime = cta->livetime();
  410 + // Now split GTIs depending on rotation angle change:
  411 + // We now move in time steps of 1 seconds and monitor the
  412 + // change of rotation angle. If we exceed the delta_rot_angle
  413 + // parameter, a new GTI is created and stored in splitted_gti
  414 + GGti splitted_gti;
404 415  
405   - // Loop over all bins in background cube
406   - for (int i = 0; i < eventcube.size(); ++i) {
  416 + // Split GTIs only if range of angles is larger than split criterion,
  417 + // we can skip the splitting and use the observation GTIs
  418 + if (pnt.rotation_range() * gammalib::rad2deg > delta_rot_angle) {
407 419  
408   - // Get event bin
409   - GCTAEventBin* bin = eventcube[i];
  420 + // Retrieve start and stop times
  421 + GTime time_start = obs_gti.tstart();
  422 + GTime time = time_start;
  423 + GTime time_stop = obs_gti.tstop();
410 424  
411   - // Skip if energy is not contained within RoI or the energy
412   - // boundaries of the observation. Note that the contains() method
413   - // tests on bin centre value.
414   - if (!roi.contains(*bin) || !obs_ebounds.contains(bin->energy())) {
415   - continue;
416   - }
  425 + // Set time step to 1 second
  426 + double delta_time = 1.0;
  427 +
  428 + // Retrieve rotation angle from pointing and start time
  429 + double angle0 = pnt.rotation_angle(time_start);
  430 +
  431 + // Scan over time and check for rotation angle change
  432 + while(time <= time_stop) {
  433 +
  434 + // Increment time
  435 + time += delta_time;
  436 +
  437 + // Get next rotation angle with time
  438 + double angle1 = pnt.rotation_angle(time);
  439 +
  440 + // Compute difference in angles in degrees
  441 + double diff = std::abs(angle0 - angle1) * gammalib::rad2deg;
  442 +
  443 + // Check if difference exceeds criterion
  444 + if (diff > delta_rot_angle) {
  445 +
  446 + // Create new GTI
  447 + splitted_gti.append(time_start, time);
  448 +
  449 + // Set time and angle to last values and proceed
  450 + time_start = time;
  451 + angle0 = angle1;
  452 + }
  453 +
  454 + } //endwhile: move through time
  455 +
  456 + } // endif: splitting was necessary
  457 + else {
  458 +
  459 + // Recover original GTIs
  460 + splitted_gti = obs_gti;
  461 + }
  462 +
  463 + // Loop over GTIs
  464 + for (int igti = 0; igti < splitted_gti.size(); ++igti) {
  465 +
  466 + // Create new single GTI from splitted GTIs
  467 + GGti gti = GGti(splitted_gti.tstart(igti), splitted_gti.tstop(igti));
  468 +
  469 + // Set GTI of actual observations as the GTI of the event cube
  470 + eventcube.gti(gti);
  471 +
  472 + // Compute time at GTI centre
  473 + GTime tmean = gti.tstart() + 0.5 * (gti.tstop() - gti.tstart());
  474 +
  475 + // Get observation livetime of GTI
  476 + double livetime = gti.ontime() * cta->deadc(tmean);
  477 +
  478 + // Loop over all bins in background cube
  479 + for (int ibin = 0; ibin < eventcube.size(); ++ibin) {
  480 +
  481 + // Get event bin
  482 + GCTAEventBin* bin = eventcube[ibin];
  483 +
  484 + // Skip if energy is not contained within RoI or the energy
  485 + // boundaries of the observation. Note that the contains() method
  486 + // tests on bin centre value.
  487 + if (!roi.contains(*bin) || !obs_ebounds.contains(bin->energy())) {
  488 + continue;
  489 + }
  490 +
  491 + // Compute model value for event bin. The model value is
  492 + // given in counts/MeV/s/sr.
  493 + double model = obs.models().eval(*bin, *cta);
417 494  
418   - // Compute model value for event bin. The model value is
419   - // given in counts/MeV/s/sr.
420   - double model = obs.models().eval(*bin, *cta);
  495 + // Multiply by livetime to get the correct weighting for
  496 + // each observation. We divide by the total livetime later
  497 + // to get the background model in units of counts/MeV/s/sr.
  498 + model *= livetime;
421 499  
422   - // Multiply by livetime to get the correct weighting for
423   - // each observation. We divide by the total livetime later
424   - // to get the background model in units of counts/MeV/s/sr.
425   - model *= livetime;
  500 + // Add existing number of counts
  501 + model += bin->counts();
426 502  
427   - // Add existing number of counts
428   - model += bin->counts();
  503 + // Store cumulated value (units: counts/MeV/sr)
  504 + bin->counts(model);
429 505  
430   - // Store cumulated value (units: counts/MeV/sr)
431   - bin->counts(model);
  506 + } // endfor: looped over all bins
432 507  
433   - } // endfor: looped over all bins
  508 + // Accumulate livetime
  509 + total_livetime += livetime;
434 510  
435   - // Accumulate livetime
436   - total_livetime += livetime;
  511 + } // endfor: loop over chunks of GTIs
437 512  
438 513 } // endfor: looped over all observations
439 514  
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @0b61c2e
... ... @@ -319,7 +319,7 @@ double GCTAModelAeffBackground::eval(const GEvent&amp; event,
319 319 }
320 320  
321 321 // Set DETX and DETY in instrument direction
322   - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir());
  322 + GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir(), event.time());
323 323  
324 324 // Set theta and phi from instrument coordinates
325 325 double theta = std::sqrt(inst_dir.detx() * inst_dir.detx() +
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @0b61c2e
... ... @@ -319,7 +319,7 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
319 319 }
320 320  
321 321 // Set DETX and DETY in instrument direction
322   - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir());
  322 + GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir(), event.time());
323 323  
324 324 // Evaluate function
325 325 double logE = event.energy().log10TeV();
... ...
inst/cta/src/GCTAObservation.cpp
  View file @0b61c2e
... ... @@ -1031,6 +1031,10 @@ void GCTAObservation::read(const GFits&amp; fits)
1031 1031 const GFitsHDU& hdu = *fits.at(extname);
1032 1032 read_attributes(hdu);
1033 1033  
  1034 + // Compute rotation angles between instrument and sky system
  1035 + // as function of time
  1036 + m_pointing.compute_rotation_angles(events);
  1037 +
1034 1038 }
1035 1039  
1036 1040 // ... otherwise we have a binned observation
... ...
inst/cta/src/GCTAPointing.cpp
  View file @0b61c2e
... ... @@ -40,6 +40,8 @@
40 40 #include "GFitsTable.hpp"
41 41 #include "GCTAPointing.hpp"
42 42 #include "GCTAInstDir.hpp"
  43 +#include "GCTAEventList.hpp"
  44 +#include "GCTAEventAtom.hpp"
43 45  
44 46 /* __ Method name definitions ____________________________________________ */
45 47 #define G_DIR_HORIZ "GCTAPointing::dir_horiz(GTime&)"
... ... @@ -254,12 +256,16 @@ void GCTAPointing::dir(const GSkyDir&amp; dir)
254 256 * @brief Get instrument direction from sky direction
255 257 *
256 258 * @param[in] skydir Sky direction.
  259 + * @param[in] time Time
257 260 * @return CTA instrument direction.
258 261 *
259   - * Returns instrument direction as function of a sky direction
  262 + * Returns instrument direction as function of a sky direction and time
260 263 ***************************************************************************/
261   -GCTAInstDir GCTAPointing::instdir(const GSkyDir& skydir) const
  264 +GCTAInstDir GCTAPointing::instdir(const GSkyDir& skydir, const GTime& time) const
262 265 {
  266 + // Get rotation angle
  267 + double rotangle = time.secs() != 0.0 ? rotation_angle(time) : 0.0;
  268 +
263 269 #if defined(G_USE_VECTORS)
264 270 // Compute rotation matrix
265 271 update();
... ... @@ -279,13 +285,13 @@ GCTAInstDir GCTAPointing::instdir(const GSkyDir&amp; skydir) const
279 285  
280 286 // Check if theta and phi are defined
281 287 if (theta > 0.0 ) {
282   - double phi = std::asin(inst[1] / std::sin(theta));
  288 + double phi = std::asin(inst[1] / std::sin(theta)) - rotangle;
283 289 detx = theta * std::cos(phi);
284 290 dety = theta * std::sin(phi);
285 291 }
286 292 #else
287 293 double theta = m_dir.dist(skydir);
288   - double phi = m_dir.posang(skydir);
  294 + double phi = m_dir.posang(skydir) - rotangle;
289 295 double detx(0.0);
290 296 double dety(0.0);
291 297 if (theta > 0.0 ) {
... ... @@ -308,12 +314,16 @@ GCTAInstDir GCTAPointing::instdir(const GSkyDir&amp; skydir) const
308 314 * @brief Get sky direction direction from instrument direction
309 315 *
310 316 * @param[in] instdir Instrument direction.
  317 + * @param[in] time Time
311 318 * @return Sky direction.
312 319 *
313   - * Returns sky direction as function of an instrument direction
  320 + * Returns sky direction as function of an instrument direction and time
314 321 ***************************************************************************/
315   -GSkyDir GCTAPointing::skydir(const GCTAInstDir& instdir)const
  322 +GSkyDir GCTAPointing::skydir(const GCTAInstDir& instdir, const GTime& time)const
316 323 {
  324 + // Get rotation angle
  325 + double rotangle = time.secs() != 0.0 ? rotation_angle(time) : 0.0;
  326 +
317 327 // Compute rotation matrix
318 328 update();
319 329  
... ... @@ -323,7 +333,7 @@ GSkyDir GCTAPointing::skydir(const GCTAInstDir&amp; instdir)const
323 333  
324 334 // Convert to polar coordinates
325 335 double theta = std::sqrt(inst_x * inst_x + inst_y * inst_y);
326   - double phi = std::atan2(inst_y, inst_x);
  336 + double phi = std::atan2(inst_y, inst_x) + rotangle;
327 337 double sin_phi = std::sin(phi);
328 338 double cos_phi = std::cos(phi);
329 339 double sin_theta = std::sin(theta);
... ... @@ -603,6 +613,121 @@ std::string GCTAPointing::print(const GChatter&amp; chatter) const
603 613 return result;
604 614 }
605 615  
  616 +/***********************************************************************//**
  617 + * @brief Compute rotation angles between sky and instrument directions
  618 + * as a function of time
  619 + *
  620 + * @param[in] events CTA Event list
  621 + *
  622 + * Computes the rotation angles between sky and instrument coordinate system
  623 + * as a function of time
  624 + ***************************************************************************/
  625 +void GCTAPointing::compute_rotation_angles(const GCTAEventList* events)
  626 +{
  627 + // Continue only if DETX and DETY information were found
  628 + if (events->has_detxy()) {
  629 +
  630 + // Clear rotation angle node array
  631 + m_time_nodes.clear();
  632 + m_rotangles.clear();
  633 +
  634 + // Initialise node array vectors
  635 + std::vector<double> nodes;
  636 +
  637 + // Initialise minimum and maximum rotation angles for later use
  638 + double rot_min(1e3), rot_max(-1e3);
  639 +
  640 + // Loop over events
  641 + for(int i = 0; i < events->size(); ++i) {
  642 +
  643 + // Get event atom
  644 + const GCTAEventAtom* atom = events->operator[](i);
  645 +
  646 + // Check for valid event
  647 + if (atom != NULL) {
  648 +
  649 + // Retrieve DETX and DETY values
  650 + double detx = atom->dir().detx();
  651 + double dety = atom->dir().dety();
  652 +
  653 + // Skip event if both are zero
  654 + if(detx == 0.0 && dety == 0.0) {
  655 + continue;
  656 + }
  657 +
  658 + // Compute rotation angle in instrument system
  659 + double instang = std::atan2(dety, detx);
  660 +
  661 + // Compute position angle in sky system
  662 + double posang = m_dir.posang(atom->dir().dir());
  663 +
  664 + // Compute relative rotation angle between two systems
  665 + double rot_angle = posang - instang;
  666 +
  667 + // Map rotation angles into continuous range
  668 + // This prevents jumps at rotation angles close to zero
  669 + // Todo: make sure that angles never jump from -180 to 180 degrees
  670 + // Otherwise rotation angle range is not correct
  671 + if (rot_angle < -gammalib::pi) {
  672 + rot_angle += gammalib::twopi;
  673 + }
  674 + else if (rot_angle > gammalib::pi) {
  675 + rot_angle -= gammalib::twopi;
  676 + }
  677 +
  678 + // Update min and max rotation angles
  679 + if (rot_angle < rot_min) {
  680 + rot_min = rot_angle;
  681 + }
  682 + if (rot_angle > rot_max) {
  683 + rot_max = rot_angle;
  684 + }
  685 +
  686 + // Update rotation angle range
  687 + m_rot_range = rot_max - rot_min;
  688 +
  689 + // Set node to time in seconds
  690 + nodes.push_back(atom->time().secs());
  691 +
  692 + // Append to rotation angle vector
  693 + m_rotangles.push_back(rot_angle);
  694 +
  695 + } // endif: atom was available
  696 +
  697 + } // endfor: loop over event list
  698 +
  699 + // Assign nodes to array
  700 + m_time_nodes.nodes(nodes);
  701 +
  702 + } // endif: events had detector coordinates information
  703 +
  704 +}
  705 +
  706 +
  707 +/***********************************************************************//**
  708 + * @brief Returns rotation angle between sky and instrument directions
  709 + * as a function of time
  710 + *
  711 + * @param[in] time Time at which rotation angle should be returned
  712 + * @return Rotation angle as function of time [radians]
  713 + *
  714 + * Returns the rotation angles between sky and instrument coordinate system
  715 + * as a function of time
  716 + * todo: Should we throw an exception if time is out of range of m_time_nodes?
  717 + ***************************************************************************/
  718 +double GCTAPointing::rotation_angle(const GTime& time) const
  719 +{
  720 + // Compute rotation angle via interpolation
  721 + double rotation_angle = 0.0;
  722 +
  723 + // Interpolation if possible
  724 + if (m_time_nodes.size() > 1) {
  725 + rotation_angle = m_time_nodes.interpolate(time.secs(), m_rotangles);
  726 + }
  727 +
  728 + // Return rotation angle
  729 + return (rotation_angle);
  730 +}
606 731  
607 732 /*==========================================================================
608 733 = =
... ... @@ -625,6 +750,9 @@ void GCTAPointing::init_members(void)
625 750 m_table_alt.clear();
626 751 m_table_tmin.clear();
627 752 m_table_tmax.clear();
  753 + m_time_nodes.clear();
  754 + m_rotangles.clear();
  755 + m_rot_range = 0.0;
628 756 m_reference.clear();
629 757  
630 758 // Initialise cache
... ... @@ -653,6 +781,9 @@ void GCTAPointing::copy_members(const GCTAPointing&amp; pnt)
653 781 m_table_alt = pnt.m_table_alt;
654 782 m_table_tmin = pnt.m_table_tmin;
655 783 m_table_tmax = pnt.m_table_tmax;
  784 + m_time_nodes = pnt.m_time_nodes;
  785 + m_rotangles = pnt.m_rotangles;
  786 + m_rot_range = pnt.m_rot_range;
656 787 m_reference = pnt.m_reference;
657 788  
658 789 // Copy cache
... ...