GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTACubeExposure.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTACubeExposure.cpp - CTA cube analysis exposure class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-2018 by Chia-Chun Lu *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTACubeExposure.cpp
23  * @brief CTA cube analysis exposure class implementation
24  * @author Chia-Chun Lu
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GLog.hpp"
34 #include "GObservations.hpp"
35 #include "GSkyRegionCircle.hpp"
36 #include "GCTALib.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_FILL_CUBE "GCTACubeExposure::fill_cube(GCTAObservation&, GLog*)"
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 
47 /* __ Constants __________________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Void constructor
58  ***************************************************************************/
60 {
61  // Initialise class members
62  init_members();
63 
64  // Return
65  return;
66 }
67 
68 
69 /***********************************************************************//**
70  * @brief Copy constructor
71  *
72  * @param[in] cube Exposure cube.
73  ***************************************************************************/
75 {
76  // Initialise class members
77  init_members();
78 
79  // Copy members
80  copy_members(cube);
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief File constructor
89  *
90  * @param[in] filename Exposure cube filename.
91  *
92  * Construct exposure cube by loading the information from an exposure cube
93  * file.
94  ***************************************************************************/
96 {
97  // Initialise class members
98  init_members();
99 
100  // Load exposure cube from file
101  load(filename);
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Event cube constructor
110  *
111  * @param[in] cube Event cube.
112  *
113  * Construct exposure cube using the same binning and sky projection that is
114  * used for the event cube.
115  ***************************************************************************/
117 {
118  // Initialise class members
119  init_members();
120 
121  // Store energy boundaries
122  m_energies.set(cube.ebounds());
123 
124  // Set GNodeArray used for interpolation
125  set_eng_axis();
126 
127  // Set exposure cube to event cube
128  m_cube = cube.counts();
130 
131  // Set all exposure cube pixels to zero as we want to have a clean map
132  // upon construction
133  m_cube = 0.0;
134 
135  // Return
136  return;
137 
138 }
139 
140 
141 /***********************************************************************//**
142  * @brief Exposure cube constructor
143  *
144  * @param[in] wcs World Coordinate System.
145  * @param[in] coords Coordinate System (CEL or GAL).
146  * @param[in] x X coordinate of sky map centre (deg).
147  * @param[in] y Y coordinate of sky map centre (deg).
148  * @param[in] dx Pixel size in x direction at centre (deg/pixel).
149  * @param[in] dy Pixel size in y direction at centre (deg/pixel).
150  * @param[in] nx Number of pixels in x direction.
151  * @param[in] ny Number of pixels in y direction.
152  * @param[in] energies Energies.
153  *
154  * Constructs an exposure cube by specifying the sky map grid and the
155  * energies.
156  ***************************************************************************/
157 GCTACubeExposure::GCTACubeExposure(const std::string& wcs,
158  const std::string& coords,
159  const double& x,
160  const double& y,
161  const double& dx,
162  const double& dy,
163  const int& nx,
164  const int& ny,
165  const GEnergies& energies)
166 {
167  // Initialise class members
168  init_members();
169 
170  // Store energies
172 
173  // Set GNodeArray used for interpolation
174  set_eng_axis();
175 
176  // Create sky map
177  m_cube = GSkyMap(wcs, coords, x, y, dx, dy, nx, ny, m_energies.size());
178 
179  // Return
180  return;
181 }
182 
183 
184 /***********************************************************************//**
185  * @brief Destructor
186  ***************************************************************************/
188 {
189  // Free members
190  free_members();
191 
192  // Return
193  return;
194 }
195 
196 
197 /*==========================================================================
198  = =
199  = Operators =
200  = =
201  ==========================================================================*/
202 
203 /***********************************************************************//**
204  * @brief Assignment operator
205  *
206  * @param[in] cube Exposure cube.
207  * @return Exposure cube.
208  ***************************************************************************/
210 {
211  // Execute only if object is not identical
212  if (this != &cube) {
213 
214  // Free members
215  free_members();
216 
217  // Initialise private members
218  init_members();
219 
220  // Copy members
221  copy_members(cube);
222 
223  } // endif: object was not identical
224 
225  // Return this object
226  return *this;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Return exposure (in units of cm2 s)
232  *
233  * @param[in] dir Coordinate of the true photon position.
234  * @param[in] energy Energy of the true photon.
235  * @return Exposure (in units of cm2 s)
236  ***************************************************************************/
237 double GCTACubeExposure::operator()(const GSkyDir& dir, const GEnergy& energy) const
238 {
239  // Set indices and weighting factors for interpolation
240  update(energy.log10TeV());
241 
242  // Perform interpolation
243  double exposure = m_wgt_left * m_cube(dir, m_inx_left) +
245 
246  // Make sure that exposure does not become negative
247  if (exposure < 0.0) {
248  exposure = 0.0;
249  }
250 
251  // Return exposure
252  return exposure;
253 }
254 
255 
256 /*==========================================================================
257  = =
258  = Public methods =
259  = =
260  ==========================================================================*/
261 
262 /***********************************************************************//**
263  * @brief Clear instance
264  *
265  * This method properly resets the object to an initial state.
266  ***************************************************************************/
268 {
269  // Free class members
270  free_members();
271 
272  // Initialise members
273  init_members();
274 
275  // Return
276  return;
277 }
278 
279 
280 /***********************************************************************//**
281  * @brief Clone exposure cube
282  *
283  * @return Deep copy of exposure cube instance.
284  ***************************************************************************/
286 {
287  return new GCTACubeExposure(*this);
288 }
289 
290 
291 /***********************************************************************//**
292  * @brief Set exposure cube for one CTA observation
293  *
294  * @param[in] obs CTA observation.
295  *
296  * Set the exposure cube for one CTA observation. The cube pixel values are
297  * computed as product of the effective area and the livetime.
298  ***************************************************************************/
300 {
301  // Clear GTIs, reset livetime and exposure cube pixels
302  m_gti.clear();
303  m_livetime = 0.0;
304  m_cube = 0.0;
305 
306  // Fill exposure cube
307  fill_cube(obs);
308 
309  // Return
310  return;
311 }
312 
313 
314 /***********************************************************************//**
315  * @brief Fill exposure cube from observation container
316  *
317  * @param[in] obs Observation container.
318  * @param[in] log Pointer to logger.
319  *
320  * Set the exposure cube by summing the exposure for all CTA observations in
321  * an observation container. The cube pixel values are computed as the sum
322  * over the products of the effective area and the livetime.
323  ***************************************************************************/
325 {
326  // Clear GTIs, reset livetime and exposure cube pixels
327  m_gti.clear();
328  m_livetime = 0.0;
329  m_cube = 0.0;
330 
331  // Loop over all observations in container
332  for (int i = 0; i < obs.size(); ++i) {
333 
334  // Get observation and continue only if it is a CTA observation
335  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>
336  (obs[i]);
337 
338  // Skip observation if it's not CTA
339  if (cta == NULL) {
340  if (log != NULL) {
341  *log << "Skipping ";
342  *log << obs[i]->instrument();
343  *log << " observation ";
344  *log << "\"" << obs[i]->name() << "\"";
345  *log << " (id=" << obs[i]->id() << ")" << std::endl;
346  }
347  continue;
348  }
349 
350  // Skip observation if we have a binned observation
351  if (cta->eventtype() == "CountsCube") {
352  if (log != NULL) {
353  *log << "Skipping binned ";
354  *log << cta->instrument();
355  *log << " observation ";
356  *log << "\"" << cta->name() << "\"";
357  *log << " (id=" << cta->id() << ")" << std::endl;
358  }
359  continue;
360  }
361 
362  // Fill exposure cube
363  fill_cube(*cta, log);
364 
365  } // endfor: looped over observations
366 
367  // Return
368  return;
369 }
370 
371 
372 /***********************************************************************//**
373  * @brief Read exposure cube from FITS object
374  *
375  * @param[in] fits FITS object.
376  *
377  * Read the exposure cube from a FITS object.
378  ***************************************************************************/
379 void GCTACubeExposure::read(const GFits& fits)
380 {
381  // Clear object
382  clear();
383 
384  // Get HDUs
385  const GFitsImage& hdu_expcube = *fits.image("Primary");
386  const GFitsTable& hdu_energies = *fits.table(gammalib::extname_energies);
387  const GFitsTable& hdu_gti = *fits.table(gammalib::extname_gti);
388 
389  // Read cube
390  m_cube.read(hdu_expcube);
391 
392  // Read cube attributes
393  read_attributes(hdu_expcube);
394 
395  // Read energies
396  m_energies.read(hdu_energies);
397 
398  // Read GTIs
399  m_gti.read(hdu_gti);
400 
401  // Set energy node array
402  set_eng_axis();
403 
404  // Return
405  return;
406 }
407 
408 
409 /***********************************************************************//**
410  * @brief Write CTA exposure cube into FITS file.
411  *
412  * @param[in] fits FITS file.
413  *
414  * Writes the exposure cube image, the energies and the Good Time Intervals
415  * into the FITS file.
416  ***************************************************************************/
418 {
419  // Write cube
420  m_cube.write(fits);
421 
422  // Get last HDU and write attributes
423  if (fits.size() > 0) {
424  GFitsHDU& hdu = *fits[fits.size()-1];
425  write_attributes(hdu);
426  }
427 
428  // Write energies
429  m_energies.write(fits);
430 
431  // Write GTIs
432  m_gti.write(fits);
433 
434  // Return
435  return;
436 }
437 
438 
439 /***********************************************************************//**
440  * @brief Load exposure cube from FITS file
441  *
442  * @param[in] filename Performance table file name.
443  *
444  * Loads the exposure cube from a FITS file into the object.
445  ***************************************************************************/
446 void GCTACubeExposure::load(const GFilename& filename)
447 {
448  // Put into OpenMP criticial zone
449  #pragma omp critical(GCTACubeExposure_load)
450  {
451 
452  // Open FITS file
453  GFits fits(filename);
454 
455  // Read PSF cube
456  read(fits);
457 
458  // Close FITS file
459  fits.close();
460 
461  } // end of OpenMP critical zone
462 
463  // Store filename
465 
466  // Return
467  return;
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Save exposure cube into FITS file
473  *
474  * @param[in] filename Exposure cube FITS file name.
475  * @param[in] clobber Overwrite existing file?
476  *
477  * Save the exposure cube into a FITS file.
478  ***************************************************************************/
479 void GCTACubeExposure::save(const GFilename& filename, const bool& clobber) const
480 {
481  // Put into OpenMP criticial zone
482  #pragma omp critical(GCTACubeExposure_save)
483  {
484 
485  // Create FITS file
486  GFits fits;
487 
488  // Write exposure cube
489  write(fits);
490 
491  // Save FITS file
492  fits.saveto(filename, clobber);
493 
494  // Close Edisp file
495  fits.close();
496 
497  } // end of OpenMP critical zone
498 
499  // Store filename
501 
502  // Return
503  return;
504 }
505 
506 
507 /***********************************************************************//**
508  * @brief Print exposure cube information
509  *
510  * @param[in] chatter Chattiness.
511  * @return String containing exposure cube information.
512  *
513  * @todo Add content
514  ***************************************************************************/
515 std::string GCTACubeExposure::print(const GChatter& chatter) const
516 {
517  // Initialise result string
518  std::string result;
519 
520  // Continue only if chatter is not silent
521  if (chatter != SILENT) {
522 
523  // Append header
524  result.append("=== GCTACubeExposure ===");
525 
526  // Append information
527  result.append("\n"+gammalib::parformat("Filename")+m_filename);
528  result.append("\n"+gammalib::parformat("Livetime"));
529  result.append(gammalib::str(m_livetime)+" sec");
530 
531  // Append energies
532  if (m_energies.size() > 0) {
533  result.append("\n"+m_energies.print(chatter));
534  }
535  else {
536  result.append("\n"+gammalib::parformat("Energies") +
537  "not defined");
538  }
539 
540  // Append GTIs
541  if (m_gti.size() > 0) {
542  result.append("\n"+m_gti.print(chatter));
543  }
544  else {
545  result.append("\n"+gammalib::parformat("Good Time Intervals") +
546  "not defined");
547  }
548 
549  // Append skymap definition
550  result.append("\n"+m_cube.print(chatter));
551 
552  } // endif: chatter was not silent
553 
554  // Return result
555  return result;
556 }
557 
558 
559 /*==========================================================================
560  = =
561  = Private methods =
562  = =
563  ==========================================================================*/
564 
565 /***********************************************************************//**
566  * @brief Initialise class members
567  ***************************************************************************/
569 {
570  // Initialise members
571  m_filename.clear();
572  m_cube.clear();
573  m_energies.clear();
574  m_elogmeans.clear();
575  m_gti.clear();
576  m_livetime = 0.0;
577 
578  // Initialise cache
579  m_inx_left = 0;
580  m_inx_right = 0;
581  m_wgt_left = 0.0;
582  m_wgt_right = 0.0;
583 
584  // Set CTA time reference for GTIs
585  m_gti.reference(GTimeReference(G_CTA_MJDREF, "s", "TT", "LOCAL"));
586 
587  // Return
588  return;
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief Copy class members
594  *
595  * @param[in] cube Exposure cube
596  ***************************************************************************/
598 {
599  // Copy members
600  m_filename = cube.m_filename;
601  m_cube = cube.m_cube;
602  m_energies = cube.m_energies;
603  m_elogmeans = cube.m_elogmeans;
604  m_gti = cube.m_gti;
605  m_livetime = cube.m_livetime;
606 
607  // Copy cache
608  m_inx_left = cube.m_inx_left;
609  m_inx_right = cube.m_inx_right;
610  m_wgt_left = cube.m_wgt_left;
611  m_wgt_right = cube.m_wgt_right;
612 
613  // Return
614  return;
615 }
616 
617 /***********************************************************************//**
618  * @brief Delete class members
619  ***************************************************************************/
621 {
622  // Return
623  return;
624 }
625 
626 
627 /***********************************************************************//**
628  * @brief Fill exposure cube for one observation
629  *
630  * @param[in] obs Observation.
631  * @param[in] log Pointer to logger.
632  *
633  * @exception GException::invalid_value
634  * No RoI or response found in CTA observation.
635  *
636  * Fill the exposure cube from one CTA observations. The exposure cube pixel
637  * values are computed as the product of the effective area and the livetime.
638  ***************************************************************************/
640 {
641  // Only continue if we have an event list
642  if (obs.eventtype() == "EventList") {
643 
644  // Extract pointing direction, energy boundaries and ROI from
645  // observation
646  GSkyDir pnt = obs.pointing().dir();
647  GEbounds obs_ebounds = obs.ebounds();
648  GCTARoi roi = obs.roi();
649 
650  // Check for RoI sanity
651  if (!roi.is_valid()) {
652  std::string msg = "No RoI information found in "+obs.instrument()+
653  " observation \""+obs.name()+"\". Run ctselect "
654  "to specify an RoI for this observation.";
656  }
657 
658  // Convert RoI into a circular region for overlap checking
659  GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
660 
661  // Extract response from observation
662  const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
663  (obs.response());
664  if (rsp == NULL) {
665  std::string msg = "No valid instrument response function found in "+
666  obs.instrument()+" observation \""+obs.name()+
667  "\". Please specify the instrument response "
668  "function for this observation.";
670  }
671 
672  // Skip observation if livetime is zero
673  if (obs.livetime() == 0.0) {
674  if (log != NULL) {
675  *log << "Skipping unbinned ";
676  *log << obs.instrument();
677  *log << " observation ";
678  *log << "\"" << obs.name() << "\"";
679  *log << " (id=" << obs.id() << ") due to zero livetime";
680  *log << std::endl;
681  }
682  }
683 
684  // Skip observation if observation is outside the bounds of the cube
685  else if (!m_cube.overlaps(roi_reg)) {
686  if (log != NULL) {
687  *log << "Skipping unbinned ";
688  *log << obs.instrument();
689  *log << " observation ";
690  *log << "\"" << obs.name() << "\"";
691  *log << " (id=" << obs.id() << ") since it does not overlap ";
692  *log << "with the exposure cube.";
693  *log << std::endl;
694  }
695  }
696 
697  // ... otherwise continue
698  else {
699 
700  // Announce observation usage
701  if (log != NULL) {
702  *log << "Including ";
703  *log << obs.instrument();
704  *log << " observation \"" << obs.name();
705  *log << "\" (id=" << obs.id() << ")";
706  *log << " in exposure cube computation." << std::endl;
707  }
708 
709  // Loop over all pixels in sky map
710  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
711 
712  // Get pixel sky direction
713  GSkyDir dir = m_cube.inx2dir(pixel);
714 
715  // Skip pixel if it is outside the RoI
716  if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
717  continue;
718  }
719 
720  // Compute theta angle with respect to pointing direction in
721  // radians
722  double theta = pnt.dist(dir);
723 
724  // Loop over all exposure cube energies
725  for (int iebin = 0; iebin < m_energies.size(); ++iebin){
726 
727  // Get logE/TeV
728  double logE = m_energies[iebin].log10TeV();
729 
730  // Add to exposure cube (effective area * livetime)
731  m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
732  obs.livetime();
733 
734  } // endfor: looped over energy bins
735 
736  } // endfor: looped over all pixels
737 
738  // If GTI is empty then set its time reference from the
739  // observation. From then on we keep that time reference
740  if (m_gti.is_empty()) {
741  m_gti.reference(obs.gti().reference());
742  }
743 
744  // Append GTIs and increment livetime
745  m_gti.extend(obs.gti());
746  m_livetime += obs.livetime();
747 
748  } // endelse: livetime was not zero
749 
750  } // endif: observation contained an event list
751 
752  // Return
753  return;
754 }
755 
756 
757 /***********************************************************************//**
758  * @brief Update 1D cache
759  *
760  * @param[in] logE Log10 energy in TeV.
761  *
762  * Updates the 1D interpolation cache. The interpolation cache is composed
763  * of two indices and weights that define 2 data values of the 2D skymap
764  * that are used for linear interpolation.
765  *
766  * @todo Write down formula
767  ***************************************************************************/
768 void GCTACubeExposure::update(const double& logE) const
769 {
770  // Set value for node array
771  m_elogmeans.set_value(logE);
772 
773  // Set indices and weighting factors for interpolation
778 
779  // Return
780  return;
781 }
782 
783 
784 /***********************************************************************//**
785  * @brief Set nodes for a logarithmic (base 10) energy axis
786  *
787  * Set axis nodes so that each node is the logarithm of the energy values.
788  ***************************************************************************/
790 {
791  // Get number of bins
792  int bins = m_energies.size();
793 
794  // Clear node array
795  m_elogmeans.clear();
796 
797  // Compute nodes
798  for (int i = 0; i < bins; ++i) {
799 
800  // Get logE/TeV
801  m_elogmeans.append(m_energies[i].log10TeV());
802 
803  } // endfor: looped over energy bins
804 
805  // Return
806  return;
807 }
808 
809 
810 /***********************************************************************//**
811  * @brief Read exposure attributes
812  *
813  * @param[in] hdu FITS HDU.
814  *
815  * Reads CTA exposure attributes from the HDU.
816  ***************************************************************************/
818 {
819  // Read mandatory attributes
820  m_livetime = (hdu.has_card("LIVETIME")) ? hdu.real("LIVETIME") : 0.0;
821 
822  // Return
823  return;
824 }
825 
826 
827 /***********************************************************************//**
828  * @brief Write attributes to exposure extension
829  *
830  * @param[in] hdu FITS HDU.
831  ***************************************************************************/
833 {
834  // Compute some attributes
835  double tstart = m_gti.tstart().convert(m_gti.reference());
836  double tstop = m_gti.tstop().convert(m_gti.reference());
837  double telapse = m_gti.telapse();
838  double ontime = m_gti.ontime();
839  double deadc = (ontime > 0.0 && m_livetime > 0.0) ?
840  m_livetime / ontime : 1.0;
841  std::string utc_obs = m_gti.tstart().utc();
842  std::string utc_end = m_gti.tstop().utc();
843  std::string date_obs = utc_obs.substr(0, 10);
844  std::string time_obs = utc_obs.substr(11, 8);
845  std::string date_end = utc_end.substr(0, 10);
846  std::string time_end = utc_end.substr(11, 8);
847 
848  // Set observation information
849  hdu.card("CREATOR", "GammaLib", "Program which created the file");
850  hdu.card("TELESCOP", "unknown", "Telescope");
851  hdu.card("OBS_ID", "unknown", "Observation identifier");
852  hdu.card("DATE-OBS", date_obs, "Observation start date");
853  hdu.card("TIME-OBS", time_obs, "Observation start time");
854  hdu.card("DATE-END", date_end, "Observation end date");
855  hdu.card("TIME-END", time_end, "Observation end time");
856 
857  // Set observation time information
858  hdu.card("TSTART", tstart, "[s] Mission time of start of observation");
859  hdu.card("TSTOP", tstop, "[s] Mission time of end of observation");
860  m_gti.reference().write(hdu);
861  hdu.card("TELAPSE", telapse, "[s] Mission elapsed time");
862  hdu.card("ONTIME", ontime, "[s] Total good time including deadtime");
863  hdu.card("LIVETIME", m_livetime, "[s] Total livetime");
864  hdu.card("DEADC", deadc, "Deadtime correction factor");
865  hdu.card("TIMEDEL", 1.0, "Time resolution");
866 
867  // Return
868  return;
869 }
std::string print(const GChatter &chatter=NORMAL) const
Print energy container information.
Definition: GEnergies.cpp:764
void clear(void)
Clear energy container.
Definition: GEnergies.cpp:225
CTA exposure cube class.
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
double operator()(const GSkyDir &dir, const GEnergy &energy) const
Return exposure (in units of cm2 s)
void read(const GFitsTable &table)
Read energies from FITS table.
Definition: GEnergies.cpp:681
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
double m_livetime
Livetime (sec)
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
#define G_CTA_MJDREF
Reference of CTA time frame.
Definition: GCTALib.hpp:92
void write(GFits &file, const std::string &extname=gammalib::extname_energies) const
Write energies into FITS object.
Definition: GEnergies.cpp:725
void fill(const GObservations &obs, GLog *log=NULL)
Fill exposure cube from observation container.
void write(GFits &file) const
Write CTA exposure cube into FITS file.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:1944
#define G_FILL_CUBE
GEbounds ebounds(void) const
Get energy boundaries.
void clear(void)
Clear instance.
void read(const GFits &fits)
Read exposure cube from FITS object.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
virtual void response(const GResponse &rsp)
Set response function.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void counts(const GSkyMap &counts)
Set event cube counts from sky map.
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
void write_attributes(GFitsHDU &hdu) const
Write attributes to exposure extension.
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
Gammalib tools definition.
GCTACubeExposure & operator=(const GCTACubeExposure &exp)
Assignment operator.
FITS file class.
Definition: GFits.hpp:63
const std::string extname_gti
Definition: GGti.hpp:44
Interface for the circular sky region class.
Collection of CTA support header files.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
void init_members(void)
Initialise class members.
std::string print(const GChatter &chatter=NORMAL) const
Print Good Time Intervals.
Definition: GGti.cpp:988
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
void dir(const GSkyDir &dir)
Set sky direction.
GCTACubeExposure * clone(void) const
Clone exposure cube.
void update(const double &logE) const
Update 1D cache.
Energy container class.
Definition: GEnergies.hpp:60
std::string print(const GChatter &chatter=NORMAL) const
Print exposure cube information.
void fill_cube(const GCTAObservation &obs, GLog *log=NULL)
Fill exposure cube for one observation.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
void id(const std::string &id)
Set observation identifier.
void load(const GFilename &filename)
Load exposure cube from FITS file.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2557
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
void reference(const GTimeReference &ref)
Set time reference for Good Time Intervals.
Definition: GGti.hpp:256
Information logger interface definition.
Definition: GLog.hpp:62
CTA event bin container class.
virtual double livetime(void) const
Return livetime.
GCTACubeExposure(void)
Void constructor.
void clear(void)
Clear Good Time Intervals.
Definition: GGti.cpp:237
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
Energy boundaries container class.
Definition: GEbounds.hpp:60
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1321
GCTARoi roi(void) const
Get Region of Interest.
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2467
Filename class.
Definition: GFilename.hpp:62
GGti gti(void) const
Get Good Time Intervals.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
const GEnergies & energies(void) const
Return energies.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
virtual std::string instrument(void) const
Return instrument name.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
const GFilename & filename(void) const
Return exposure cube filename.
GSkyMap m_cube
Average Exposure cube.
double m_wgt_right
Weight of right node.
GFilename m_filename
Filename.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void read_attributes(const GFitsHDU &hdu)
Read exposure attributes.
GChatter
Definition: GTypemaps.hpp:33
double deadc(void) const
Return deadtime correction.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
int size(void) const
Return number of observations in container.
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:206
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:193
Observation container class.
void copy_members(const GCTACubeExposure &exp)
Copy class members.
CTA instrument response function class.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
const std::string extname_energies
Definition: GEnergies.hpp:44
void set(const GEbounds &ebounds)
Set energies from energy boundaries.
Definition: GEnergies.cpp:416
void name(const std::string &name)
Set observation name.
void free_members(void)
Delete class members.
void save(const GFilename &filename, const bool &clobber=false) const
Save exposure cube into FITS file.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
Information logger class definition.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void set(const GCTAObservation &obs)
Set exposure cube for one CTA observation.
void read(const GFitsTable &table)
Read Good Time Intervals and time reference from FITS table.
Definition: GGti.cpp:639
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
const double & ontime(void) const
Return ontime.
int m_inx_right
Index of right node.
std::string eventtype(void) const
Return event type string.
bool is_valid(void) const
Checks if RoI is valid.
Definition: GCTARoi.hpp:151
void set_eng_axis(void)
Set nodes for a logarithmic (base 10) energy axis.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
Implements a time reference.
int m_inx_left
Index of left node.
virtual ~GCTACubeExposure(void)
Destructor.
GNodeArray m_elogmeans
Mean energy for the Exposure cube.
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
const double & telapse(void) const
Returns elapsed time.
Definition: GGti.hpp:223
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
bool is_empty(void) const
Signal if there are no Good Time Intervals.
Definition: GGti.hpp:165
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Observations container class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
void write(GFits &fits, const std::string &extname=gammalib::extname_gti) const
Write Good Time Intervals and time reference into FITS object.
Definition: GGti.cpp:682
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
std::string utc(const int &precision=0) const
Return time as string in UTC time system.
Definition: GTime.cpp:464
Circular sky region class interface definition.
void extend(const GGti &gti)
Append Good Time Intervals.
Definition: GGti.cpp:511
const double & ontime(void) const
Returns ontime.
Definition: GGti.hpp:239
GEnergies m_energies
Energy values for the Exposure cube.
double convert(const GTimeReference &ref) const
Return time in specified reference.
Definition: GTime.cpp:671
GGti m_gti
Good time interval for the Exposure cube.
Mathematical function definitions.
double m_wgt_left
Weight of left node.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413