GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSPIResponse.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSPIResponse.cpp - INTEGRAL/SPI response class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2020-2022 by Juergen Knoedlseder *
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 GSPIResponse.hpp
23  * @brief INTEGRAL/SPI instrument response function class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <string>
32 #include <typeinfo>
33 #include <algorithm>
34 #include "GException.hpp"
35 #include "GTools.hpp"
36 #include "GEbounds.hpp"
37 #include "GEnergies.hpp"
38 #include "GNodeArray.hpp"
39 #include "GFits.hpp"
40 #include "GFitsBinTable.hpp"
41 #include "GFitsTableShortCol.hpp"
42 #include "GFitsTableDoubleCol.hpp"
43 #include "GSource.hpp"
45 #include "GSPITools.hpp"
46 #include "GSPIResponse.hpp"
47 #include "GSPIObservation.hpp"
48 #include "GSPIEventCube.hpp"
49 #include "GSPIEventBin.hpp"
50 
51 /* __ Method name definitions ____________________________________________ */
52 #define G_IRF "GSPIResponse::irf(GEvent&, GSource&, GObservation&)"
53 #define G_NROI "GSPIResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
54  "GObservation&)"
55 #define G_EBOUNDS "GSPIResponse::ebounds(GEnergy&)"
56 #define G_SET "GSPIResponse::set(const GSPIObservation& obs)"
57 #define G_LOAD_IRF "GSPIResponse::load_irf(GFilename&)"
58 #define G_LOAD_IRFS "GSPIResponse::load_irfs(int&)"
59 
60 /* __ Macros _____________________________________________________________ */
61 
62 /* __ Coding definitions _________________________________________________ */
63 
64 /* __ Debug definitions __________________________________________________ */
65 //#define G_DEBUG_COMPUTE_IRF //!< Debug compute_irf() method
66 
67 /* __ Constants __________________________________________________________ */
68 
69 
70 /*==========================================================================
71  = =
72  = Constructors/destructors =
73  = =
74  ==========================================================================*/
75 
76 /***********************************************************************//**
77  * @brief Void constructor
78  *
79  * Creates an empty INTEGRAL/SPI response.
80  ***************************************************************************/
82 {
83  // Initialise members
84  init_members();
85 
86  // Return
87  return;
88 }
89 
90 
91 /***********************************************************************//**
92  * @brief Response constructor
93  *
94  * Constructs an INTEGRAL/SPI response for an observation using a response
95  * group filename.
96  *
97  * This constructor simply stores the file name of a response group, the
98  * actual loading will be done using the set() method.
99  ***************************************************************************/
101 {
102  // Initialise members
103  init_members();
104 
105  // Store response group filename
106  m_rspname = rspname;
107 
108  // Return
109  return;
110 }
111 
112 
113 /***********************************************************************//**
114  * @brief Copy constructor
115  *
116  * @param[in] rsp INTEGRAL/SPI response.
117  **************************************************************************/
119 {
120  // Initialise members
121  init_members();
122 
123  // Copy members
124  copy_members(rsp);
125 
126  // Return
127  return;
128 }
129 
130 
131 /***********************************************************************//**
132  * @brief Destructor
133  ***************************************************************************/
135 {
136  // Free members
137  free_members();
138 
139  // Return
140  return;
141 }
142 
143 
144 /*==========================================================================
145  = =
146  = Operators =
147  = =
148  ==========================================================================*/
149 
150 /***********************************************************************//**
151  * @brief Assignment operator
152  *
153  * @param[in] rsp INTEGRAL/SPI response.
154  * @return INTEGRAL/SPI response.
155  *
156  * Assign INTEGRAL/SPI response to this object. The assignment performs
157  * a deep copy of all information, hence the original object from which the
158  * assignment was made can be destroyed after this operation without any loss
159  * of information.
160  ***************************************************************************/
162 {
163  // Execute only if object is not identical
164  if (this != &rsp) {
165 
166  // Copy base class members
167  this->GResponse::operator=(rsp);
168 
169  // Free members
170  free_members();
171 
172  // Initialise members
173  init_members();
174 
175  // Copy members
176  copy_members(rsp);
177 
178  } // endif: object was not identical
179 
180  // Return this object
181  return *this;
182 }
183 
184 
185 /*==========================================================================
186  = =
187  = Public methods =
188  = =
189  ==========================================================================*/
190 
191 /***********************************************************************//**
192  * @brief Clear instance
193  *
194  * Clears INTEGRAL/SPI response by resetting all class members to an initial
195  * state. Any information that was present before will be lost.
196  ***************************************************************************/
198 {
199  // Free class members (base and derived classes, derived class first)
200  free_members();
201  this->GResponse::free_members();
202 
203  // Initialise members
204  this->GResponse::init_members();
205  init_members();
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Clone instance
214  *
215  * @return Pointer to deep copy of INTEGRAL/SPI response.
216  ***************************************************************************/
218 {
219  return new GSPIResponse(*this);
220 }
221 
222 
223 /***********************************************************************//**
224  * @brief Return value of INTEGRAL/SPI instrument response for a photon
225  *
226  * @param[in] event Observed event.
227  * @param[in] photon Incident photon.
228  * @param[in] obs Observation.
229  * @return Instrument response \f$(cm^2 sr^{-1})\f$
230  *
231  * @exception GException::invalid_argument
232  * Observation is not a INTEGRAL/SPI observation.
233  *
234  * Returns the instrument response function for a given observed photon
235  * direction as function of the assumed true photon direction. The result
236  * is given by
237  *
238  * \f[
239  * R(p'|p) =
240  * \f]
241  *
242  * @todo Write down formula
243  * @todo Describe in detail how the response is computed.
244  ***************************************************************************/
245 double GSPIResponse::irf(const GEvent& event,
246  const GPhoton& photon,
247  const GObservation& obs) const
248 {
249  // Initialise IRF
250  double irf = 0.0;
251 
252  // Extract INTEGRAL/SPI observation
253  const GSPIObservation* spi_obs = dynamic_cast<const GSPIObservation*>(&obs);
254  if (spi_obs == NULL) {
255  std::string cls = std::string(typeid(&obs).name());
256  std::string msg = "Observation of type \""+cls+"\" is not an "
257  "INTEGRAL/SPI observation. Please specify an "
258  "INTEGRAL/SPI observation as argument.";
260  }
261 
262  // Extract INTEGRAL/SPI event cube
263  const GSPIEventCube* cube = dynamic_cast<const GSPIEventCube*>(spi_obs->events());
264  if (cube == NULL) {
265  std::string msg = "INTEGRAL/SPI observation does not contain a valid "
266  "event cube. Please specify an observation with an "
267  "event cube as argument.";
269  }
270 
271  // Extract INTEGRAL/SPI event bin
272  const GSPIEventBin* bin = dynamic_cast<const GSPIEventBin*>(&event);
273  if (bin == NULL) {
274  std::string cls = std::string(typeid(&event).name());
275  std::string msg = "Event of type \""+cls+"\" is not an INTEGRAL/SPI "
276  "event. Please specify an INTEGRAL/SPI event as "
277  "argument.";
279  }
280 
281  // Continue only if livetime for event is positive
282  if (bin->livetime() > 0.0) {
283 
284  // Get energy bins of event and photon. Since only the photo peak
285  // response is supported so far, both indices need to be identical.
286  // Continue only if this is the case.
287  if (cube->ebounds().index(photon.energy()) == bin->iebin()) {
288 
289  // Get IRF value for photo peak
290  irf = irf_value(photon.dir(), *bin, 0);
291 
292  // Compile option: Check for NaN/Inf
293  #if defined(G_NAN_CHECK)
295  std::cout << "*** ERROR: GSPIResponse::irf:";
296  std::cout << " NaN/Inf encountered";
297  std::cout << " (irf=" << irf;
298  std::cout << ")";
299  std::cout << std::endl;
300  }
301  #endif
302 
303  } // endif: photon energy in same energy bin as event energy
304 
305  } // endif: livetime of event was positive
306 
307  // Return IRF value
308  return irf;
309 }
310 
311 
312 /***********************************************************************//**
313  * @brief Return value of INTEGRAL/SPI instrument response for sky direction
314  * and event bin
315  *
316  * @param[in] srcDir Sky direction.
317  * @param[in] bin INTEGRAL/SPI event bin.
318  * @param[in] ireg IRF region (0: photo peak).
319  * @return Instrument response \f$(cm^2 sr^{-1})\f$
320  *
321  * Returns the instrument response function for a given sky direction and
322  * event bin. The value of the IRF is bilinearly interpolated from the
323  * pre-computed IRFs cube that is stored in the class.
324  ***************************************************************************/
325 double GSPIResponse::irf_value(const GSkyDir& srcDir,
326  const GSPIEventBin& bin,
327  const int& ireg) const
328 {
329  // Initialise IRF value
330  double irf = 0.0;
331 
332  // Convert sky direction to zenith angle. Continue only is zenith angle
333  // is below maximum zenith angle
334  double zenith = this->zenith(bin.ipt(), srcDir);
335  if (zenith < m_max_zenith) {
336 
337  // Convert sky direction to azimuth angle
338  double azimuth = this->azimuth(bin.ipt(), srcDir);
339 
340  // Compute pixel
341  double xpix = (zenith * std::cos(azimuth) - m_wcs_xmin) / m_wcs_xbin;
342  double ypix = (zenith * std::sin(azimuth) - m_wcs_ymin) / m_wcs_ybin;
343 
344  // Continue only if pixel is within IRF
345  if (xpix > 0.0 && xpix < m_wcs_xpix_max &&
346  ypix > 0.0 && ypix < m_wcs_ypix_max) {
347 
348  // Get number of pixels in X direction
349  int nx = m_irfs.nx();
350  int ndet = m_irfs.shape()[0];
351  int nreg = m_irfs.shape()[1];
352 
353  // Get IRF detector for event
354  int idet = irf_detid(bin.dir().detid());
355  int ieng = bin.iebin();
356 
357  // Get map
358  int map = idet + (ireg + ieng * nreg) * ndet;
359 
360  // Get 4 nearest neighbours
361  int ix_left = int(xpix);
362  int ix_right = ix_left + 1;
363  int iy_top = int(ypix);
364  int iy_bottom = iy_top + 1;
365 
366  // Get weighting factors
367  double wgt_right = xpix - double(ix_left);
368  double wgt_left = 1.0 - wgt_right;
369  double wgt_bottom = ypix - double(iy_top);
370  double wgt_top = 1.0 - wgt_bottom;
371 
372  // Get indices of 4 nearest neighbours
373  int inx_1 = ix_left + iy_top * nx;
374  int inx_2 = ix_right + iy_top * nx;
375  int inx_3 = ix_left + iy_bottom * nx;
376  int inx_4 = ix_right + iy_bottom * nx;
377 
378  // Get weights of 4 nearest neighbours
379  double wgt_1 = wgt_left * wgt_top;
380  double wgt_2 = wgt_right * wgt_top;
381  double wgt_3 = wgt_left * wgt_bottom;
382  double wgt_4 = wgt_right * wgt_bottom;
383 
384  // Compute IRF
385  irf = m_irfs(inx_1, map) * wgt_1;
386  irf += m_irfs(inx_2, map) * wgt_2;
387  irf += m_irfs(inx_3, map) * wgt_3;
388  irf += m_irfs(inx_4, map) * wgt_4;
389 
390  // Make sure that IRF does not get negative
391  if (irf < 0.0) {
392  irf = 0.0;
393  }
394 
395  } // endif: zenith angle was valid
396 
397  } // endif: pixel was within IRF
398 
399  // Return IRF value
400  return irf;
401 }
402 
403 
404 /***********************************************************************//**
405  * @brief Return integral of event probability for a given sky model over ROI
406  *
407  * @param[in] model Sky model.
408  * @param[in] obsEng Observed photon energy.
409  * @param[in] obsTime Observed photon arrival time.
410  * @param[in] obs Observation.
411  * @return 0.0
412  *
413  * @exception GException::feature_not_implemented
414  * Method is not implemented.
415  ***************************************************************************/
416 double GSPIResponse::nroi(const GModelSky& model,
417  const GEnergy& obsEng,
418  const GTime& obsTime,
419  const GObservation& obs) const
420 {
421  // Method is not implemented
422  std::string msg = "Spatial integration of sky model over the data space "
423  "is not implemented.";
425 
426  // Return Npred
427  return (0.0);
428 }
429 
430 
431 /***********************************************************************//**
432  * @brief Return true energy boundaries for a specific observed energy
433  *
434  * @param[in] obsEnergy Observed Energy.
435  * @return True energy boundaries for given observed energy.
436  *
437  * @exception GException::feature_not_implemented
438  * Method is not implemented.
439  *
440  * @todo Implement this method if you need energy dispersion.
441  ***************************************************************************/
442 GEbounds GSPIResponse::ebounds(const GEnergy& obsEnergy) const
443 {
444  // Initialise an empty boundary object
446 
447  // Throw an exception
448  std::string msg = "Energy dispersion not implemented.";
450 
451  // Return energy boundaries
452  return ebounds;
453 }
454 
455 
456 /***********************************************************************//**
457  * @brief Set response for a specific observation
458  *
459  * @param[in] obs INTEGRAL/SPI observation.
460  * @param[in] energy Line energy
461  *
462  * Set response for a specific INTEGRAL/SPI observation.
463  *
464  * If the @p energy argument is set to a positive value, the IRF will be
465  * computed for the specified line energy.
466  ***************************************************************************/
467 void GSPIResponse::set(const GSPIObservation& obs, const GEnergy& energy)
468 {
469  // Reset response information
470  m_detids.clear();
471  m_energies.clear();
472  m_ebounds.clear();
473  m_irfs.clear();
474  m_has_wcs = false;
475 
476  // Continue only if the observation contains an event cube
477  const GSPIEventCube* cube = dynamic_cast<const GSPIEventCube*>(obs.events());
478  if (cube != NULL) {
479 
480  // Set requested detector identifiers
481  set_detids(cube);
482 
483  // Set coordination transformation cache
484  set_cache(cube);
485 
486  // Load IRFs for photo peak region (region 0)
487  load_irfs(0);
488 
489  // Get dimension of pre-computed IRF
490  int npix = m_irfs.npix();
491  int ndet = m_irfs.shape()[0];
492  int nreg = m_irfs.shape()[1];
493  int neng = cube->naxis(2);
494 
495  // Initialise pre-computed IRF
496  GSkyMap irfs = m_irfs;
497  irfs.nmaps(ndet * nreg * neng);
498  irfs.shape(ndet, nreg, neng);
499  irfs = 0.0;
500 
501  // Pre-compute IRFs for all energy bins
502  for (int ieng = 0; ieng < neng; ++ieng) {
503 
504  // Get energy boundaries. If line argument is true then use the
505  double emin = cube->ebounds().emin(ieng).keV();
506  double emax = cube->ebounds().emax(ieng).keV();
507 
508  // If line energy is specified then compute a line response for
509  // the energy bin that overlaps with the line energy
510  m_energy_keV = energy.keV();
511  if (m_energy_keV > 0.0) {
512  if (m_energy_keV < emin || m_energy_keV > emax) {
513  continue;
514  }
515  emin = m_energy_keV;
516  emax = m_energy_keV;
517  }
518 
519  // Pre-compute IRF for this energy bin
520  GSkyMap irf = compute_irf(emin, emax);
521 
522  // Copy over IRF
523  for (int idet = 0; idet < ndet; ++idet) {
524  for (int ireg = 0; ireg < nreg; ++ireg) {
525  int map_irf = idet + ireg * ndet;
526  int map_irfs = idet + (ireg + ieng * nreg) * ndet;
527  for (int i = 0; i < npix; ++i) {
528  irfs(i, map_irfs) = irf(i, map_irf);
529  }
530  }
531  }
532 
533  // Append energy boundary
534  m_ebounds.append(GEnergy(emin, "keV"), GEnergy(emax, "keV"));
535 
536  } // endfor: looped over all energy bins
537 
538  // Replace IRFs
539  m_irfs = irfs;
540 
541  // Clear energies
542  m_energies.clear();
543 
544  } // endif: observation contained an event cube
545 
546  // Return
547  return;
548 }
549 
550 
551 /***********************************************************************//**
552  * @brief Read SPI response from FITS object
553  *
554  * @param[in] fits FITS object.
555  *
556  * Reads the SPI response from FITS object.
557  ***************************************************************************/
558 void GSPIResponse::read(const GFits& fits)
559 {
560  // Get IRF image
561  const GFitsImage* image = fits.image(0);
562 
563  // Read IRFs
564  m_irfs.read(*image);
565 
566  // Set WCS image limits
567  set_wcs(image);
568 
569  // Read detector identifiers
570  read_detids(fits);
571 
572  // Read energies
573  read_energies(fits);
574 
575  // Return
576  return;
577 }
578 
579 
580 /***********************************************************************//**
581  * @brief Write SPI response into FITS object
582  *
583  * @param[in,out] fits FITS object.
584  *
585  * Writes the SPI response into FITS object.
586  ***************************************************************************/
587 void GSPIResponse::write(GFits& fits) const
588 {
589  // Write IRFs as sky map
590  m_irfs.write(fits);
591 
592  // Write detector identifiers
593  write_detids(fits);
594 
595  // Write energies
596  write_energies(fits);
597 
598  // Return
599  return;
600 }
601 
602 
603 /***********************************************************************//**
604  * @brief Load SPI response from file
605  *
606  * @param[in] filename Response file name.
607  *
608  * Loads SPI response from response file.
609  ***************************************************************************/
610 void GSPIResponse::load(const GFilename& filename)
611 {
612  // Open FITS file
613  GFits fits(filename);
614 
615  // Read response
616  read(fits);
617 
618  // Close FITS file
619  fits.close();
620 
621  // Return
622  return;
623 }
624 
625 
626 /***********************************************************************//**
627  * @brief Save SPI response into file
628  *
629  * @param[in] filename Response file name.
630  * @param[in] clobber Overwrite existing FITS file?
631  *
632  * Saves SPI response into response file.
633  ***************************************************************************/
634 void GSPIResponse::save(const GFilename& filename, const bool& clobber) const
635 {
636  // Creat FITS file
637  GFits fits;
638 
639  // Write response
640  write(fits);
641 
642  // Save FITS file
643  fits.saveto(filename, clobber);
644 
645  // Close FITS file
646  fits.close();
647 
648  // Return
649  return;
650 }
651 
652 
653 /***********************************************************************//**
654  * @brief Print INTEGRAL/SPI response information
655  *
656  * @param[in] chatter Chattiness.
657  * @return String containing INTEGRAL/SPI response information.
658  ***************************************************************************/
659 std::string GSPIResponse::print(const GChatter& chatter) const
660 {
661  // Initialise result string
662  std::string result;
663 
664  // Continue only if chatter is not silent
665  if (chatter != SILENT) {
666 
667  // Extract information
668  int nx = m_irfs.nx();
669  int ny = m_irfs.ny();
670  int ndet = (m_irfs.shape().size() > 0) ? m_irfs.shape()[0] : 0;
671  int nreg = (m_irfs.shape().size() > 1) ? m_irfs.shape()[1] : 0;
672  int neng = (m_irfs.shape().size() > 2) ? m_irfs.shape()[2] : 0;
673 
674  // Append header
675  result.append("=== GSPIResponse ===");
676 
677  // Append information
678  result.append("\n"+gammalib::parformat("Response group filename"));
679  result.append(m_rspname.url());
680  result.append("\n"+gammalib::parformat("Response map pixels"));
681  result.append(gammalib::str(nx)+" * "+gammalib::str(ny));
682  result.append("\n"+gammalib::parformat("X axis range"));
683  result.append("["+gammalib::str(m_wcs_xmin * gammalib::rad2deg));
684  result.append(","+gammalib::str(m_wcs_xmax * gammalib::rad2deg));
685  result.append("] deg");
686  result.append("\n"+gammalib::parformat("Y axis range"));
687  result.append("["+gammalib::str(m_wcs_ymin * gammalib::rad2deg));
688  result.append(","+gammalib::str(m_wcs_ymax * gammalib::rad2deg));
689  result.append("] deg");
690  result.append("\n"+gammalib::parformat("Maximum zenith angle"));
691  result.append(gammalib::str(m_max_zenith * gammalib::rad2deg)+" deg");
692  result.append("\n"+gammalib::parformat("Number of detectors"));
693  result.append(gammalib::str(ndet));
694  result.append("\n"+gammalib::parformat("Number of regions"));
695  result.append(gammalib::str(nreg));
696  result.append("\n"+gammalib::parformat("Number of energies"));
697  result.append(gammalib::str(neng));
698  if (m_energy_keV > 0.0) {
699  result.append("\n"+gammalib::parformat("Line IRF energy"));
700  result.append(gammalib::str(m_energy_keV)+" keV");
701  }
702  else {
703  result.append("\n"+gammalib::parformat("Continuum IRF gamma"));
704  result.append(gammalib::str(m_gamma));
705  result.append("\n"+gammalib::parformat("Continnum IRF log(E) step"));
706  result.append(gammalib::str(m_dlogE));
707  }
708 
709  } // endif: chatter was not silent
710 
711  // Return result
712  return result;
713 }
714 
715 
716 /*==========================================================================
717  = =
718  = Private methods =
719  = =
720  ==========================================================================*/
721 
722 /***********************************************************************//**
723  * @brief Initialise class members
724  ***************************************************************************/
726 {
727  // Initialise members
728  m_rspname.clear();
729  m_detids.clear();
730  m_energies.clear();
731  m_ebounds.clear();
732  m_irfs.clear();
733  m_energy_keV = 0.0;
734  m_dlogE = 0.03;
735  m_gamma = 2.0;
736 
737  // Initialise cache
738  m_spix.clear();
739  m_posang.clear();
740  m_has_wcs = false;
741  m_wcs_xmin = 0.0;
742  m_wcs_ymin = 0.0;
743  m_wcs_xmax = 0.0;
744  m_wcs_ymax = 0.0;
745  m_wcs_xbin = 0.0;
746  m_wcs_ybin = 0.0;
747  m_wcs_xpix_max = 0.0;
748  m_wcs_ypix_max = 0.0;
750 
751  // Return
752  return;
753 }
754 
755 
756 /***********************************************************************//**
757  * @brief Copy class members
758  *
759  * @param[in] rsp INTEGRAL/SPI response function.
760  ***************************************************************************/
762 {
763  // Copy members
764  m_rspname = rsp.m_rspname;
765  m_detids = rsp.m_detids;
766  m_energies = rsp.m_energies;
767  m_ebounds = rsp.m_ebounds;
768  m_irfs = rsp.m_irfs;
770  m_dlogE = rsp.m_dlogE;
771  m_gamma = rsp.m_gamma;
772 
773  // Copy cache
774  m_spix = rsp.m_spix;
775  m_posang = rsp.m_posang;
776  m_has_wcs = rsp.m_has_wcs;
777  m_wcs_xmin = rsp.m_wcs_xmin;
778  m_wcs_ymin = rsp.m_wcs_ymin;
779  m_wcs_xmax = rsp.m_wcs_xmax;
780  m_wcs_ymax = rsp.m_wcs_ymax;
781  m_wcs_xbin = rsp.m_wcs_xbin;
782  m_wcs_ybin = rsp.m_wcs_ybin;
786 
787  // Return
788  return;
789 }
790 
791 
792 /***********************************************************************//**
793  * @brief Delete class members
794  ***************************************************************************/
796 {
797  // Return
798  return;
799 }
800 
801 
802 /***********************************************************************//**
803  * @brief Read detector identifiers from FITS object
804  *
805  * @param[in] fits FITS object.
806  *
807  * Read the detector identifiers from a FITS file into the m_detids member.
808  ***************************************************************************/
810 {
811  // Clear detids vector
812  m_detids.clear();
813 
814  // Get DETID table extension
815  const GFitsTable& table = *fits.table("DETIDS");
816  const GFitsTableCol* col_detid = table["DET_ID"];
817 
818  // Get number of detector identifiers
819  int num = table.nrows();
820 
821  // Reserve space for detector identifiers
822  m_detids.reserve(num);
823 
824  // Extract values
825  for (int i = 0; i < num; ++i) {
826  m_detids.push_back(col_detid->integer(i));
827  }
828 
829  // Return
830  return;
831 }
832 
833 
834 /***********************************************************************//**
835  * @brief Read energies from FITS object
836  *
837  * @param[in] fits FITS object.
838  *
839  * Read the energies from a FITS file. If the FITS file contains energy
840  * boundaries then read them, otherwise read the energy vector.
841  ***************************************************************************/
843 {
844  // Clear energies and energy boundary vectors
845  m_energies.clear();
846  m_ebounds.clear();
847 
848  // If FITS file contains energy boundaries then the IRF was precomputed
849  // and hence we read the energy boundaries from the FITS file
851 
852  // Get reference to FITS table
853  const GFitsTable& table = *fits.table(gammalib::extname_ebounds);
854 
855  // Read energy boundaries
856  m_ebounds.read(table);
857 
858  // Set response attributes
859  m_energy_keV = (table.has_card("ENERGY")) ? table.real("ENERGY") : 0.0;
860  m_dlogE = (table.has_card("DLOGE")) ? table.real("DLOGE") : 0.03;
861  m_gamma = (table.has_card("GAMMA")) ? table.real("GAMMA") : 2.0;
862 
863  } // endif: energy boundaries detected
864 
865  // ... otherwise we have a response group and we read the energies of
866  // the response group
867  else {
868 
869  // Get energies table extension
870  const GFitsTable& table = *fits.table(gammalib::extname_energies);
871  const GFitsTableCol* col_energy = table["ENERGY"];
872 
873  // Get number of energies
874  int num = table.nrows();
875 
876  // Reserve space for energies
877  m_energies.reserve(num);
878 
879  // Get energy unit
880  std::string unit = "keV";
881  if (table.has_card("TUNIT1")) {
882  unit = table.string("TUNIT1");
883  }
884 
885  // Extract values
886  for (int i = 0; i < num; ++i) {
887 
888  // Get energy
889  GEnergy energy = GEnergy(col_energy->real(i), unit);
890 
891  // Store as keV
892  m_energies.append(energy.keV());
893 
894  } // endfor: looped over all values
895 
896  } // endelse: energies extension was required
897 
898  // Return
899  return;
900 }
901 
902 
903 /***********************************************************************//**
904  * @brief Write detector identifiers into FITS object
905  *
906  * @param[in,out] fits FITS object.
907  *
908  * Writes the detector identifiers stored into the m_detids member into a
909  * FITS object.
910  ***************************************************************************/
912 {
913  // Get number of detector identifiers
914  int num = m_detids.size();
915 
916  // Create columns
917  GFitsTableShortCol col_detid("DET_ID", num);
918 
919  // Fill energy column in units of keV
920  for (int i = 0; i < num; ++i) {
921  col_detid(i) = m_detids[i];
922  }
923 
924  // Create energies table
925  GFitsBinTable table(num);
926  table.append(col_detid);
927  table.extname("DETIDS");
928 
929  // Append detector identifiers to FITS file
930  fits.append(table);
931 
932  // Return
933  return;
934 }
935 
936 
937 /***********************************************************************//**
938  * @brief Write energies into FITS object
939  *
940  * @param[in,out] fits FITS object.
941  *
942  * Writes the energy nodes stored into the m_energies member into a FITS
943  * object. The energies are written in units of keV.
944  ***************************************************************************/
946 {
947  // If there are energy boundaries then the IRF was precomputed and
948  // hence we store the energy boundaries in the FITS file
949  if (!m_ebounds.is_empty()) {
950 
951  // Write energy boundaries
952  m_ebounds.write(fits);
953 
954  // Recover the FITS table to write some keywords
956 
957  // Write keywords
958  table.card("ENERGY", m_energy_keV, "[keV] IRF line energy (0 for continuum IRF)");
959  table.card("DLOGE", m_dlogE, "Logarithmic step size for continuum IRF");
960  table.card("GAMMA", m_gamma, "Power-law index for continuum IRF");
961 
962  } // endif: IRF was precomputed
963 
964  // ... otherwise we write the energies
965  else {
966 
967  // Get number of energies
968  int num = m_energies.size();
969 
970  // Create columns
971  GFitsTableDoubleCol col_energy("ENERGY", num);
972 
973  // Fill energy column in units of keV
974  for (int i = 0; i < num; ++i) {
975  col_energy(i) = m_energies[i];
976  }
977  col_energy.unit("keV");
978 
979  // Create energies table
980  GFitsBinTable table(num);
981  table.append(col_energy);
983 
984  // Append energies table to FITS file
985  fits.append(table);
986 
987  } // endelse: energies written
988 
989  // Return
990  return;
991 }
992 
993 
994 /***********************************************************************//**
995  * @brief Load IRF as sky map
996  *
997  * @param[in] irfname IRF file name.
998  *
999  * Loads an IRF FITS file as a sky map. The sky map is returned in ARC
1000  * projection that is a zenith equidistant projection, which is the
1001  * projection that is natively used to store the IRFs.
1002  ***************************************************************************/
1004 {
1005  // Open IRF FITS file
1006  GFits fits(irfname);
1007 
1008  // Access IRF FITS image
1009  GFitsImage* image = fits.image("SPI.-IRF.-RSP");
1010 
1011  // Get image attributes
1012  int naxis1 = image->integer("NAXIS1");
1013  int naxis2 = image->integer("NAXIS2");
1014  int naxis3 = image->integer("NAXIS3");
1015  int naxis4 = image->integer("NAXIS4");
1016  double crval2 = image->real("CRVAL2");
1017  double crval3 = image->real("CRVAL3");
1018  double crpix2 = image->real("CRPIX2");
1019  double crpix3 = image->real("CRPIX3");
1020  double cdelt2 = image->real("CDELT2");
1021  double cdelt3 = image->real("CDELT3");
1022 
1023  // Derive image limits. Limits are taken at the pixel centres since
1024  // we want to use them for bilinear interpolation. This means that
1025  // we will throw away half a pixel at the edge of the IRFs.
1026  double wcs_xmin = (crval2 - (crpix2-1.0) * cdelt2) * gammalib::deg2rad;
1027  double wcs_ymin = (crval3 - (crpix3-1.0) * cdelt3) * gammalib::deg2rad;
1028  double wcs_xbin = cdelt2 * gammalib::deg2rad;
1029  double wcs_ybin = cdelt3 * gammalib::deg2rad;
1030  double wcs_xmax = wcs_xmin + double(naxis2-1) * wcs_xbin;
1031  double wcs_ymax = wcs_ymin + double(naxis3-1) * wcs_ybin;
1032  double wcs_xpix_max = double(naxis2-1);
1033  double wcs_ypix_max = double(naxis3-1);
1034 
1035  // If no image limits exists so far then store them for fast IRF access
1036  // that does not depend on the actual IRF projection (we just use here
1037  // the sky map as a convient container)
1038  if (!m_has_wcs) {
1039  m_has_wcs = true;
1040  m_wcs_xmin = wcs_xmin;
1041  m_wcs_ymin = wcs_ymin;
1042  m_wcs_xbin = wcs_xbin;
1043  m_wcs_ybin = wcs_ybin;
1044  m_wcs_xmax = wcs_xmax;
1045  m_wcs_ymax = wcs_ymax;
1046  m_wcs_xpix_max = wcs_xpix_max;
1047  m_wcs_ypix_max = wcs_ypix_max;
1048  }
1049 
1050  // ... otherwise check if the limits are consistent
1051  else {
1052  if ((std::abs(wcs_xmin - m_wcs_xmin) > 1.0e-6) ||
1053  (std::abs(wcs_ymin - m_wcs_ymin) > 1.0e-6) ||
1054  (std::abs(wcs_xbin - m_wcs_xbin) > 1.0e-6) ||
1055  (std::abs(wcs_ybin - m_wcs_ybin) > 1.0e-6) ||
1056  (std::abs(wcs_xmax - m_wcs_xmax) > 1.0e-6) ||
1057  (std::abs(wcs_ymax - m_wcs_ymax) > 1.0e-6) ||
1058  (std::abs(wcs_xpix_max - m_wcs_xpix_max) > 1.0e-6) ||
1059  (std::abs(wcs_ypix_max - m_wcs_ypix_max) > 1.0e-6)) {
1060  std::string msg = "Inconsistent IRFs encountered in file \""+
1061  irfname.url()+"\". Please specify a response "
1062  "group where all IRFs have the same definition.";
1064  }
1065  }
1066 
1067  // Set maximum zenith angle
1070 
1071  // Compute sky map attributes
1072  int nmap = naxis1 * naxis4;
1073 
1074  // Initialise IRF in celestial coordinates. We use an ARC projection as
1075  // this is the native projection of the IRF files.
1076  GSkyMap irf("ARC", "CEL", crval2, crval3, cdelt2, cdelt3, naxis2, naxis3, nmap);
1077 
1078  // Fill sky map
1079  for (int ix = 0; ix < naxis2; ++ix) {
1080  for (int iy = 0; iy < naxis3; ++iy) {
1081  for (int idet = 0; idet < naxis1; ++idet) {
1082  for (int ireg = 0; ireg < naxis4; ++ireg) {
1083  int index = ix + iy * naxis2;
1084  int map = idet + ireg * naxis1;
1085  irf(index,map) = image->pixel(idet, ix, iy, ireg);
1086  }
1087  }
1088  }
1089  }
1090 
1091  // Set sky map shape
1092  irf.shape(naxis1, naxis4);
1093 
1094  // Close FITS file
1095  fits.close();
1096 
1097  // Return IRF
1098  return irf;
1099 }
1100 
1101 
1102 /***********************************************************************//**
1103  * @brief Load Instrument Response Functions
1104  *
1105  * @param[in] region IRF region (-1 = load all regions)
1106  *
1107  * The method requires that the required detector identifiers were previously
1108  * setup using the set_detids() method.
1109  ***************************************************************************/
1110 void GSPIResponse::load_irfs(const int& region)
1111 {
1112  // Initialise results
1113  m_energies.clear();
1114  m_irfs.clear();
1115 
1116  // Determine number of requested detectors. Continue only if there are
1117  // required detectors.
1118  int num_det = m_detids.size();
1119  if (num_det > 0) {
1120 
1121  // Open response group
1122  GFits rsp(m_rspname);
1123 
1124  // Get response grouping table
1125  const GFitsTable* grp = gammalib::spi_hdu(rsp, "SPI.-IRF.-RSP-IDX");
1126  if (grp == NULL) {
1127  std::string msg = "No response grouping table found in file \""+
1128  m_rspname.url()+"\". Please specify a valid "
1129  "response grouping file.";
1131  }
1132 
1133  // Determine number of response members
1134  int num_irfs = gammalib::spi_num_hdus(rsp, "SPI.-IRF.-RSP");
1135 
1136  // Setup node array for response energies
1137  for (int i_irf = 0; i_irf < num_irfs; ++i_irf) {
1138  m_energies.append((*grp)["ENERGY"]->real(i_irf)); // in keV
1139  }
1140 
1141  // Loop over all IRFs
1142  for (int i_irf = 0; i_irf < num_irfs; ++i_irf) {
1143 
1144  // Get IRF filename
1145  std::string irfname = m_rspname.path() +
1146  (*grp)["MEMBER_LOCATION"]->string(i_irf);
1147 
1148  // Load IRF
1149  GSkyMap irf = load_irf(irfname);
1150 
1151  // Determine number of requested regions
1152  int num_regions = (region == -1) ? irf.shape()[1] : 1;
1153 
1154  // If this is the first IRF then allocate IRFs
1155  if (i_irf == 0) {
1156 
1157  // Set IRFs from IRF
1158  m_irfs = irf;
1159 
1160  // Allocate sufficient maps and reshape maps
1161  m_irfs.nmaps(num_det * num_regions * num_irfs);
1162  m_irfs.shape(num_det, num_regions, num_irfs);
1163 
1164  // Initialise maps
1165  m_irfs = 0.0;
1166 
1167  } // endif: this was the first IRF
1168 
1169  // Extract relevant part or IRF
1170  int nx = irf.nx();
1171  int ny = irf.ny();
1172  int ndet = irf.shape()[0];
1173 
1174  // Loop over requested regions
1175  for (int i_region = 0; i_region < num_regions; ++i_region) {
1176 
1177  // Loop over requested detectors
1178  for (int i_det = 0; i_det < num_det; ++i_det) {
1179 
1180  // Set map indices in source IRF and target IRFs
1181  int map_irf = m_detids[i_det] + i_region * ndet;
1182  int map_irfs = i_det + (i_region + i_irf * num_regions) * num_det;
1183 
1184  // Copy IRF pixels
1185  for (int i = 0, iy = 0; iy < ny; ++iy) {
1186  for (int ix = 0; ix < nx; ++ix, ++i) {
1187  m_irfs(i, map_irfs) = irf(i, map_irf);
1188  }
1189  }
1190 
1191  } // endfor: looped over requested detectors
1192 
1193  } // endfor: looped over requested regions
1194 
1195  } // endfor: looped over all IRFs
1196 
1197  // Close FITS file
1198  rsp.close();
1199 
1200  } // endif: there were requested detectors
1201 
1202  // Return
1203  return;
1204 }
1205 
1206 
1207 /***********************************************************************//**
1208  * @brief Compute as sky map
1209  *
1210  * @param[in] emin Minimum energy (keV).
1211  * @param[in] emax Maximum energy (keV).
1212  *
1213  * Compute IRF for an energy band specified by an energy interval. If the
1214  * width of the energy interval is zero a line IRF will be computed.
1215  ***************************************************************************/
1216 GSkyMap GSPIResponse::compute_irf(const double& emin, const double& emax) const
1217 {
1218  // Get IRF dimension
1219  int npix = m_irfs.npix();
1220  int ndet = m_irfs.shape()[0];
1221  int nreg = m_irfs.shape()[1];
1222 
1223  // Debug: print dimensions
1224  #if defined(G_DEBUG_COMPUTE_IRF)
1225  std::cout << "GSPIResponse::compute_irf:" << std::endl;
1226  std::cout << "- emin = " << emin << " keV" << std::endl;
1227  std::cout << "- emax = " << emax << " keV" << std::endl;
1228  std::cout << "- npix = " << npix << std::endl;
1229  std::cout << "- ndet = " << ndet << std::endl;
1230  std::cout << "- nreg = " << nreg << std::endl;
1231  #endif
1232 
1233  // Initialise IRF
1234  GSkyMap irf = m_irfs;
1235  irf.nmaps(ndet * nreg);
1236  irf.shape(ndet, nreg);
1237  irf = 0.0;
1238 
1239  // Case A: Integrate response over energy band
1240  if (emax > emin) {
1241 
1242  // Get logarithmic energy limits
1243  double log_emin = std::log10(emin);
1244  double log_emax = std::log10(emax);
1245  double log_ewidth = log_emax - log_emin;
1246 
1247  // Set energy weighting factors
1248  double beta = 1.0 - m_gamma;
1249  double wgt0 = 1.0 / irf_weight(beta, log_emin, log_emax);
1250 
1251  // Determine number of integration intervals and logarithmic step
1252  // size
1253  int num_int = (m_dlogE > DBL_MIN) ? int(log_ewidth / m_dlogE + 0.5) : 1;
1254  if (num_int < 1) {
1255  num_int = 1;
1256  }
1257  double log_estep = log_ewidth / num_int;
1258 
1259  // Debug: print weighting and internal binning
1260  #if defined(G_DEBUG_COMPUTE_IRF)
1261  std::cout << "- log_emin = " << log_emin << std::endl;
1262  std::cout << "- log_emax = " << log_emax << std::endl;
1263  std::cout << "- log_ewidth = " << log_ewidth << std::endl;
1264  std::cout << "- beta = " << beta << std::endl;
1265  std::cout << "- wgt0 = " << wgt0 << std::endl;
1266  std::cout << "- num_int = " << num_int << std::endl;
1267  std::cout << "- log_estep = " << log_estep << std::endl;
1268  double wgt_check = 0.0;
1269  #endif
1270 
1271  // Loop over integration intervals
1272  for (int i_int = 0; i_int < num_int; ++i_int) {
1273 
1274  // Determine integration interval width and midpoint
1275  double log_emin_bin = log_emin + i_int * log_estep;
1276  double log_emax_bin = log_emin_bin + log_estep;
1277  double energy = std::exp(0.5 * gammalib::ln10 *
1278  (log_emin_bin + log_emax_bin));
1279 
1280  // Get response weighting factor
1281  double wgt = irf_weight(beta, log_emin_bin, log_emax_bin) * wgt0;
1282 
1283  // Get IRFs bracketing the mean energy and the IRF weights
1284  m_energies.set_value(energy);
1285  int irf_low = m_energies.inx_left();
1286  int irf_up = m_energies.inx_right();
1287  double wgt_low = m_energies.wgt_left() * wgt;
1288  double wgt_up = m_energies.wgt_right() * wgt;
1289 
1290  // Debug: add weight for final check
1291  #if defined(G_DEBUG_COMPUTE_IRF)
1292  wgt_check += wgt_low + wgt_up;
1293  #endif
1294 
1295  // Add contribution to all IRF pixels
1296  for (int idet = 0; idet < ndet; ++idet) {
1297  for (int ireg = 0; ireg < nreg; ++ireg) {
1298  int map = idet + ireg * ndet;
1299  int map_low = map + irf_low * ndet * nreg;
1300  int map_up = map + irf_up * ndet * nreg;
1301  for (int i = 0; i < npix; ++i) {
1302  irf(i, map) += wgt_low * m_irfs(i, map_low) +
1303  wgt_up * m_irfs(i, map_up);
1304  }
1305  }
1306  }
1307 
1308  } // endfor: looped over integration intervals
1309 
1310  // Debug: print weighting check
1311  #if defined(G_DEBUG_COMPUTE_IRF)
1312  std::cout << "- wgt_check = " << wgt_check << " (should be 1)" << std::endl;
1313  #endif
1314 
1315  } // endif: integrated response over energy band
1316 
1317  // Case B: compute response for line energy
1318  else {
1319 
1320  // Get IRFs bracketing the minimum energy and the IRF weights
1321  m_energies.set_value(emin);
1322  int irf_low = m_energies.inx_left();
1323  int irf_up = m_energies.inx_right();
1324  double wgt_low = m_energies.wgt_left();
1325  double wgt_up = m_energies.wgt_right();
1326 
1327  // Compute all IRF pixels
1328  for (int idet = 0; idet < ndet; ++idet) {
1329  for (int ireg = 0; ireg < nreg; ++ireg) {
1330  int map = idet + ireg * ndet;
1331  int map_low = map + irf_low * ndet * nreg;
1332  int map_up = map + irf_up * ndet * nreg;
1333  for (int i = 0; i < npix; ++i) {
1334  irf(i, map) = wgt_low * m_irfs(i, map_low) +
1335  wgt_up * m_irfs(i, map_up);
1336  }
1337  }
1338  }
1339 
1340  // Debug: print weighting check
1341  #if defined(G_DEBUG_COMPUTE_IRF)
1342  std::cout << "- wgt_low = " << wgt_low << std::endl;
1343  std::cout << "- wgt_up = " << wgt_up << std::endl;
1344  std::cout << "- wgt_low+wgt_up = " << wgt_low+wgt_up;
1345  std::cout << " (should be 1)" << std::endl;
1346  #endif
1347 
1348  } // endelse: compute response for line energy
1349 
1350  // Return IRF
1351  return irf;
1352 }
1353 
1354 
1355 /***********************************************************************//**
1356  * @brief Set IRF image limits
1357  *
1358  * @param[in] image IRF FITS image.
1359  *
1360  * Computes the IRF image limits.
1361  ***************************************************************************/
1363 {
1364  // Get image attributes
1365  int naxis1 = image->integer("NAXIS1");
1366  int naxis2 = image->integer("NAXIS2");
1367  double crval1 = image->real("CRVAL1");
1368  double crval2 = image->real("CRVAL2");
1369  double crpix1 = image->real("CRPIX1");
1370  double crpix2 = image->real("CRPIX2");
1371  double cdelt1 = image->real("CDELT1");
1372  double cdelt2 = image->real("CDELT2");
1373 
1374  // Derive image limits. Limits are taken at the pixel centres since
1375  // we want to use them for bilinear interpolation. This means that
1376  // we will throw away half a pixel at the edge of the IRFs.
1377  m_wcs_xmin = (crval1 - (crpix1-1.0) * cdelt1) * gammalib::deg2rad;
1378  m_wcs_ymin = (crval2 - (crpix2-1.0) * cdelt2) * gammalib::deg2rad;
1379  m_wcs_xbin = cdelt1 * gammalib::deg2rad;
1380  m_wcs_ybin = cdelt2 * gammalib::deg2rad;
1381  m_wcs_xmax = m_wcs_xmin + double(naxis1-1) * m_wcs_xbin;
1382  m_wcs_ymax = m_wcs_ymin + double(naxis2-1) * m_wcs_ybin;
1383  m_wcs_xpix_max = double(naxis1-1);
1384  m_wcs_ypix_max = double(naxis2-1);
1385 
1386  // Set maximum zenith angle
1389 
1390  // Signal that image limits exist
1391  m_has_wcs = true;
1392 
1393  // Return
1394  return;
1395 }
1396 
1397 
1398 /***********************************************************************//**
1399  * @brief Set vector of detector identifiers used by the observation
1400  *
1401  * @param[in] cube INTEGRAL/SPI event cube.
1402  *
1403  * Sets the vector of detector identifiers that is used by the observation.
1404  * The method scans the detector identifiers for all pointings and builds
1405  * a vector of unique detector identifiers in the order they were encountered
1406  * in the event cube.
1407  ***************************************************************************/
1409 {
1410  // Clear vector of detector identifiers
1411  m_detids.clear();
1412 
1413  // Extract relevant event cube dimensions
1414  int npt = cube->naxis(0);
1415  int ndet = cube->naxis(1);
1416 
1417  // Loop over all pointings and detectors
1418  for (int ipt = 0; ipt < npt; ++ipt) {
1419  for (int idet = 0; idet < ndet; ++idet) {
1420 
1421  // Get IRF detector identifier
1422  int detid = irf_detid(cube->dir(ipt, idet).detid());
1423 
1424  // Push IRF detector identifier on vector if it does not yet exist
1425  std::vector<int>::iterator it = find(m_detids.begin(), m_detids.end(), detid);
1426  if (it == m_detids.end()) {
1427  m_detids.push_back(detid);
1428  }
1429 
1430  } // endfor: looped over detectors
1431  } // endfor: looped over pointings
1432 
1433  // Return
1434  return;
1435 }
1436 
1437 
1438 /***********************************************************************//**
1439  * @brief Set computation cache
1440  *
1441  * @param[in] cube INTEGRAL/SPI event cube.
1442  *
1443  * Setup of two vectors for fast coordinate transformation into the
1444  * instrument system. The first vector m_spix stores the SPI pointing
1445  * direction (the X direction) while the second vector stores the position
1446  * angle in celestial coordinates of the SPI Y direction.
1447  ***************************************************************************/
1449 {
1450  // Extract relevant event cube dimensions
1451  int npt = cube->naxis(0);
1452 
1453  // Clear vectors of SPI X direction and position angles
1454  m_spix.clear();
1455  m_spix.reserve(npt);
1456  m_posang.clear();
1457  m_posang.reserve(npt);
1458 
1459  // Loop over all pointings
1460  for (int ipt = 0; ipt < npt; ++ipt) {
1461 
1462  // Compute position angle is celestial coordinates
1463  double posang = cube->spi_x(ipt).posang(cube->spi_z(ipt)) + gammalib::pihalf;
1464 
1465  // Store angles
1466  m_spix.push_back(cube->spi_x(ipt));
1467  m_posang.push_back(posang);
1468 
1469  } // endfor: looped over pointings
1470 
1471  // Return
1472  return;
1473 }
1474 
1475 
1476 /***********************************************************************//**
1477  * @brief Convert detector identifier into IRF detector identifier
1478  *
1479  * @param[in] detid SPI event detector identifier.
1480  * @return IRF detector identifier
1481  *
1482  * Converts a SPI event detector identifier into an IRF detector identifier.
1483  * TODO: Describe how this is done and why
1484  ***************************************************************************/
1485 int GSPIResponse::irf_detid(const int& detid) const
1486 {
1487  // Initialise irf detector identifier
1488  int irf_detid = detid;
1489 
1490  // Put detector identifier in the relevant range
1491  if (irf_detid >= 123) {
1492  irf_detid -= 123;
1493  }
1494  else if (irf_detid >= 104) {
1495  irf_detid -= 104;
1496  }
1497  else if (irf_detid >= 85) {
1498  irf_detid -= 85;
1499  }
1500 
1501  // Return irf detector identifier
1502  return irf_detid;
1503 }
1504 
1505 
1506 /***********************************************************************//**
1507  * @brief Compute weight of logarithmic energy bin
1508  *
1509  * @param[in] beta 1-gamma.
1510  * @param[in] emin Logarithmic minimum energy.
1511  * @param[in] emax Logarithmic maximum energy.
1512  ***************************************************************************/
1513 double GSPIResponse::irf_weight(const double& beta,
1514  const double& emin,
1515  const double& emax) const
1516 {
1517  // Initialise weight
1518  double weight;
1519 
1520  // Assign weight
1521  if (std::abs(beta) < DBL_MIN) {
1522  weight = (emax - emin);
1523  }
1524  else {
1525  weight = std::exp(beta * gammalib::ln10 * emax) -
1526  std::exp(beta * gammalib::ln10 * emin);
1527  }
1528 
1529  // Return weight
1530  return weight;
1531 }
#define G_LOAD_IRFS
void unit(const std::string &unit)
Set column unit.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
Sky map class.
Definition: GSkyMap.hpp:89
double m_wcs_xmin
Minimum X value (radians)
double m_wcs_ymax
Maximum Y value (radians)
FITS table double column class interface definition.
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
double azimuth(const int &ipt, const GSkyDir &dir) const
Return azimuth angle of sky direction for pointing in radians.
virtual void clear(void)
Clear instance.
GSkyMap m_irfs
IRFs stored as sky map.
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:379
const GSkyDir & spi_z(const int &ipt) const
Return SPI Z direction.
Point source spatial model class interface definition.
void set_cache(const GSPIEventCube *cube)
Set computation cache.
#define G_IRF
const GSkyDir & spi_x(const int &ipt) const
Return SPI X direction (pointing direction)
double m_max_zenith
Maximum zenith angle (radians)
const GSPIInstDir & dir(const int &ipt, const int &idet) const
Return instrument direction.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
std::vector< GSkyDir > m_spix
SPI pointing direction.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:301
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
const double & livetime(void) const
Return livetime of event bin.
double m_energy_keV
IRF line energy (optional)
Abstract interface for the event classes.
Definition: GEvent.hpp:71
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
virtual GEbounds ebounds(const GEnergy &obsEnergy) const
Return true energy boundaries for a specific observed energy.
virtual ~GSPIResponse(void)
Destructor.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void reserve(const int &num)
Reserves space for nodes in node array.
Definition: GNodeArray.hpp:220
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
INTEGRAL/SPI event bin container class definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print INTEGRAL/SPI response information.
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:761
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
void write_detids(GFits &fits) const
Write detector identifiers into FITS object.
void detid(const int &detid)
Set detector identifier.
FITS file class.
Definition: GFits.hpp:63
int spi_num_hdus(const GFits &fits, const std::string &extname)
Return number of HDU versions.
Definition: GSPITools.cpp:166
FITS file class interface definition.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
virtual const GSPIInstDir & dir(void) const
Return instrument direction.
double m_wcs_ybin
Y value bin size (radians)
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:787
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2664
Source class definition.
const std::vector< int > & shape(void) const
Returns shape of maps.
Definition: GSkyMap.hpp:403
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return value of INTEGRAL/SPI instrument response for a photon.
void write(GFits &file, const std::string &extname=gammalib::extname_ebounds, const std::string &unit="keV") const
Write energy boundaries into FITS object.
Definition: GEbounds.cpp:816
int irf_detid(const int &detid) const
Convert detector identifier into IRF detector identifier.
const double ln10
Definition: GMath.hpp:46
double m_dlogE
Logarithmic energy step for IRF band.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
void write_energies(GFits &fits) const
Write energies into FITS object.
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
void set_wcs(const GFitsImage *image)
Set IRF image limits.
INTEGRAL/SPI event bin container class.
virtual GSPIResponse * clone(void) const
Clone instance.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
double zenith(const int &ipt, const GSkyDir &dir) const
Return zenith angle of sky direction for pointing in radians.
std::string path(void) const
Return access path.
Definition: GFilename.hpp:217
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
Node array class interface definition.
void read_detids(const GFits &fits)
Read detector identifiers from FITS object.
void copy_members(const GSPIResponse &rsp)
Copy class members.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
#define G_EBOUNDS
Energy boundaries container class.
Definition: GEbounds.hpp:60
double m_wcs_xpix_max
Maximum X pixel index.
const GFilename & rspname(void) const
Get response group file name.
void set_detids(const GSPIEventCube *cube)
Set vector of detector identifiers used by the observation.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:391
double m_wcs_xmax
Maximum X value (radians)
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2697
Filename class.
Definition: GFilename.hpp:62
Energy container class definition.
Abstract interface for FITS table column.
void read_energies(const GFits &fits)
Read energies from FITS object.
const GFitsTable * spi_hdu(const GFits &fits, const std::string &extname, const int &extver=1)
Return FITS table.
Definition: GSPITools.cpp:57
double m_wcs_xbin
X value bin size (radians)
const int & ipt(void) const
Return event bin pointing index.
bool is_empty(void) const
Signal if there are no energy boundaries.
Definition: GEbounds.hpp:175
virtual double pixel(const int &ix) const =0
virtual GSPIResponse & operator=(const GSPIResponse &rsp)
Assignment operator.
const std::string extname_ebounds
Definition: GEbounds.hpp:44
std::vector< int > m_detids
Vector of detector IDs.
GSkyMap compute_irf(const double &emin, const double &emax) const
Compute as sky map.
double irf_weight(const double &beta, const double &emin, const double &emax) const
Compute weight of logarithmic energy bin.
void save(const GFilename &filename, const bool &clobber=false) const
Save SPI response into file.
void load_irfs(const int &region)
Load Instrument Response Functions.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GNodeArray m_energies
Node array of IRF energies.
double m_wcs_ymin
Minimum Y value (radians)
GChatter
Definition: GTypemaps.hpp:33
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
GSkyMap load_irf(const GFilename &rspname) const
Load IRF as sky map.
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
double irf_value(const GSkyDir &srcDir, const GSPIEventBin &bin, const int &ireg) const
Return value of INTEGRAL/SPI instrument response for sky direction and event bin. ...
Abstract observation base class.
INTEGRAL/SPI instrument response function class definition.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:268
double m_wcs_ypix_max
Maximum Y pixel index.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
void free_members(void)
Delete class members.
Definition: GResponse.cpp:835
const std::string extname_energies
Definition: GEnergies.hpp:44
std::vector< double > m_posang
Position angle of Y axis (CEL, radians)
INTEGRAL/SPI observation class.
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:139
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1151
const int & iebin(void) const
Return event bin energy index.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:368
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
virtual int integer(const int &row, const int &inx=0) const =0
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
Sky model class.
Definition: GModelSky.hpp:122
virtual double real(const int &row, const int &inx=0) const =0
INTEGRAL/SPI observation class definition.
double keV(void) const
Return energy in keV.
Definition: GEnergy.cpp:306
FITS binary table class.
GSPIResponse(void)
Void constructor.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
virtual GEvents * events(void)
Return events.
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
Energy boundaries class interface definition.
INTEGRAL/SPI event bin class.
Exception handler interface definition.
void read(const GFits &fits)
Read SPI response from FITS object.
GFilename m_rspname
File name of response group.
GEbounds m_ebounds
Energy bounaries of IRF.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
INTEGRAL/SPI instrument response function class.
FITS binary table class definition.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:363
double m_gamma
Power-law spectral index for IRF band.
bool m_has_wcs
Has WCS information.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
FITS table short integer column class interface definition.
virtual int naxis(const int &axis) const
Return number of bins in axis.
FITS table short integer column.
void init_members(void)
Initialise class members.
const double rad2deg
Definition: GMath.hpp:44
void load(const GFilename &filename)
Load SPI response from file.
Abstract instrument response base class.
Definition: GResponse.hpp:77
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
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:347
void write(GFits &fits) const
Write SPI response into FITS object.
#define G_LOAD_IRF
INTEGRAL/SPI event bin class definition.
virtual double nroi(const GModelSky &model, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return integral of event probability for a given sky model over ROI.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
void free_members(void)
Delete class members.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
Definition of INTEGRAL/SPI tools.
FITS table double column.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1295
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
#define G_NROI
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void set(const GSPIObservation &obs, const GEnergy &energy=GEnergy())
Set response for a specific observation.