GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTACubePsf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTACubePsf.cpp - CTA cube analysis point spread function 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 GCTACubePsf.cpp
23  * @brief CTA cube analysis point spread function 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 "GCTACubePsf.hpp"
37 #include "GCTAObservation.hpp"
38 #include "GCTAResponseIrf.hpp"
39 #include "GCTAEventList.hpp"
40 #include "GCTAEventCube.hpp"
41 
42 /* __ Method name definitions ____________________________________________ */
43 #define G_SET "GCTACubePsf::set(GCTAObservation&)"
44 #define G_FILL "GCTACubePsf::fill(GObservations&, GLog*)"
45 #define G_FILL_CUBE "GCTACubePsf::fill_cube(GCTAObservation&, GSkyMap*, "\
46  "GLog*)"
47 
48 /* __ Macros _____________________________________________________________ */
49 
50 /* __ Coding definitions _________________________________________________ */
51 #define G_SMOOTH_PSF //!< Guarantee no singularities
52 #define G_QUADRATIC_BINNING //!< Use quadratic delta spacing
53 
54 /* __ Debug definitions __________________________________________________ */
55 
56 /* __ Constants __________________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise class members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief Copy constructor
80  *
81  * @param[in] cube Point spread function.
82  ***************************************************************************/
84 {
85  // Initialise class members
86  init_members();
87 
88  // Copy members
89  copy_members(cube);
90 
91  // Return
92  return;
93 }
94 
95 
96 /***********************************************************************//**
97  * @brief File constructor
98  *
99  * @param[in] filename PSF cube filename.
100  *
101  * Construct PSF cube by loading the information from a PSF cube file.
102  ***************************************************************************/
104 {
105  // Initialise class members
106  init_members();
107 
108  // Load PSF cube from file
109  load(filename);
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief Event cube constructor
118  *
119  * @param[in] cube Event cube.
120  * @param[in] dmax Maximum delta (deg).
121  * @param[in] ndbins Number of delta bins.
122  *
123  * Construct PSF cube using the same binning and sky projection that is
124  * used for the event cube.
125  ***************************************************************************/
126 GCTACubePsf::GCTACubePsf(const GCTAEventCube& cube, const double& dmax,
127  const int& ndbins)
128 {
129  // Initialise class members
130  init_members();
131 
132  // Store energy boundaries
133  m_energies.set(cube.ebounds());
134 
135  // Set GNodeArray used for interpolation
136  set_eng_axis();
137 
138  // Set delta node array
139  m_deltas.clear();
140  for (int i = 0; i < ndbins; ++i) {
141  #if defined(G_QUADRATIC_BINNING)
142  double binsize = std::sqrt(dmax) / double(ndbins);
143  double delta = binsize * (double(i) + 0.5); // avoid central singularity
144  delta *= delta;
145  #else
146  double binsize = dmax / double(ndbins);
147  double delta = binsize * (double(i) + 0.5); // avoid central singularity
148  #endif
149  m_deltas.append(delta);
150  }
151 
152  // Set delta node array for computation
153  set_delta_axis();
154 
155  // Compute number of sky maps
156  int nmaps = m_energies.size() * m_deltas.size();
157 
158  // Set PSF cube to event cube
159  m_cube = cube.counts();
160 
161  // Set appropriate number of skymaps
162  m_cube.nmaps(nmaps);
163 
164  // Set cube shape
166 
167  // Set all PSF cube pixels to zero as we want to have a clean map
168  // upon construction
169  m_cube = 0.0;
170 
171  // Return
172  return;
173 
174 }
175 
176 
177 /***********************************************************************//**
178  * @brief PSF cube constructor
179  *
180  * @param[in] wcs World Coordinate System.
181  * @param[in] coords Coordinate System (CEL or GAL).
182  * @param[in] x X coordinate of sky map centre (deg).
183  * @param[in] y Y coordinate of sky map centre (deg).
184  * @param[in] dx Pixel size in x direction at centre (deg/pixel).
185  * @param[in] dy Pixel size in y direction at centre (deg/pixel).
186  * @param[in] nx Number of pixels in x direction.
187  * @param[in] ny Number of pixels in y direction.
188  * @param[in] energies Energies.
189  * @param[in] dmax Maximum delta (deg).
190  * @param[in] ndbins Number of delta bins.
191  *
192  * Constructs a PSF cube by specifying the sky map grid and the energies.
193  ***************************************************************************/
194 GCTACubePsf::GCTACubePsf(const std::string& wcs,
195  const std::string& coords,
196  const double& x,
197  const double& y,
198  const double& dx,
199  const double& dy,
200  const int& nx,
201  const int& ny,
202  const GEnergies& energies,
203  const double& dmax,
204  const int& ndbins)
205 {
206  // Initialise class members
207  init_members();
208 
209  // Store energies
211 
212  // Set energy node array
213  set_eng_axis();
214 
215  // Set delta node array
216  m_deltas.clear();
217  for (int i = 0; i < ndbins; ++i) {
218  #if defined(G_QUADRATIC_BINNING)
219  double binsize = std::sqrt(dmax) / double(ndbins);
220  double delta = binsize * (double(i) + 0.5); // avoid central singularity
221  delta *= delta;
222  #else
223  double binsize = dmax / double(ndbins);
224  double delta = binsize * (double(i) + 0.5); // avoid central singularity
225  #endif
226  m_deltas.append(delta);
227  }
228 
229  // Set delta node array for computation
230  set_delta_axis();
231 
232  // Compute number of sky maps
233  int nmaps = m_energies.size() * m_deltas.size();
234 
235  // Create sky map
236  m_cube = GSkyMap(wcs, coords, x, y, dx, dy, nx, ny, nmaps);
237 
238  // Set cube shape
240 
241  // Return
242  return;
243 }
244 
245 
246 /***********************************************************************//**
247  * @brief Destructor
248  ***************************************************************************/
250 {
251  // Free members
252  free_members();
253 
254  // Return
255  return;
256 }
257 
258 
259 /*==========================================================================
260  = =
261  = Operators =
262  = =
263  ==========================================================================*/
264 
265 /***********************************************************************//**
266  * @brief Assignment operator
267  *
268  * @param[in] cube Mean PSF cube.
269  * @return Mean PSF cube.
270  ***************************************************************************/
272 {
273  // Execute only if object is not identical
274  if (this != &cube) {
275 
276  // Free members
277  free_members();
278 
279  // Initialise private members
280  init_members();
281 
282  // Copy members
283  copy_members(cube);
284 
285  } // endif: object was not identical
286 
287  // Return this object
288  return *this;
289 }
290 
291 
292 /***********************************************************************//**
293  * @brief Return point spread function (in units of sr^-1)
294  *
295  * @param[in] dir Coordinate of the true photon position.
296  * @param[in] delta Angular separation between true and measured photon
297  * directions (rad).
298  * @param[in] energy Energy of the true photon.
299  * @return point spread function (in units of sr^-1)
300  *
301  * Returns the point spread function for a given angular separation in units
302  * of sr^-1 for a given energy and coordinate.
303  ***************************************************************************/
304 double GCTACubePsf::operator()(const GSkyDir& dir,
305  const double& delta,
306  const GEnergy& energy) const
307 {
308  // Update indices and weighting factors for interpolation
309  update(delta, energy.log10TeV());
310 
311  // Perform bi-linear interpolation
312  double psf = m_wgt1 * m_cube(dir, m_inx1) +
313  m_wgt2 * m_cube(dir, m_inx2) +
314  m_wgt3 * m_cube(dir, m_inx3) +
315  m_wgt4 * m_cube(dir, m_inx4);
316 
317  // Make sure that PSF does not become negative
318  if (psf < 0.0) {
319  psf = 0.0;
320  }
321 
322  // Return PSF
323  return psf;
324 }
325 
326 
327 /*==========================================================================
328  = =
329  = Public methods =
330  = =
331  ==========================================================================*/
332 
333 /***********************************************************************//**
334  * @brief Clear instance
335  *
336  * This method properly resets the object to an initial state.
337  ***************************************************************************/
339 {
340  // Free class members
341  free_members();
342 
343  // Initialise members
344  init_members();
345 
346  // Return
347  return;
348 }
349 
350 
351 /***********************************************************************//**
352  * @brief Clone instance
353  *
354  * @return Deep copy of mean PSF instance.
355  ***************************************************************************/
357 {
358  return new GCTACubePsf(*this);
359 }
360 
361 
362 /***********************************************************************//**
363  * @brief Set PSF cube from one CTA observation
364  *
365  * @param[in] obs CTA observation.
366  *
367  * Sets the PSF cube for one CTA observation.
368  ***************************************************************************/
370 {
371  // Clear PSF cube
372  clear_cube();
373 
374  // Fill PSF cube
375  fill_cube(obs);
376 
377  // Compile option: guarantee smooth Psf
378  #if defined(G_SMOOTH_PSF)
379  set_to_smooth();
380  #endif
381 
382  // Return
383  return;
384 }
385 
386 
387 /***********************************************************************//**
388  * @brief Fill PSF cube from observation container
389  *
390  * @param[in] obs Observation container.
391  * @param[in] log Pointer towards logger.
392  *
393  * Sets the PSF cube from all CTA observations in the observation container.
394  ***************************************************************************/
396 {
397  // Clear PSF cube
398  clear_cube();
399 
400  // Initialise skymap for exposure weight accumulation
401  GSkyMap exposure(m_cube);
402 
403  // Loop over all observations in container
404  for (int i = 0; i < obs.size(); ++i) {
405 
406  // Get observation and continue only if it is a CTA observation
407  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>
408  (obs[i]);
409 
410  // Skip observation if it's not CTA
411  if (cta == NULL) {
412  if (log != NULL) {
413  *log << "Skipping ";
414  *log << obs[i]->instrument();
415  *log << " observation ";
416  *log << "\"" << obs[i]->name() << "\"";
417  *log << " (id=" << obs[i]->id() << ")" << std::endl;
418  }
419  continue;
420  }
421 
422  // Skip observation if we have a binned observation
423  if (cta->eventtype() == "CountsCube") {
424  if (log != NULL) {
425  *log << "Skipping binned ";
426  *log << cta->instrument();
427  *log << " observation ";
428  *log << "\"" << cta->name() << "\"";
429  *log << " (id=" << cta->id() << ")" << std::endl;
430  }
431  continue;
432  }
433 
434  // Fill PSF cube
435  fill_cube(*cta, &exposure, log);
436 
437  } // endfor: looped over observations
438 
439  // Compute mean PSF cube by dividing though the weights
440  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
441  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
442  if (exposure(pixel, iebin) > 0.0) {
443  double norm = 1.0 / exposure(pixel, iebin);
444  for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
445  int imap = offset(idelta, iebin);
446  m_cube(pixel, imap) *= norm;
447  }
448  }
449  else {
450  for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
451  int imap = offset(idelta, iebin);
452  m_cube(pixel, imap) = 0.0;
453  }
454  }
455  }
456  }
457 
458  // Compile option: guarantee smooth Psf
459  #if defined(G_SMOOTH_PSF)
460  set_to_smooth();
461  #endif
462 
463  // Return
464  return;
465 }
466 
467 
468 /***********************************************************************//**
469  * @brief Read PSF cube from FITS object
470  *
471  * @param[in] fits FITS object.
472  *
473  * Read the PSF cube from a FITS object.
474  ***************************************************************************/
475 void GCTACubePsf::read(const GFits& fits)
476 {
477  // Clear object
478  clear();
479 
480  // Get HDUs
481  const GFitsImage& hdu_psfcube = *fits.image("Primary");
482  const GFitsTable& hdu_energies = *fits.table(gammalib::extname_energies);
483  const GFitsTable& hdu_deltas = *fits.table(gammalib::extname_deltas);
484 
485  // Read cube
486  m_cube.read(hdu_psfcube);
487 
488  // Read energies
489  m_energies.read(hdu_energies);
490 
491  // Read delta nodes
492  m_deltas.read(hdu_deltas);
493 
494  // Set delta node array for computation
495  set_delta_axis();
496 
497  // Set energy node array
498  set_eng_axis();
499 
500  // Compile option: guarantee smooth Psf
501  #if defined(G_SMOOTH_PSF)
502  set_to_smooth();
503  #endif
504 
505  // Return
506  return;
507 }
508 
509 
510 /***********************************************************************//**
511  * @brief Write CTA PSF cube into FITS object.
512  *
513  * @param[in] fits FITS object.
514  *
515  * Write the CTA PSF cube into a FITS object.
516  ***************************************************************************/
517 void GCTACubePsf::write(GFits& fits) const
518 {
519  // Write cube
520  m_cube.write(fits);
521 
522  // Write energies
523  m_energies.write(fits);
524 
525  // Write delta nodes
527 
528  // Set the nodes unit to "deg"
529  (*fits.table(gammalib::extname_deltas))["Value"]->unit("deg");
530 
531  // Return
532  return;
533 }
534 
535 
536 /***********************************************************************//**
537  * @brief Load PSF cube from FITS file
538  *
539  * @param[in] filename Performance table file name.
540  *
541  * Loads the PSF cube from a FITS file into the object.
542  ***************************************************************************/
543 void GCTACubePsf::load(const GFilename& filename)
544 {
545  // Put into OpenMP criticial zone
546  #pragma omp critical(GCTACubePsf_load)
547  {
548 
549  // Open FITS file
550  GFits fits(filename);
551 
552  // Read PSF cube
553  read(fits);
554 
555  // Close FITS file
556  fits.close();
557 
558  // Store filename
560 
561  } // end of OpenMP critical zone
562 
563  // Return
564  return;
565 }
566 
567 
568 /***********************************************************************//**
569  * @brief Save PSF cube into FITS file
570  *
571  * @param[in] filename PSF cube FITS file name.
572  * @param[in] clobber Overwrite existing file?
573  *
574  * Save the PSF cube into a FITS file.
575  ***************************************************************************/
576 void GCTACubePsf::save(const GFilename& filename, const bool& clobber) const
577 {
578  // Put into OpenMP criticial zone
579  #pragma omp critical(GCTACubePsf_save)
580  {
581 
582  // Create FITS file
583  GFits fits;
584 
585  // Write PSF cube
586  write(fits);
587 
588  // Save FITS file
589  fits.saveto(filename, clobber);
590 
591  // Close Edisp file
592  fits.close();
593 
594  } // end of OpenMP critical zone
595 
596  // Store filename
598 
599  // Return
600  return;
601 }
602 
603 
604 /***********************************************************************//**
605  * @brief Print PSF cube information
606  *
607  * @param[in] chatter Chattiness.
608  * @return String containing PSF cube information.
609  ***************************************************************************/
610 std::string GCTACubePsf::print(const GChatter& chatter) const
611 {
612  // Initialise result string
613  std::string result;
614 
615  // Continue only if chatter is not silent
616  if (chatter != SILENT) {
617 
618  // Append header
619  result.append("=== GCTACubePsf ===");
620 
621  // Append information
622  result.append("\n"+gammalib::parformat("Filename")+m_filename);
623 
624  // Append energies
625  if (m_energies.size() > 0) {
626  result.append("\n"+m_energies.print(chatter));
627  }
628  else {
629  result.append("\n"+gammalib::parformat("Energies") +
630  "not defined");
631  }
632 
633  // Append number of delta bins
634  result.append("\n"+gammalib::parformat("Number of delta bins") +
636 
637  // Append delta range
638  result.append("\n"+gammalib::parformat("Delta range"));
639  if (m_deltas.size() > 0) {
640  result.append(gammalib::str(m_deltas[0]));
641  result.append(" - ");
642  result.append(gammalib::str(m_deltas[m_deltas.size()-1]));
643  result.append(" deg");
644  }
645  else {
646  result.append("not defined");
647  }
648 
649  // Append skymap definition
650  result.append("\n"+m_cube.print(chatter));
651 
652  } // endif: chatter was not silent
653 
654  // Return result
655  return result;
656 }
657 
658 
659 /*==========================================================================
660  = =
661  = Private methods =
662  = =
663  ==========================================================================*/
664 
665 /***********************************************************************//**
666  * @brief Initialise class members
667  ***************************************************************************/
669 {
670  // Initialise members
671  m_filename.clear();
672  m_cube.clear();
673  m_energies.clear();
674  m_elogmeans.clear();
675  m_deltas.clear();
677  m_quadratic_binning = false;
678 
679  // Initialise cache
680  m_inx1 = 0;
681  m_inx2 = 0;
682  m_inx3 = 0;
683  m_inx4 = 0;
684  m_wgt1 = 0.0;
685  m_wgt2 = 0.0;
686  m_wgt3 = 0.0;
687  m_wgt4 = 0.0;
688 
689  // Return
690  return;
691 }
692 
693 
694 /***********************************************************************//**
695  * @brief Copy class members
696  *
697  * @param[in] cube PSF cube.
698  ***************************************************************************/
700 {
701  // Copy members
702  m_filename = cube.m_filename;
703  m_cube = cube.m_cube;
704  m_energies = cube.m_energies;
705  m_elogmeans = cube.m_elogmeans;
706  m_deltas = cube.m_deltas;
709 
710  // Copy cache
711  m_inx1 = cube.m_inx1;
712  m_inx2 = cube.m_inx2;
713  m_inx3 = cube.m_inx3;
714  m_inx4 = cube.m_inx4;
715  m_wgt1 = cube.m_wgt1;
716  m_wgt2 = cube.m_wgt2;
717  m_wgt3 = cube.m_wgt3;
718  m_wgt4 = cube.m_wgt4;
719 
720  // Return
721  return;
722 }
723 
724 /***********************************************************************//**
725  * @brief Delete class members
726  ***************************************************************************/
728 {
729  // Return
730  return;
731 }
732 
733 
734 /***********************************************************************//**
735  * @brief Clear all pixels in the PSF cube
736  ***************************************************************************/
738 {
739  // Loop over all maps
740  for (int imap = 0; imap < m_cube.nmaps(); ++imap) {
741 
742  // Loop over all pixels in sky map
743  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
744 
745  // Reset cube value to zero
746  m_cube(pixel, imap) = 0.0;
747 
748  } // endfor: looped over all pixels
749 
750  } // endfor: looped over maps
751 
752  // Return
753  return;
754 }
755 
756 
757 /***********************************************************************//**
758  * @brief Fill PSF cube for one observation
759  *
760  * @param[in] obs Observation.
761  * @param[in] exposure Pointer towards exposure map.
762  * @param[in] log Pointer towards logger.
763  *
764  * @exception GException::invalid_value
765  * No RoI or response found in CTA observation.
766  ***************************************************************************/
768  GSkyMap* exposure,
769  GLog* log)
770 {
771  // Only continue if we have an event list
772  if (obs.eventtype() == "EventList") {
773 
774  // Extract pointing direction, energy boundaries and ROI from
775  // observation
776  GSkyDir pnt = obs.pointing().dir();
777  GEbounds obs_ebounds = obs.ebounds();
778  GCTARoi roi = obs.roi();
779 
780  // Check for RoI sanity
781  if (!roi.is_valid()) {
782  std::string msg = "No RoI information found in "+obs.instrument()+
783  " observation \""+obs.name()+"\". Run ctselect "
784  "to specify an RoI for this observation.";
786  }
787 
788  // Convert RoI into a circular region for overlap checking
789  GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
790 
791  // Extract response from observation
792  const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
793  (obs.response());
794  if (rsp == NULL) {
795  std::string msg = "No valid instrument response function found in "+
796  obs.instrument()+" observation \""+obs.name()+
797  "\". Please specify the instrument response "
798  "function for this observation.";
800  }
801 
802  // Skip observation if livetime is zero
803  if (obs.livetime() == 0.0) {
804  if (log != NULL) {
805  *log << "Skipping unbinned ";
806  *log << obs.instrument();
807  *log << " observation ";
808  *log << "\"" << obs.name() << "\"";
809  *log << " (id=" << obs.id() << ") due to zero livetime";
810  *log << std::endl;
811  }
812  }
813 
814  // Skip observation if observation is outside the bounds of the cube
815  else if (!m_cube.overlaps(roi_reg)) {
816  if (log != NULL) {
817  *log << "Skipping unbinned ";
818  *log << obs.instrument();
819  *log << " observation ";
820  *log << "\"" << obs.name() << "\"";
821  *log << " (id=" << obs.id() << ") since it does not overlap ";
822  *log << "with the point spread function cube.";
823  *log << std::endl;
824  }
825  }
826 
827  // ... otherwise continue
828  else {
829 
830  // Announce observation usage
831  if (log != NULL) {
832  *log << "Including ";
833  *log << obs.instrument();
834  *log << " observation \"" << obs.name();
835  *log << "\" (id=" << obs.id() << ")";
836  *log << " in point spread function cube computation." << std::endl;
837  }
838 
839  // Loop over all pixels in sky map
840  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
841 
842  // Get pixel sky direction
843  GSkyDir dir = m_cube.inx2dir(pixel);
844 
845  // Skip pixel if it is outside the RoI
846  if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
847  continue;
848  }
849 
850  // Compute theta angle with respect to pointing direction in
851  // radians
852  double theta = pnt.dist(dir);
853 
854  // Loop over all energy bins
855  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
856 
857  // Get logE/TeV
858  double logE = m_energies[iebin].log10TeV();
859 
860  // Compute exposure weight. If no exposure map is
861  // specified, the weight is one. Otherwise, the weight
862  // is the effective area for this offset angle and
863  // source energy times the livetime of the observation.
864  double weight = 1.0;
865  if (exposure != NULL) {
866  weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
867  obs.livetime();
868  (*exposure)(pixel, iebin) += weight;
869  }
870 
871  // Loop over delta values
872  for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
873 
874  // Compute delta in radians
875  double delta = m_deltas[idelta] * gammalib::deg2rad;
876 
877  // Set map index
878  int imap = offset(idelta, iebin);
879 
880  // Add on PSF cube
881  m_cube(pixel, imap) +=
882  rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
883 
884  } // endfor: looped over delta bins
885 
886  } // endfor: looped over energy bins
887 
888  } // endfor: looped over all pixels
889 
890  } // endelse: livetime was not zero
891 
892  } // endif: observation contained an event list
893 
894  // Return
895  return;
896 }
897 
898 
899 /***********************************************************************//**
900  * @brief Update PSF parameter cache
901  *
902  * @param[in] delta Angular separation between true and measured photon
903  * directions (radians).
904  * @param[in] logE Log10 true photon energy (TeV).
905  *
906  * This method updates the PSF parameter cache.
907  ***************************************************************************/
908 void GCTACubePsf::update(const double& delta, const double& logE) const
909 {
910  // Set node array for delta interpolation
911  if (m_quadratic_binning) {
913  }
914  else {
915  m_deltas_cache.set_value(delta);
916  }
917 
918  // Set node array for energy interpolation
919  m_elogmeans.set_value(logE);
920 
921  // Set indices for bi-linear interpolation
926 
927  // Set weighting factors for bi-linear interpolation
932 
933  // Return
934  return;
935 }
936 
937 
938 /***********************************************************************//**
939  * @brief Set nodes for delta axis in radians
940  *
941  * Set the delta axis nodes in radians. If a quadratic binning is detected,
942  * the axis is set to sqrt(delta) so that a linear interpolation scheme
943  * can be used which is much faster than a quadratic scheme.
944  *
945  * @todo Check that none of the axis boundaries is non-positive.
946  ***************************************************************************/
948 {
949  // Initialise computation members
951  m_quadratic_binning = false;
952 
953  // Set up a linear array that represents the square root of the delta
954  // values. If this array corresponds within some numerical precision
955  // to the actual delta values, will will use these values for the
956  // interpolation (the m_quadratic_binning member will be set to true).
957  // Otherwise, the original delta values will be kept and the
958  // m_quadratic_binning member will be set to false.
959  int ndbins = m_deltas.size();
960  if (ndbins > 1) {
961  m_quadratic_binning = true;
962  double binsize = (std::sqrt(m_deltas[ndbins-1] * gammalib::deg2rad) -
964  double(ndbins-1);
965  for (int i = 0; i < ndbins; ++i) {
966  double delta = binsize * (double(i) + 0.5);
967  double diff = std::abs(delta*delta-m_deltas[i] * gammalib::deg2rad);
968  if (diff > 1.0e-6) {
969  m_quadratic_binning = false;
970  break;
971  }
972  m_deltas_cache.append(delta);
973  }
974  }
975 
976  // If we do not have square root binning then use the original values
977  // but convert them to radians
978  if (!m_quadratic_binning) {
980  for (int i = 0; i < m_deltas.size(); ++i) {
982  }
983  }
984 
985  // Return
986  return;
987 }
988 
989 
990 /***********************************************************************//**
991  * @brief Set nodes for a logarithmic (base 10) energy axis
992  *
993  *
994  * Set axis nodes so that each node is the logarithmic mean of the lower and
995  * upper energy boundary, i.e.
996  * \f[ n_i = \log \sqrt{{\rm LO}_i \times {\rm HI}_i} \f]
997  * where
998  * \f$n_i\f$ is node \f$i\f$,
999  * \f${\rm LO}_i\f$ is the lower bin boundary for bin \f$i\f$, and
1000  * \f${\rm HI}_i\f$ is the upper bin boundary for bin \f$i\f$.
1001  *
1002  * @todo Check that none of the axis boundaries is non-positive.
1003  ***************************************************************************/
1005 {
1006  // Get number of bins
1007  int bins = m_energies.size();
1008 
1009  // Clear node array
1010  m_elogmeans.clear();
1011 
1012  // Compute nodes
1013  for (int i = 0; i < bins; ++i) {
1014 
1015  // Get logE/TeV
1016  m_elogmeans.append(m_energies[i].log10TeV());
1017 
1018  } // endfor: looped over energy bins
1019 
1020  // Return
1021  return;
1022 }
1023 
1024 
1025 /***********************************************************************//**
1026  * @brief Pad the last delta bins with zero
1027  *
1028  * Zero padding of the last delta bins assures that the Psf goes to zero
1029  * without any step at the last delta value.
1030  ***************************************************************************/
1032 {
1033  // Continue only if there are delta bins
1034  if (m_deltas.size() > 2) {
1035 
1036  // Set first delta bin to value of second bin (capped Psf)
1037  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
1038  int isrc = offset(1, iebin);
1039  int idst = offset(0, iebin);
1040  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1041  m_cube(pixel, idst) = m_cube(pixel, isrc);
1042  }
1043  }
1044 
1045  // Get index of last delta bin
1046  int idelta = m_deltas.size()-1;
1047 
1048  // Pad mean PSF with zeros in the last delta bin
1049  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
1050  int imap = offset(idelta, iebin);
1051  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1052  m_cube(pixel, imap) = 0.0;
1053  }
1054  }
1055 
1056  } // endif: there were bins to pad
1057 
1058  // Return
1059  return;
1060 }
void save(const GFilename &filename, const bool &clobber=false) const
Save PSF cube into FITS file.
std::string print(const GChatter &chatter=NORMAL) const
Print energy container information.
Definition: GEnergies.cpp:764
const GEnergies & energies(void) const
Return energies.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
void clear(void)
Clear energy container.
Definition: GEnergies.cpp:225
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
CTA event list class interface definition.
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 norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
GNodeArray m_deltas
Delta bins (deg) for the PSF cube.
void update(const double &delta, const double &logE) const
Update PSF parameter cache.
#define G_FILL_CUBE
Definition: GCTACubePsf.cpp:45
void write(GFits &file, const std::string &extname=gammalib::extname_energies) const
Write energies into FITS object.
Definition: GEnergies.cpp:725
CTA cube analysis point spread function class definition.
void load(const GFilename &filename)
Load PSF cube from FITS file.
void clear(void)
Clear instance.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:1944
GCTACubePsf & operator=(const GCTACubePsf &cube)
Assignment operator.
GEbounds ebounds(void) const
Get energy boundaries.
void free_members(void)
Delete class members.
void write(GFits &file) const
Write CTA PSF cube into FITS object.
CTA instrument response function class definition.
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
GNodeArray m_deltas_cache
Internal delta bins (rad)
virtual void response(const GResponse &rsp)
Set response function.
const std::string extname_deltas
Definition: GCTACubePsf.hpp:47
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
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
Interface for the circular sky region class.
GFilename m_filename
Filename.
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.
int m_inx1
Index of upper left node.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
const std::vector< int > & shape(void) const
Returns shape of maps.
Definition: GSkyMap.hpp:390
void dir(const GSkyDir &dir)
Set sky direction.
Energy container class.
Definition: GEnergies.hpp:60
double m_wgt4
Weight of lower right node.
void id(const std::string &id)
Set observation identifier.
void write(GFits &fits, const std::string &extname="NODES") const
Write nodes into FITS object.
Definition: GNodeArray.cpp:807
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
void fill(const GObservations &obs, GLog *log=NULL)
Fill PSF cube from observation container.
Information logger interface definition.
Definition: GLog.hpp:62
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
CTA event bin container class.
virtual double livetime(void) const
Return livetime.
double operator()(const GSkyDir &dir, const double &delta, const GEnergy &energy) const
Return point spread function (in units of sr^-1)
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.
void init_members(void)
Initialise class members.
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
double m_wgt3
Weight of upper right node.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
GCTACubePsf * clone(void) const
Clone instance.
int m_inx3
Index of upper right node.
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
void read(const GFits &fits)
Read PSF cube from FITS object.
int m_inx2
Index of lower left node.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void set_eng_axis(void)
Set nodes for a logarithmic (base 10) energy axis.
GChatter
Definition: GTypemaps.hpp:33
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.
CTA observation class interface definition.
Observation container class.
CTA instrument response function class.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
GCTACubePsf(void)
Void constructor.
Definition: GCTACubePsf.cpp:68
const std::string extname_energies
Definition: GEnergies.hpp:44
void set(const GEbounds &ebounds)
Set energies from energy boundaries.
Definition: GEnergies.cpp:416
void set_delta_axis(void)
Set nodes for delta axis in radians.
void name(const std::string &name)
Set observation name.
void set(const GCTAObservation &obs)
Set PSF cube from one CTA observation.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
GEnergies m_energies
Energy values for the PSF cube.
void clear_cube(void)
Clear all pixels in the PSF cube.
Information logger class definition.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
GSkyMap m_cube
PSF cube.
CTA point spread function for cube analysis.
Definition: GCTACubePsf.hpp:62
int m_inx4
Index of lower right node.
std::string print(const GChatter &chatter=NORMAL) const
Print PSF cube information.
CTA event bin container class interface definition.
std::string eventtype(void) const
Return event type string.
bool is_valid(void) const
Checks if RoI is valid.
Definition: GCTARoi.hpp:151
void fill_cube(const GCTAObservation &obs, GSkyMap *exposure=NULL, GLog *log=NULL)
Fill PSF cube for one observation.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
GNodeArray m_elogmeans
Mean log10TeV energy for the PSF cube.
void read(const GFitsTable &table)
Read nodes from FITS table.
Definition: GNodeArray.cpp:770
void set_to_smooth(void)
Pad the last delta bins with zero.
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 int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
virtual ~GCTACubePsf(void)
Destructor.
const GFilename & filename(void) const
Return exposure cube filename.
Observations container class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
int offset(const int &idelta, const int &iebin) const
Return map offset.
Circular sky region class interface definition.
bool m_quadratic_binning
Internal binning is linear.
double m_wgt2
Weight of lower left node.
void copy_members(const GCTACubePsf &cube)
Copy class members.
double m_wgt1
Weight of upper left node.
Mathematical function definitions.
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