Loading [MathJax]/extensions/MathZoom.js
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-2024 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 <vector>
33#include <typeinfo>
34#include <algorithm>
35#include "GException.hpp"
36#include "GTools.hpp"
37#include "GEbounds.hpp"
38#include "GEnergies.hpp"
39#include "GNodeArray.hpp"
40#include "GVector.hpp"
41#include "GMatrix.hpp"
42#include "GIntegrals.hpp"
43#include "GFits.hpp"
44#include "GFitsBinTable.hpp"
45#include "GFitsTableShortCol.hpp"
46#include "GFitsTableDoubleCol.hpp"
47#include "GSource.hpp"
48#include "GModelSky.hpp"
49#include "GModelSpatialPointSource.hpp"
50#include "GModelSpatialRadial.hpp"
51#include "GModelSpatialElliptical.hpp"
52#include "GModelSpatialDiffuse.hpp"
53#include "GModelSpatialDiffuseMap.hpp"
54#include "GSPITools.hpp"
55#include "GSPIResponse.hpp"
56#include "GSPIObservation.hpp"
57#include "GSPIEventCube.hpp"
58#include "GSPIEventBin.hpp"
59#include "spi_helpers_response_vector.hpp"
60
61/* __ Method name definitions ____________________________________________ */
62#define G_IRF "GSPIResponse::irf(GEvent&, GSource&, GObservation&)"
63#define G_NROI "GSPIResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
64 "GObservation&)"
65#define G_EBOUNDS "GSPIResponse::ebounds(GEnergy&)"
66#define G_SET "GSPIResponse::set(const GSPIObservation& obs)"
67#define G_LOAD_IRF "GSPIResponse::load_irf(GFilename&)"
68#define G_LOAD_IRFS "GSPIResponse::load_irfs(int&)"
69
70/* __ Macros _____________________________________________________________ */
71
72/* __ Coding definitions _________________________________________________ */
73#define G_IRF_RADIAL_PROJECTION //!< Use IRF projection for radial response
74#define G_IRF_ELLIPTICAL_PROJECTION //!< Use IRF projection for elliptical response
75#define G_IRF_DIFFUSE_PROJECTION //!< Use IRF projection for diffuse response
76
77/* __ Debug definitions __________________________________________________ */
78//#define G_DEBUG_COMPUTE_IRF //!< Debug compute_irf() method
79//#define G_DEBUG_IRF_SKYDIR //!< Debug irf_skydir method
80//#define G_DEBUG_IRF_TIMING //!< Time response computation
81
82/* __ Constants __________________________________________________________ */
83
84
85/*==========================================================================
86 = =
87 = Constructors/destructors =
88 = =
89 ==========================================================================*/
90
91/***********************************************************************//**
92 * @brief Void constructor
93 *
94 * Creates an empty INTEGRAL/SPI response.
95 ***************************************************************************/
96GSPIResponse::GSPIResponse(void) : GResponse()
97{
98 // Initialise members
99 init_members();
100
101 // Return
102 return;
103}
104
105
106/***********************************************************************//**
107 * @brief Response constructor
108 *
109 * Constructs an INTEGRAL/SPI response for an observation using a response
110 * group filename.
111 *
112 * This constructor simply stores the file name of a response group, the
113 * actual loading will be done using the set() method.
114 ***************************************************************************/
115GSPIResponse::GSPIResponse(const GFilename& rspname) : GResponse()
116{
117 // Initialise members
118 init_members();
119
120 // Store response group filename
121 m_rspname = rspname;
122
123 // Return
124 return;
125}
126
127
128/***********************************************************************//**
129 * @brief Copy constructor
130 *
131 * @param[in] rsp INTEGRAL/SPI response.
132 **************************************************************************/
133GSPIResponse::GSPIResponse(const GSPIResponse& rsp) : GResponse(rsp)
134{
135 // Initialise members
136 init_members();
137
138 // Copy members
139 copy_members(rsp);
140
141 // Return
142 return;
143}
144
145
146/***********************************************************************//**
147 * @brief Destructor
148 ***************************************************************************/
149GSPIResponse::~GSPIResponse(void)
150{
151 // Free members
152 free_members();
153
154 // Return
155 return;
156}
157
158
159/*==========================================================================
160 = =
161 = Operators =
162 = =
163 ==========================================================================*/
164
165/***********************************************************************//**
166 * @brief Assignment operator
167 *
168 * @param[in] rsp INTEGRAL/SPI response.
169 * @return INTEGRAL/SPI response.
170 *
171 * Assign INTEGRAL/SPI response to this object. The assignment performs
172 * a deep copy of all information, hence the original object from which the
173 * assignment was made can be destroyed after this operation without any loss
174 * of information.
175 ***************************************************************************/
176GSPIResponse& GSPIResponse::operator=(const GSPIResponse& rsp)
177{
178 // Execute only if object is not identical
179 if (this != &rsp) {
180
181 // Copy base class members
182 this->GResponse::operator=(rsp);
183
184 // Free members
185 free_members();
186
187 // Initialise members
188 init_members();
189
190 // Copy members
191 copy_members(rsp);
192
193 } // endif: object was not identical
194
195 // Return this object
196 return *this;
197}
198
199
200/*==========================================================================
201 = =
202 = Public methods =
203 = =
204 ==========================================================================*/
205
206/***********************************************************************//**
207 * @brief Clear instance
208 *
209 * Clears INTEGRAL/SPI response by resetting all class members to an initial
210 * state. Any information that was present before will be lost.
211 ***************************************************************************/
212void GSPIResponse::clear(void)
213{
214 // Free class members (base and derived classes, derived class first)
215 free_members();
216 this->GResponse::free_members();
217
218 // Initialise members
219 this->GResponse::init_members();
220 init_members();
221
222 // Return
223 return;
224}
225
226
227/***********************************************************************//**
228 * @brief Clone instance
229 *
230 * @return Pointer to deep copy of INTEGRAL/SPI response.
231 ***************************************************************************/
232GSPIResponse* GSPIResponse::clone(void) const
233{
234 return new GSPIResponse(*this);
235}
236
237
238/***********************************************************************//**
239 * @brief Return value of INTEGRAL/SPI instrument response for a photon
240 *
241 * @param[in] event Observed event.
242 * @param[in] photon Incident photon.
243 * @param[in] obs Observation.
244 * @return Instrument response \f$(cm^2 sr^{-1})\f$
245 *
246 * @exception GException::invalid_argument
247 * Observation is not a INTEGRAL/SPI observation.
248 *
249 * Returns the instrument response function for a given observed photon
250 * direction as function of the assumed true photon direction. The result
251 * is given by
252 *
253 * \f[
254 * R(p'|p) =
255 * \f]
256 *
257 * @todo Write down formula
258 * @todo Describe in detail how the response is computed.
259 ***************************************************************************/
260double GSPIResponse::irf(const GEvent& event,
261 const GPhoton& photon,
262 const GObservation& obs) const
263{
264 // Initialise IRF
265 double irf = 0.0;
266
267 // Extract INTEGRAL/SPI observation
268 const GSPIObservation* spi_obs = dynamic_cast<const GSPIObservation*>(&obs);
269 if (spi_obs == NULL) {
270 std::string cls = std::string(typeid(&obs).name());
271 std::string msg = "Observation of type \""+cls+"\" is not an "
272 "INTEGRAL/SPI observation. Please specify an "
273 "INTEGRAL/SPI observation as argument.";
274 throw GException::invalid_argument(G_IRF, msg);
275 }
276
277 // Extract INTEGRAL/SPI event cube
278 const GSPIEventCube* cube = dynamic_cast<const GSPIEventCube*>(spi_obs->events());
279 if (cube == NULL) {
280 std::string msg = "INTEGRAL/SPI observation does not contain a valid "
281 "event cube. Please specify an observation with an "
282 "event cube as argument.";
283 throw GException::invalid_argument(G_IRF, msg);
284 }
285
286 // Extract INTEGRAL/SPI event bin
287 const GSPIEventBin* bin = dynamic_cast<const GSPIEventBin*>(&event);
288 if (bin == NULL) {
289 std::string cls = std::string(typeid(&event).name());
290 std::string msg = "Event of type \""+cls+"\" is not an INTEGRAL/SPI "
291 "event. Please specify an INTEGRAL/SPI event as "
292 "argument.";
293 throw GException::invalid_argument(G_IRF, msg);
294 }
295
296 // Continue only if livetime for event is positive
297 if (bin->livetime() > 0.0) {
298
299 // Get energy bins of event and photon. Since only the photo peak
300 // response is supported so far, both indices need to be identical.
301 // Continue only if this is the case.
302 if (cube->ebounds().index(photon.energy()) == bin->iebin()) {
303
304 // Get IRF value for photo peak
305 irf = irf_value(photon.dir(), *bin, 0);
306
307 // Compile option: Check for NaN/Inf
308 #if defined(G_NAN_CHECK)
309 if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
310 std::cout << "*** ERROR: GSPIResponse::irf:";
311 std::cout << " NaN/Inf encountered";
312 std::cout << " (irf=" << irf;
313 std::cout << ")";
314 std::cout << std::endl;
315 }
316 #endif
317
318 } // endif: photon energy in same energy bin as event energy
319
320 } // endif: livetime of event was positive
321
322 // Return IRF value
323 return irf;
324}
325
326
327/***********************************************************************//**
328 * @brief Return value of INTEGRAL/SPI instrument response for sky direction
329 * and event bin
330 *
331 * @param[in] srcDir Sky direction.
332 * @param[in] bin INTEGRAL/SPI event bin.
333 * @param[in] ireg IRF region (0 = photo peak).
334 * @return Instrument response \f$(cm^2 sr^{-1})\f$
335 *
336 * Returns the instrument response function for a given sky direction and
337 * event bin. The value of the IRF is bilinearly interpolated from the
338 * pre-computed IRFs cube that is stored in the class.
339 ***************************************************************************/
340double GSPIResponse::irf_value(const GSkyDir& srcDir,
341 const GSPIEventBin& bin,
342 const int& ireg) const
343{
344 // Initialise IRF value
345 double irf = 0.0;
346
347 // Convert sky direction to zenith angle. Continue only if zenith
348 // angle is smaller than the maximum zenith angle.
349 double zenith = this->zenith(bin.ipt(), srcDir);
350 if (zenith < m_max_zenith) {
351
352 // Convert sky direction to azimuth angle
353 double azimuth = this->azimuth(bin.ipt(), srcDir);
354
355 // Compute pixel
356 double xpix = (zenith * std::cos(azimuth) - m_wcs_xmin) / m_wcs_xbin;
357 double ypix = (zenith * std::sin(azimuth) - m_wcs_ymin) / m_wcs_ybin;
358
359 // Continue only if pixel is within IRF
360 if (xpix > 0.0 && xpix < m_wcs_xpix_max &&
361 ypix > 0.0 && ypix < m_wcs_ypix_max) {
362
363 // Get number of pixels in X direction
364 int nx = m_irfs.nx();
365 int ndet = m_irfs.shape()[0];
366 int nreg = m_irfs.shape()[1];
367
368 // Get IRF detector for event
369 int idet = irf_detid(bin.dir().detid());
370 int ieng = bin.iebin();
371
372 // Get map
373 int map = idet + (ireg + ieng * nreg) * ndet;
374
375 // Get 4 nearest neighbours
376 int ix_left = int(xpix);
377 int ix_right = ix_left + 1;
378 int iy_top = int(ypix);
379 int iy_bottom = iy_top + 1;
380
381 // Get weighting factors
382 double wgt_right = xpix - double(ix_left);
383 double wgt_left = 1.0 - wgt_right;
384 double wgt_bottom = ypix - double(iy_top);
385 double wgt_top = 1.0 - wgt_bottom;
386
387 // Get indices of 4 nearest neighbours
388 int inx_1 = ix_left + iy_top * nx;
389 int inx_2 = ix_right + iy_top * nx;
390 int inx_3 = ix_left + iy_bottom * nx;
391 int inx_4 = ix_right + iy_bottom * nx;
392
393 // Get weights of 4 nearest neighbours
394 double wgt_1 = wgt_left * wgt_top;
395 double wgt_2 = wgt_right * wgt_top;
396 double wgt_3 = wgt_left * wgt_bottom;
397 double wgt_4 = wgt_right * wgt_bottom;
398
399 // Compute IRF
400 irf = m_irfs(inx_1, map) * wgt_1;
401 irf += m_irfs(inx_2, map) * wgt_2;
402 irf += m_irfs(inx_3, map) * wgt_3;
403 irf += m_irfs(inx_4, map) * wgt_4;
404
405 // Make sure that IRF does not get negative
406 if (irf < 0.0) {
407 irf = 0.0;
408 }
409
410 } // endif: zenith angle was valid
411
412 } // endif: pixel was within IRF
413
414 // Return IRF value
415 return irf;
416}
417
418
419/***********************************************************************//**
420 * @brief Return integral of event probability for a given sky model over ROI
421 *
422 * @param[in] model Sky model.
423 * @param[in] obsEng Observed photon energy.
424 * @param[in] obsTime Observed photon arrival time.
425 * @param[in] obs Observation.
426 * @return 0.0
427 *
428 * @exception GException::feature_not_implemented
429 * Method is not implemented.
430 ***************************************************************************/
431double GSPIResponse::nroi(const GModelSky& model,
432 const GEnergy& obsEng,
433 const GTime& obsTime,
434 const GObservation& obs) const
435{
436 // Method is not implemented
437 std::string msg = "Spatial integration of sky model over the data space "
438 "is not implemented.";
439 throw GException::feature_not_implemented(G_NROI, msg);
440
441 // Return Npred
442 return (0.0);
443}
444
445
446/***********************************************************************//**
447 * @brief Return true energy boundaries for a specific observed energy
448 *
449 * @param[in] obsEnergy Observed Energy.
450 * @return True energy boundaries for given observed energy.
451 *
452 * @exception GException::feature_not_implemented
453 * Method is not implemented.
454 *
455 * @todo Implement this method if you need energy dispersion.
456 ***************************************************************************/
457GEbounds GSPIResponse::ebounds(const GEnergy& obsEnergy) const
458{
459 // Initialise an empty boundary object
460 GEbounds ebounds;
461
462 // Throw an exception
463 std::string msg = "Energy dispersion not implemented.";
464 throw GException::feature_not_implemented(G_EBOUNDS, msg);
465
466 // Return energy boundaries
467 return ebounds;
468}
469
470
471/***********************************************************************//**
472 * @brief Set response for a specific observation
473 *
474 * @param[in] obs INTEGRAL/SPI observation.
475 * @param[in] energy Line energy
476 *
477 * Sets the response for a specific INTEGRAL/SPI observation. This is the
478 * high-level method for loading and pre-computing SPI IRFs for a specific
479 * energy binning. The method loads the IRF from a response group file and
480 * pre-computes the IRFs for all energy bins of the SPI observation. If
481 * @p energy is positive, IRFs will be pre-computed taking a line energy
482 * specified by @p energy, otherwise the IRFs will be pre-computed taking
483 * the energy interval of the energy bins (continuum IRFs).
484 ***************************************************************************/
485void GSPIResponse::set(const GSPIObservation& obs, const GEnergy& energy)
486{
487 // Reset response information
488 m_detids.clear();
489 m_energies.clear();
490 m_ebounds.clear();
491 m_irfs.clear();
492 m_has_wcs = false;
493
494 // Continue only if the observation contains an event cube
495 const GSPIEventCube* cube = dynamic_cast<const GSPIEventCube*>(obs.events());
496 if (cube != NULL) {
497
498 // Set requested detector identifiers
499 set_detids(cube);
500
501 // Set coordination transformation cache
502 set_cache(cube);
503
504 // Load IRFs for photo peak region (region 0)
505 load_irfs(0);
506
507 // Get dimension of pre-computed IRF
508 int npix = m_irfs.npix();
509 int ndet = m_irfs.shape()[0];
510 int nreg = m_irfs.shape()[1];
511 int neng = cube->naxis(2);
512
513 // Initialise pre-computed IRF
514 GSkyMap irfs = m_irfs;
515 irfs.nmaps(ndet * nreg * neng);
516 irfs.shape(ndet, nreg, neng);
517 irfs = 0.0;
518
519 // Pre-compute IRFs for all energy bins
520 for (int ieng = 0; ieng < neng; ++ieng) {
521
522 // Get energy boundaries
523 double emin = cube->ebounds().emin(ieng).keV();
524 double emax = cube->ebounds().emax(ieng).keV();
525
526 // If line energy is specified then compute a line response for
527 // the energy bin that overlaps with the line energy
528 m_energy_keV = energy.keV();
529 if (m_energy_keV > 0.0) {
530 if (m_energy_keV < emin || m_energy_keV > emax) {
531 continue;
532 }
533 emin = m_energy_keV;
534 emax = m_energy_keV;
535 }
536
537 // Pre-compute IRF for this energy bin
538 GSkyMap irf = compute_irf(emin, emax);
539
540 // Copy over IRF
541 for (int idet = 0; idet < ndet; ++idet) {
542 for (int ireg = 0; ireg < nreg; ++ireg) {
543 int map_irf = idet + ireg * ndet;
544 int map_irfs = idet + (ireg + ieng * nreg) * ndet;
545 for (int i = 0; i < npix; ++i) {
546 irfs(i, map_irfs) = irf(i, map_irf);
547 }
548 }
549 }
550
551 // Append energy boundary
552 m_ebounds.append(GEnergy(emin, "keV"), GEnergy(emax, "keV"));
553
554 } // endfor: looped over all energy bins
555
556 // Replace IRFs
557 m_irfs = irfs;
558
559 // Clear energies
560 m_energies.clear();
561
562 } // endif: observation contained an event cube
563
564 // Return
565 return;
566}
567
568
569/***********************************************************************//**
570 * @brief Read SPI response from FITS object
571 *
572 * @param[in] fits FITS object.
573 *
574 * Reads the SPI response from FITS object.
575 ***************************************************************************/
576void GSPIResponse::read(const GFits& fits)
577{
578 // Get IRF image
579 const GFitsImage* image = fits.image(0);
580
581 // Read IRFs
582 m_irfs.read(*image);
583
584 // Set WCS image limits
585 set_wcs(image);
586
587 // Read detector identifiers
588 read_detids(fits);
589
590 // Read energies
591 read_energies(fits);
592
593 // Return
594 return;
595}
596
597
598/***********************************************************************//**
599 * @brief Write SPI response into FITS object
600 *
601 * @param[in,out] fits FITS object.
602 *
603 * Writes the SPI response into FITS object.
604 ***************************************************************************/
605void GSPIResponse::write(GFits& fits) const
606{
607 // Write IRFs as sky map
608 m_irfs.write(fits);
609
610 // Write detector identifiers
611 write_detids(fits);
612
613 // Write energies
614 write_energies(fits);
615
616 // Return
617 return;
618}
619
620
621/***********************************************************************//**
622 * @brief Load SPI response from file
623 *
624 * @param[in] filename Response file name.
625 *
626 * Loads SPI response from response file.
627 ***************************************************************************/
628void GSPIResponse::load(const GFilename& filename)
629{
630 // Open FITS file
631 GFits fits(filename);
632
633 // Read response
634 read(fits);
635
636 // Close FITS file
637 fits.close();
638
639 // Return
640 return;
641}
642
643
644/***********************************************************************//**
645 * @brief Save SPI response into file
646 *
647 * @param[in] filename Response file name.
648 * @param[in] clobber Overwrite existing FITS file?
649 *
650 * Saves SPI response into response file.
651 ***************************************************************************/
652void GSPIResponse::save(const GFilename& filename, const bool& clobber) const
653{
654 // Creat FITS file
655 GFits fits;
656
657 // Write response
658 write(fits);
659
660 // Save FITS file
661 fits.saveto(filename, clobber);
662
663 // Close FITS file
664 fits.close();
665
666 // Return
667 return;
668}
669
670
671/***********************************************************************//**
672 * @brief Print INTEGRAL/SPI response information
673 *
674 * @param[in] chatter Chattiness.
675 * @return String containing INTEGRAL/SPI response information.
676 ***************************************************************************/
677std::string GSPIResponse::print(const GChatter& chatter) const
678{
679 // Initialise result string
680 std::string result;
681
682 // Continue only if chatter is not silent
683 if (chatter != SILENT) {
684
685 // Extract information
686 int nx = m_irfs.nx();
687 int ny = m_irfs.ny();
688 int ndet = (m_irfs.shape().size() > 0) ? m_irfs.shape()[0] : 0;
689 int nreg = (m_irfs.shape().size() > 1) ? m_irfs.shape()[1] : 0;
690 int neng = (m_irfs.shape().size() > 2) ? m_irfs.shape()[2] : 0;
691
692 // Append header
693 result.append("=== GSPIResponse ===");
694
695 // Append information
696 result.append("\n"+gammalib::parformat("Response group filename"));
697 result.append(m_rspname.url());
698 result.append("\n"+gammalib::parformat("Response map pixels"));
699 result.append(gammalib::str(nx)+" * "+gammalib::str(ny));
700 result.append("\n"+gammalib::parformat("X axis range"));
701 result.append("["+gammalib::str(m_wcs_xmin * gammalib::rad2deg));
702 result.append(","+gammalib::str(m_wcs_xmax * gammalib::rad2deg));
703 result.append("] deg");
704 result.append("\n"+gammalib::parformat("Y axis range"));
705 result.append("["+gammalib::str(m_wcs_ymin * gammalib::rad2deg));
706 result.append(","+gammalib::str(m_wcs_ymax * gammalib::rad2deg));
707 result.append("] deg");
708 result.append("\n"+gammalib::parformat("Maximum zenith angle"));
709 result.append(gammalib::str(m_max_zenith * gammalib::rad2deg)+" deg");
710 result.append("\n"+gammalib::parformat("Number of detectors"));
711 result.append(gammalib::str(ndet));
712 result.append("\n"+gammalib::parformat("Number of regions"));
713 result.append(gammalib::str(nreg));
714 result.append("\n"+gammalib::parformat("Number of energies"));
715 result.append(gammalib::str(neng));
716 if (m_energy_keV > 0.0) {
717 result.append("\n"+gammalib::parformat("Line IRF energy"));
718 result.append(gammalib::str(m_energy_keV)+" keV");
719 }
720 else {
721 result.append("\n"+gammalib::parformat("Continuum IRF gamma"));
722 result.append(gammalib::str(m_gamma));
723 result.append("\n"+gammalib::parformat("Continnum IRF log(E) step"));
724 result.append(gammalib::str(m_dlogE));
725 }
726
727 } // endif: chatter was not silent
728
729 // Return result
730 return result;
731}
732
733
734/*==========================================================================
735 = =
736 = Private methods =
737 = =
738 ==========================================================================*/
739
740/***********************************************************************//**
741 * @brief Initialise class members
742 ***************************************************************************/
743void GSPIResponse::init_members(void)
744{
745 // Initialise members
746 m_rspname.clear();
747 m_detids.clear();
748 m_energies.clear();
749 m_ebounds.clear();
750 m_irfs.clear();
751 m_energy_keV = 0.0;
752 m_dlogE = 0.03;
753 m_gamma = 2.0;
754
755 // Initialise cache
756 m_spix.clear();
757 m_posang.clear();
758 m_has_wcs = false;
759 m_wcs_xmin = 0.0;
760 m_wcs_ymin = 0.0;
761 m_wcs_xmax = 0.0;
762 m_wcs_ymax = 0.0;
763 m_wcs_xbin = 0.0;
764 m_wcs_ybin = 0.0;
765 m_wcs_xpix_max = 0.0;
766 m_wcs_ypix_max = 0.0;
767 m_max_zenith = 180.0 * gammalib::deg2rad;
768
769 // Return
770 return;
771}
772
773
774/***********************************************************************//**
775 * @brief Copy class members
776 *
777 * @param[in] rsp INTEGRAL/SPI response function.
778 ***************************************************************************/
779void GSPIResponse::copy_members(const GSPIResponse& rsp)
780{
781 // Copy members
782 m_rspname = rsp.m_rspname;
783 m_detids = rsp.m_detids;
784 m_energies = rsp.m_energies;
785 m_ebounds = rsp.m_ebounds;
786 m_irfs = rsp.m_irfs;
787 m_energy_keV = rsp.m_energy_keV;
788 m_dlogE = rsp.m_dlogE;
789 m_gamma = rsp.m_gamma;
790
791 // Copy cache
792 m_spix = rsp.m_spix;
793 m_posang = rsp.m_posang;
794 m_has_wcs = rsp.m_has_wcs;
795 m_wcs_xmin = rsp.m_wcs_xmin;
796 m_wcs_ymin = rsp.m_wcs_ymin;
797 m_wcs_xmax = rsp.m_wcs_xmax;
798 m_wcs_ymax = rsp.m_wcs_ymax;
799 m_wcs_xbin = rsp.m_wcs_xbin;
800 m_wcs_ybin = rsp.m_wcs_ybin;
801 m_wcs_xpix_max = rsp.m_wcs_xpix_max;
802 m_wcs_ypix_max = rsp.m_wcs_ypix_max;
803 m_max_zenith = rsp.m_max_zenith;
804
805 // Return
806 return;
807}
808
809
810/***********************************************************************//**
811 * @brief Delete class members
812 ***************************************************************************/
813void GSPIResponse::free_members(void)
814{
815 // Return
816 return;
817}
818
819
820/***********************************************************************//**
821 * @brief Read detector identifiers from FITS object
822 *
823 * @param[in] fits FITS object.
824 *
825 * Read the detector identifiers from a FITS file into the m_detids member.
826 ***************************************************************************/
827void GSPIResponse::read_detids(const GFits& fits)
828{
829 // Clear detids vector
830 m_detids.clear();
831
832 // Get DETID table extension
833 const GFitsTable& table = *fits.table("DETIDS");
834 const GFitsTableCol* col_detid = table["DET_ID"];
835
836 // Get number of detector identifiers
837 int num = table.nrows();
838
839 // Reserve space for detector identifiers
840 m_detids.reserve(num);
841
842 // Extract values
843 for (int i = 0; i < num; ++i) {
844 m_detids.push_back(col_detid->integer(i));
845 }
846
847 // Return
848 return;
849}
850
851
852/***********************************************************************//**
853 * @brief Read energies from FITS object
854 *
855 * @param[in] fits FITS object.
856 *
857 * Read the energies from a FITS file. If the FITS file contains energy
858 * boundaries then read them, otherwise read the energy vector.
859 ***************************************************************************/
860void GSPIResponse::read_energies(const GFits& fits)
861{
862 // Clear energies and energy boundary vectors
863 m_energies.clear();
864 m_ebounds.clear();
865
866 // If FITS file contains energy boundaries then the IRF was precomputed
867 // and hence we read the energy boundaries from the FITS file
868 if (fits.contains(gammalib::extname_ebounds)) {
869
870 // Get reference to FITS table
871 const GFitsTable& table = *fits.table(gammalib::extname_ebounds);
872
873 // Read energy boundaries
874 m_ebounds.read(table);
875
876 // Set response attributes
877 m_energy_keV = (table.has_card("ENERGY")) ? table.real("ENERGY") : 0.0;
878 m_dlogE = (table.has_card("DLOGE")) ? table.real("DLOGE") : 0.03;
879 m_gamma = (table.has_card("GAMMA")) ? table.real("GAMMA") : 2.0;
880
881 } // endif: energy boundaries detected
882
883 // ... otherwise we have a response group and we read the energies of
884 // the response group
885 else {
886
887 // Get energies table extension
888 const GFitsTable& table = *fits.table(gammalib::extname_energies);
889 const GFitsTableCol* col_energy = table["ENERGY"];
890
891 // Get number of energies
892 int num = table.nrows();
893
894 // Reserve space for energies
895 m_energies.reserve(num);
896
897 // Get energy unit
898 std::string unit = "keV";
899 if (table.has_card("TUNIT1")) {
900 unit = table.string("TUNIT1");
901 }
902
903 // Extract values
904 for (int i = 0; i < num; ++i) {
905
906 // Get energy
907 GEnergy energy = GEnergy(col_energy->real(i), unit);
908
909 // Store as keV
910 m_energies.append(energy.keV());
911
912 } // endfor: looped over all values
913
914 } // endelse: energies extension was required
915
916 // Return
917 return;
918}
919
920
921/***********************************************************************//**
922 * @brief Write detector identifiers into FITS object
923 *
924 * @param[in,out] fits FITS object.
925 *
926 * Writes the detector identifiers stored into the m_detids member into a
927 * FITS object.
928 ***************************************************************************/
929void GSPIResponse::write_detids(GFits& fits) const
930{
931 // Get number of detector identifiers
932 int num = m_detids.size();
933
934 // Create columns
935 GFitsTableShortCol col_detid("DET_ID", num);
936
937 // Fill energy column in units of keV
938 for (int i = 0; i < num; ++i) {
939 col_detid(i) = m_detids[i];
940 }
941
942 // Create energies table
943 GFitsBinTable table(num);
944 table.append(col_detid);
945 table.extname("DETIDS");
946
947 // Append detector identifiers to FITS file
948 fits.append(table);
949
950 // Return
951 return;
952}
953
954
955/***********************************************************************//**
956 * @brief Write energies into FITS object
957 *
958 * @param[in,out] fits FITS object.
959 *
960 * Writes the energy nodes stored into the m_energies member into a FITS
961 * object. The energies are written in units of keV.
962 ***************************************************************************/
963void GSPIResponse::write_energies(GFits& fits) const
964{
965 // If there are energy boundaries then the IRF was precomputed and
966 // hence we store the energy boundaries in the FITS file
967 if (!m_ebounds.is_empty()) {
968
969 // Write energy boundaries
970 m_ebounds.write(fits);
971
972 // Recover the FITS table to write some keywords
973 GFitsTable& table = *fits.table(gammalib::extname_ebounds);
974
975 // Write keywords
976 table.card("ENERGY", m_energy_keV, "[keV] IRF line energy (0 for continuum IRF)");
977 table.card("DLOGE", m_dlogE, "Logarithmic step size for continuum IRF");
978 table.card("GAMMA", m_gamma, "Power-law index for continuum IRF");
979
980 } // endif: IRF was precomputed
981
982 // ... otherwise we write the energies
983 else {
984
985 // Get number of energies
986 int num = m_energies.size();
987
988 // Create columns
989 GFitsTableDoubleCol col_energy("ENERGY", num);
990
991 // Fill energy column in units of keV
992 for (int i = 0; i < num; ++i) {
993 col_energy(i) = m_energies[i];
994 }
995 col_energy.unit("keV");
996
997 // Create energies table
998 GFitsBinTable table(num);
999 table.append(col_energy);
1000 table.extname(gammalib::extname_energies);
1001
1002 // Append energies table to FITS file
1003 fits.append(table);
1004
1005 } // endelse: energies written
1006
1007 // Return
1008 return;
1009}
1010
1011
1012/***********************************************************************//**
1013 * @brief Load IRF as sky map
1014 *
1015 * @param[in] irfname IRF file name.
1016 *
1017 * Loads an IRF FITS file as a sky map. The sky map is returned in ARC
1018 * projection that is a zenith equidistant projection, which is the
1019 * projection that is natively used to store the IRFs.
1020 ***************************************************************************/
1021GSkyMap GSPIResponse::load_irf(const GFilename& irfname) const
1022{
1023 // Open IRF FITS file
1024 GFits fits(irfname);
1025
1026 // Access IRF FITS image
1027 GFitsImage* image = fits.image("SPI.-IRF.-RSP");
1028
1029 // Get image attributes
1030 int naxis1 = image->integer("NAXIS1");
1031 int naxis2 = image->integer("NAXIS2");
1032 int naxis3 = image->integer("NAXIS3");
1033 int naxis4 = image->integer("NAXIS4");
1034 double crval2 = image->real("CRVAL2");
1035 double crval3 = image->real("CRVAL3");
1036 double crpix2 = image->real("CRPIX2");
1037 double crpix3 = image->real("CRPIX3");
1038 double cdelt2 = image->real("CDELT2");
1039 double cdelt3 = image->real("CDELT3");
1040
1041 // Derive image limits. Limits are taken at the pixel centres since
1042 // we want to use them for bilinear interpolation. This means that
1043 // we will throw away half a pixel at the edge of the IRFs.
1044 double wcs_xmin = (crval2 - (crpix2-1.0) * cdelt2) * gammalib::deg2rad;
1045 double wcs_ymin = (crval3 - (crpix3-1.0) * cdelt3) * gammalib::deg2rad;
1046 double wcs_xbin = cdelt2 * gammalib::deg2rad;
1047 double wcs_ybin = cdelt3 * gammalib::deg2rad;
1048 double wcs_xmax = wcs_xmin + double(naxis2-1) * wcs_xbin;
1049 double wcs_ymax = wcs_ymin + double(naxis3-1) * wcs_ybin;
1050 double wcs_xpix_max = double(naxis2-1);
1051 double wcs_ypix_max = double(naxis3-1);
1052
1053 // If no image limits exists so far then store them for fast IRF access
1054 // that does not depend on the actual IRF projection (we just use here
1055 // the sky map as a convient container)
1056 if (!m_has_wcs) {
1057 m_has_wcs = true;
1058 m_wcs_xmin = wcs_xmin;
1059 m_wcs_ymin = wcs_ymin;
1060 m_wcs_xbin = wcs_xbin;
1061 m_wcs_ybin = wcs_ybin;
1062 m_wcs_xmax = wcs_xmax;
1063 m_wcs_ymax = wcs_ymax;
1064 m_wcs_xpix_max = wcs_xpix_max;
1065 m_wcs_ypix_max = wcs_ypix_max;
1066 }
1067
1068 // ... otherwise check if the limits are consistent
1069 else {
1070 if ((std::abs(wcs_xmin - m_wcs_xmin) > 1.0e-6) ||
1071 (std::abs(wcs_ymin - m_wcs_ymin) > 1.0e-6) ||
1072 (std::abs(wcs_xbin - m_wcs_xbin) > 1.0e-6) ||
1073 (std::abs(wcs_ybin - m_wcs_ybin) > 1.0e-6) ||
1074 (std::abs(wcs_xmax - m_wcs_xmax) > 1.0e-6) ||
1075 (std::abs(wcs_ymax - m_wcs_ymax) > 1.0e-6) ||
1076 (std::abs(wcs_xpix_max - m_wcs_xpix_max) > 1.0e-6) ||
1077 (std::abs(wcs_ypix_max - m_wcs_ypix_max) > 1.0e-6)) {
1078 std::string msg = "Inconsistent IRFs encountered in file \""+
1079 irfname.url()+"\". Please specify a response "
1080 "group where all IRFs have the same definition.";
1081 throw GException::invalid_value(G_LOAD_IRF, msg);
1082 }
1083 }
1084
1085 // Set maximum zenith angle
1086 m_max_zenith = (std::abs(m_wcs_xmax) > std::abs(m_wcs_ymax)) ?
1087 std::abs(m_wcs_xmax) : std::abs(m_wcs_ymax);
1088
1089 // Compute sky map attributes
1090 int nmap = naxis1 * naxis4;
1091
1092 // Initialise IRF in celestial coordinates. We use an ARC projection as
1093 // this is the native projection of the IRF files.
1094 GSkyMap irf("ARC", "CEL", crval2, crval3, cdelt2, cdelt3, naxis2, naxis3, nmap);
1095
1096 // Fill sky map
1097 for (int ix = 0; ix < naxis2; ++ix) {
1098 for (int iy = 0; iy < naxis3; ++iy) {
1099 for (int idet = 0; idet < naxis1; ++idet) {
1100 for (int ireg = 0; ireg < naxis4; ++ireg) {
1101 int index = ix + iy * naxis2;
1102 int map = idet + ireg * naxis1;
1103 irf(index,map) = image->pixel(idet, ix, iy, ireg);
1104 }
1105 }
1106 }
1107 }
1108
1109 // Set sky map shape
1110 irf.shape(naxis1, naxis4);
1111
1112 // Close FITS file
1113 fits.close();
1114
1115 // Return IRF
1116 return irf;
1117}
1118
1119
1120/***********************************************************************//**
1121 * @brief Load Instrument Response Functions
1122 *
1123 * @param[in] region IRF region (-1 = load all regions)
1124 *
1125 * The method requires that the required detector identifiers were previously
1126 * setup using the set_detids() method.
1127 ***************************************************************************/
1128void GSPIResponse::load_irfs(const int& region)
1129{
1130 // Initialise results
1131 m_energies.clear();
1132 m_irfs.clear();
1133
1134 // Determine number of requested detectors. Continue only if there are
1135 // required detectors.
1136 int num_det = m_detids.size();
1137 if (num_det > 0) {
1138
1139 // Open response group
1140 GFits rsp(m_rspname);
1141
1142 // Get response grouping table
1143 const GFitsTable* grp = gammalib::spi_hdu(rsp, "SPI.-IRF.-RSP-IDX");
1144 if (grp == NULL) {
1145 std::string msg = "No response grouping table found in file \""+
1146 m_rspname.url()+"\". Please specify a valid "
1147 "response grouping file.";
1148 throw GException::invalid_value(G_LOAD_IRFS, msg);
1149 }
1150
1151 // Determine number of response members
1152 int num_irfs = gammalib::spi_num_hdus(rsp, "SPI.-IRF.-RSP");
1153
1154 // Setup node array for response energies
1155 for (int i_irf = 0; i_irf < num_irfs; ++i_irf) {
1156 m_energies.append((*grp)["ENERGY"]->real(i_irf)); // in keV
1157 }
1158
1159 // Loop over all IRFs
1160 for (int i_irf = 0; i_irf < num_irfs; ++i_irf) {
1161
1162 // Get IRF filename
1163 std::string irfname = m_rspname.path() +
1164 (*grp)["MEMBER_LOCATION"]->string(i_irf);
1165
1166 // Load IRF
1167 GSkyMap irf = load_irf(irfname);
1168
1169 // Determine number of requested regions
1170 int num_regions = (region == -1) ? irf.shape()[1] : 1;
1171
1172 // If this is the first IRF then allocate IRFs
1173 if (i_irf == 0) {
1174
1175 // Set IRFs from IRF
1176 m_irfs = irf;
1177
1178 // Allocate sufficient maps and reshape maps
1179 m_irfs.nmaps(num_det * num_regions * num_irfs);
1180 m_irfs.shape(num_det, num_regions, num_irfs);
1181
1182 // Initialise maps
1183 m_irfs = 0.0;
1184
1185 } // endif: this was the first IRF
1186
1187 // Extract relevant part or IRF
1188 int nx = irf.nx();
1189 int ny = irf.ny();
1190 int ndet = irf.shape()[0];
1191
1192 // Loop over requested regions
1193 for (int i_region = 0; i_region < num_regions; ++i_region) {
1194
1195 // Loop over requested detectors
1196 for (int i_det = 0; i_det < num_det; ++i_det) {
1197
1198 // Set map indices in source IRF and target IRFs
1199 int map_irf = m_detids[i_det] + i_region * ndet;
1200 int map_irfs = i_det + (i_region + i_irf * num_regions) * num_det;
1201
1202 // Copy IRF pixels
1203 for (int i = 0, iy = 0; iy < ny; ++iy) {
1204 for (int ix = 0; ix < nx; ++ix, ++i) {
1205 m_irfs(i, map_irfs) = irf(i, map_irf);
1206 }
1207 }
1208
1209 } // endfor: looped over requested detectors
1210
1211 } // endfor: looped over requested regions
1212
1213 } // endfor: looped over all IRFs
1214
1215 // Close FITS file
1216 rsp.close();
1217
1218 } // endif: there were requested detectors
1219
1220 // Return
1221 return;
1222}
1223
1224
1225/***********************************************************************//**
1226 * @brief Compute as sky map
1227 *
1228 * @param[in] emin Minimum energy (keV).
1229 * @param[in] emax Maximum energy (keV).
1230 *
1231 * Compute IRF for an energy band specified by an energy interval. If the
1232 * width of the energy interval is zero a line IRF will be computed.
1233 ***************************************************************************/
1234GSkyMap GSPIResponse::compute_irf(const double& emin, const double& emax) const
1235{
1236 // Get IRF dimension
1237 int npix = m_irfs.npix();
1238 int ndet = m_irfs.shape()[0];
1239 int nreg = m_irfs.shape()[1];
1240
1241 // Debug: print dimensions
1242 #if defined(G_DEBUG_COMPUTE_IRF)
1243 std::cout << "GSPIResponse::compute_irf:" << std::endl;
1244 std::cout << "- emin = " << emin << " keV" << std::endl;
1245 std::cout << "- emax = " << emax << " keV" << std::endl;
1246 std::cout << "- npix = " << npix << std::endl;
1247 std::cout << "- ndet = " << ndet << std::endl;
1248 std::cout << "- nreg = " << nreg << std::endl;
1249 #endif
1250
1251 // Initialise IRF
1252 GSkyMap irf = m_irfs;
1253 irf.nmaps(ndet * nreg);
1254 irf.shape(ndet, nreg);
1255 irf = 0.0;
1256
1257 // Case A: Integrate response over energy band
1258 if (emax > emin) {
1259
1260 // Get logarithmic energy limits
1261 double log_emin = std::log10(emin);
1262 double log_emax = std::log10(emax);
1263 double log_ewidth = log_emax - log_emin;
1264
1265 // Set energy weighting factors
1266 double beta = 1.0 - m_gamma;
1267 double wgt0 = 1.0 / irf_weight(beta, log_emin, log_emax);
1268
1269 // Determine number of integration intervals and logarithmic step
1270 // size
1271 int num_int = (m_dlogE > DBL_MIN) ? int(log_ewidth / m_dlogE + 0.5) : 1;
1272 if (num_int < 1) {
1273 num_int = 1;
1274 }
1275 double log_estep = log_ewidth / num_int;
1276
1277 // Debug: print weighting and internal binning
1278 #if defined(G_DEBUG_COMPUTE_IRF)
1279 std::cout << "- log_emin = " << log_emin << std::endl;
1280 std::cout << "- log_emax = " << log_emax << std::endl;
1281 std::cout << "- log_ewidth = " << log_ewidth << std::endl;
1282 std::cout << "- beta = " << beta << std::endl;
1283 std::cout << "- wgt0 = " << wgt0 << std::endl;
1284 std::cout << "- num_int = " << num_int << std::endl;
1285 std::cout << "- log_estep = " << log_estep << std::endl;
1286 double wgt_check = 0.0;
1287 #endif
1288
1289 // Loop over integration intervals
1290 for (int i_int = 0; i_int < num_int; ++i_int) {
1291
1292 // Determine integration interval width and midpoint
1293 double log_emin_bin = log_emin + i_int * log_estep;
1294 double log_emax_bin = log_emin_bin + log_estep;
1295 double energy = std::exp(0.5 * gammalib::ln10 *
1296 (log_emin_bin + log_emax_bin));
1297
1298 // Get response weighting factor
1299 double wgt = irf_weight(beta, log_emin_bin, log_emax_bin) * wgt0;
1300
1301 // Get IRFs bracketing the mean energy and the IRF weights
1302 m_energies.set_value(energy);
1303 int irf_low = m_energies.inx_left();
1304 int irf_up = m_energies.inx_right();
1305 double wgt_low = m_energies.wgt_left() * wgt;
1306 double wgt_up = m_energies.wgt_right() * wgt;
1307
1308 // Debug: add weight for final check
1309 #if defined(G_DEBUG_COMPUTE_IRF)
1310 wgt_check += wgt_low + wgt_up;
1311 #endif
1312
1313 // Add contribution to all IRF pixels
1314 for (int idet = 0; idet < ndet; ++idet) {
1315 for (int ireg = 0; ireg < nreg; ++ireg) {
1316 int map = idet + ireg * ndet;
1317 int map_low = map + irf_low * ndet * nreg;
1318 int map_up = map + irf_up * ndet * nreg;
1319 for (int i = 0; i < npix; ++i) {
1320 irf(i, map) += wgt_low * m_irfs(i, map_low) +
1321 wgt_up * m_irfs(i, map_up);
1322 }
1323 }
1324 }
1325
1326 } // endfor: looped over integration intervals
1327
1328 // Debug: print weighting check
1329 #if defined(G_DEBUG_COMPUTE_IRF)
1330 std::cout << "- wgt_check = " << wgt_check << " (should be 1)" << std::endl;
1331 #endif
1332
1333 } // endif: integrated response over energy band
1334
1335 // Case B: compute response for line energy
1336 else {
1337
1338 // Get IRFs bracketing the minimum energy and the IRF weights
1339 m_energies.set_value(emin);
1340 int irf_low = m_energies.inx_left();
1341 int irf_up = m_energies.inx_right();
1342 double wgt_low = m_energies.wgt_left();
1343 double wgt_up = m_energies.wgt_right();
1344
1345 // Compute all IRF pixels
1346 for (int idet = 0; idet < ndet; ++idet) {
1347 for (int ireg = 0; ireg < nreg; ++ireg) {
1348 int map = idet + ireg * ndet;
1349 int map_low = map + irf_low * ndet * nreg;
1350 int map_up = map + irf_up * ndet * nreg;
1351 for (int i = 0; i < npix; ++i) {
1352 irf(i, map) = wgt_low * m_irfs(i, map_low) +
1353 wgt_up * m_irfs(i, map_up);
1354 }
1355 }
1356 }
1357
1358 // Debug: print weighting check
1359 #if defined(G_DEBUG_COMPUTE_IRF)
1360 std::cout << "- wgt_low = " << wgt_low << std::endl;
1361 std::cout << "- wgt_up = " << wgt_up << std::endl;
1362 std::cout << "- wgt_low+wgt_up = " << wgt_low+wgt_up;
1363 std::cout << " (should be 1)" << std::endl;
1364 #endif
1365
1366 } // endelse: compute response for line energy
1367
1368 // Return IRF
1369 return irf;
1370}
1371
1372
1373/***********************************************************************//**
1374 * @brief Set IRF image limits
1375 *
1376 * @param[in] image IRF FITS image.
1377 *
1378 * Computes the IRF image limits.
1379 ***************************************************************************/
1380void GSPIResponse::set_wcs(const GFitsImage* image)
1381{
1382 // Get image attributes
1383 int naxis1 = image->integer("NAXIS1");
1384 int naxis2 = image->integer("NAXIS2");
1385 double crval1 = image->real("CRVAL1");
1386 double crval2 = image->real("CRVAL2");
1387 double crpix1 = image->real("CRPIX1");
1388 double crpix2 = image->real("CRPIX2");
1389 double cdelt1 = image->real("CDELT1");
1390 double cdelt2 = image->real("CDELT2");
1391
1392 // Derive image limits. Limits are taken at the pixel centres since
1393 // we want to use them for bilinear interpolation. This means that
1394 // we will throw away half a pixel at the edge of the IRFs.
1395 m_wcs_xmin = (crval1 - (crpix1-1.0) * cdelt1) * gammalib::deg2rad;
1396 m_wcs_ymin = (crval2 - (crpix2-1.0) * cdelt2) * gammalib::deg2rad;
1397 m_wcs_xbin = cdelt1 * gammalib::deg2rad;
1398 m_wcs_ybin = cdelt2 * gammalib::deg2rad;
1399 m_wcs_xmax = m_wcs_xmin + double(naxis1-1) * m_wcs_xbin;
1400 m_wcs_ymax = m_wcs_ymin + double(naxis2-1) * m_wcs_ybin;
1401 m_wcs_xpix_max = double(naxis1-1);
1402 m_wcs_ypix_max = double(naxis2-1);
1403
1404 // Set maximum zenith angle
1405 m_max_zenith = (std::abs(m_wcs_xmax) > std::abs(m_wcs_ymax)) ?
1406 std::abs(m_wcs_xmax) : std::abs(m_wcs_ymax);
1407
1408 // Signal that image limits exist
1409 m_has_wcs = true;
1410
1411 // Return
1412 return;
1413}
1414
1415
1416/***********************************************************************//**
1417 * @brief Set vector of detector identifiers used by the observation
1418 *
1419 * @param[in] cube INTEGRAL/SPI event cube.
1420 *
1421 * Sets the vector of detector identifiers that is used by the observation.
1422 * The method scans the detector identifiers for all pointings and builds
1423 * a vector of unique detector identifiers in the order they were encountered
1424 * in the event cube.
1425 ***************************************************************************/
1426void GSPIResponse::set_detids(const GSPIEventCube* cube)
1427{
1428 // Clear vector of detector identifiers
1429 m_detids.clear();
1430
1431 // Extract relevant event cube dimensions
1432 int npt = cube->naxis(0);
1433 int ndet = cube->naxis(1);
1434
1435 // Loop over all pointings and detectors
1436 for (int ipt = 0; ipt < npt; ++ipt) {
1437 for (int idet = 0; idet < ndet; ++idet) {
1438
1439 // Get IRF detector identifier
1440 int detid = irf_detid(cube->dir(ipt, idet).detid());
1441
1442 // Push IRF detector identifier on vector if it does not yet exist
1443 std::vector<int>::iterator it = find(m_detids.begin(), m_detids.end(), detid);
1444 if (it == m_detids.end()) {
1445 m_detids.push_back(detid);
1446 }
1447
1448 } // endfor: looped over detectors
1449 } // endfor: looped over pointings
1450
1451 // Return
1452 return;
1453}
1454
1455
1456/***********************************************************************//**
1457 * @brief Set computation cache
1458 *
1459 * @param[in] cube INTEGRAL/SPI event cube.
1460 *
1461 * Setup of two vectors for fast coordinate transformation into the
1462 * instrument system. The first vector m_spix stores the SPI pointing
1463 * direction (the X direction) while the second vector stores the position
1464 * angle in celestial coordinates of the SPI Y direction.
1465 ***************************************************************************/
1466void GSPIResponse::set_cache(const GSPIEventCube* cube)
1467{
1468 // Extract relevant event cube dimensions
1469 int npt = cube->naxis(0);
1470
1471 // Clear vectors of SPI X direction and position angles
1472 m_spix.clear();
1473 m_spix.reserve(npt);
1474 m_posang.clear();
1475 m_posang.reserve(npt);
1476
1477 // Loop over all pointings
1478 for (int ipt = 0; ipt < npt; ++ipt) {
1479
1480 // Compute position angle in celestial coordinates
1481 double posang = cube->spi_x(ipt).posang(cube->spi_z(ipt)) + gammalib::pihalf;
1482
1483 // Store angles
1484 m_spix.push_back(cube->spi_x(ipt));
1485 m_posang.push_back(posang);
1486
1487 } // endfor: looped over pointings
1488
1489 // Return
1490 return;
1491}
1492
1493
1494/***********************************************************************//**
1495 * @brief Fill vector of INTEGRAL/SPI instrument response for IRF pixel
1496 *
1497 * @param[in] cube Pointer to event cube.
1498 * @param[in] xpix Response pixel in X direction.
1499 * @param[in] ypix Response pixel in Y direction.
1500 * @param[in] ireg IRF region (0 = photo peak).
1501 * @param[in] livetimes livetimes Pointer to livetime array (dimension m_num_det).
1502 * @param[out] irf Response vector.
1503 *
1504 * Computes a response vector for all detectors and energy bins for a given
1505 * response pixel defined by @p xpix and @p ypix, where
1506 *
1507 * xpix = (zenith * std::cos(azimuth) - m_wcs_xmin) / m_wcs_xbin
1508 * ypix = (zenith * std::sin(azimuth) - m_wcs_ymin) / m_wcs_ybin
1509 *
1510 * Checking of the validity of @p xpix and @p ypix needs to be done by the
1511 * client, makeing sure that
1512 *
1513 * xpix > 0.0 && xpix < m_wcs_xpix_max && ypix > 0.0 && ypix < m_wcs_ypix_max
1514 *
1515 * The client needs also to make sure that a vector with dimensions
1516 *
1517 * m_num_det*m_num_ebin
1518 *
1519 * is provided, and the method will fill all elements of the vector.
1520 ***************************************************************************/
1521void GSPIResponse::irf_vector(const GSPIEventCube* cube,
1522 const double& xpix,
1523 const double& ypix,
1524 const int& ireg,
1525 const double* livetimes,
1526 GVector* irf) const
1527{
1528 // Get number of pixels in X direction
1529 int nx = m_irfs.nx();
1530 int ndet = m_irfs.shape()[0];
1531 int nreg = m_irfs.shape()[1];
1532
1533 // Get 4 nearest neighbours
1534 int ix_left = int(xpix);
1535 int ix_right = ix_left + 1;
1536 int iy_top = int(ypix);
1537 int iy_bottom = iy_top + 1;
1538
1539 // Get weighting factors
1540 double wgt_right = xpix - double(ix_left);
1541 double wgt_left = 1.0 - wgt_right;
1542 double wgt_bottom = ypix - double(iy_top);
1543 double wgt_top = 1.0 - wgt_bottom;
1544
1545 // Get indices of 4 nearest neighbours
1546 int inx_1 = ix_left + iy_top * nx;
1547 int inx_2 = ix_right + iy_top * nx;
1548 int inx_3 = ix_left + iy_bottom * nx;
1549 int inx_4 = ix_right + iy_bottom * nx;
1550
1551 // Get weights of 4 nearest neighbours
1552 double wgt_1 = wgt_left * wgt_top;
1553 double wgt_2 = wgt_right * wgt_top;
1554 double wgt_3 = wgt_left * wgt_bottom;
1555 double wgt_4 = wgt_right * wgt_bottom;
1556
1557 // Initialise IRF value destination index
1558 int idst = 0;
1559
1560 // Loop over all detectors
1561 for (int idet = 0; idet < cube->m_num_det; ++idet) {
1562
1563 // Compute IRFs if livetime is positive
1564 if (livetimes[idet] > 0.0) {
1565
1566 // Get detid
1567 int detid = irf_detid(cube->m_detid[idet]);
1568
1569 // Loop over all energy bins
1570 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng, ++idst) {
1571
1572 // Get IRF map
1573 int map = detid + (ireg + ieng * nreg) * ndet;
1574
1575 // Compute IRF
1576 (*irf)[idst] = m_irfs(inx_1, map) * wgt_1 +
1577 m_irfs(inx_2, map) * wgt_2 +
1578 m_irfs(inx_3, map) * wgt_3 +
1579 m_irfs(inx_4, map) * wgt_4;
1580
1581 // Make sure that IRF is not negative
1582 if ((*irf)[idst] < 0.0) {
1583 (*irf)[idst] = 0.0;
1584 }
1585
1586 } // endfor: looped over energy bins
1587
1588 } // endif: livetime was positive
1589
1590 // ... otherwise set IRFs to zero
1591 else {
1592 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng, ++idst) {
1593 (*irf)[idst] = 0.0;
1594 }
1595 }
1596
1597 } // endfor: looped over detectors
1598
1599 // Return
1600 return;
1601}
1602
1603
1604/***********************************************************************//**
1605 * @brief Convert detector identifier into IRF detector identifier
1606 *
1607 * @param[in] detid SPI event detector identifier.
1608 * @return IRF detector identifier
1609 *
1610 * Converts a SPI event detector identifier into an IRF detector identifier.
1611 * TODO: Describe how this is done and why
1612 ***************************************************************************/
1613int GSPIResponse::irf_detid(const int& detid) const
1614{
1615 // Initialise irf detector identifier
1616 int irf_detid = detid;
1617
1618 // Put detector identifier in the relevant range
1619 if (irf_detid >= 123) {
1620 irf_detid -= 123;
1621 }
1622 else if (irf_detid >= 104) {
1623 irf_detid -= 104;
1624 }
1625 else if (irf_detid >= 85) {
1626 irf_detid -= 85;
1627 }
1628
1629 // Return irf detector identifier
1630 return irf_detid;
1631}
1632
1633
1634/***********************************************************************//**
1635 * @brief Compute weight of logarithmic energy bin
1636 *
1637 * @param[in] beta 1-gamma.
1638 * @param[in] emin Logarithmic minimum energy.
1639 * @param[in] emax Logarithmic maximum energy.
1640 ***************************************************************************/
1641double GSPIResponse::irf_weight(const double& beta,
1642 const double& emin,
1643 const double& emax) const
1644{
1645 // Initialise weight
1646 double weight;
1647
1648 // Assign weight
1649 if (std::abs(beta) < DBL_MIN) {
1650 weight = (emax - emin);
1651 }
1652 else {
1653 weight = std::exp(beta * gammalib::ln10 * emax) -
1654 std::exp(beta * gammalib::ln10 * emin);
1655 }
1656
1657 // Return weight
1658 return weight;
1659}
1660
1661
1662/***********************************************************************//**
1663 * @brief Return instrument response to point source sky model
1664 *
1665 * @param[in] model Sky model.
1666 * @param[in] obs Observation.
1667 * @param[in] gradients Gradients matrix.
1668 * @return Instrument response to point source sky model.
1669 *
1670 * Returns the instrument response to a point source sky model for all
1671 * events. The methods works on response vectors that are computed per
1672 * pointing. The method assumes that the spatial model component is of
1673 * type GModelSpatialPointSource and that the events in the observation
1674 * are of type GSPIEventCube. No checking of these assumptions is performed
1675 * by the method, and not respecting the assumptions will results in a
1676 * segfault.
1677 ***************************************************************************/
1678GVector GSPIResponse::irf_ptsrc(const GModelSky& model,
1679 const GObservation& obs,
1680 GMatrix* gradients) const
1681{
1682 // Get pointers to point source model spatial component
1683 const GModelSpatialPointSource* spatial = static_cast<const GModelSpatialPointSource*>(model.spatial());
1684 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
1685
1686 // Get number of events and IRF working vector values
1687 int nevents = obs.events()->size();
1688 int nirf = cube->m_num_det * cube->m_num_ebin;
1689
1690 // Initialise result and working vector
1691 GVector irfs(nevents);
1692 GVector wrk_irf(nirf);
1693
1694 // Get point source direction
1695 GSkyDir srcDir = spatial->dir();
1696
1697 // Loop over pointings
1698 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
1699
1700 // Convert sky direction to zenith angle. Continue only if zenith
1701 // angle is smaller than the maximum zenith angle.
1702 double zenith = this->zenith(ipt, srcDir);
1703 if (zenith < m_max_zenith) {
1704
1705 // Convert sky direction to azimuth angle
1706 double azimuth = this->azimuth(ipt, srcDir);
1707
1708 // Compute pixel
1709 double xpix = (zenith * std::cos(azimuth) - m_wcs_xmin) / m_wcs_xbin;
1710 double ypix = (zenith * std::sin(azimuth) - m_wcs_ymin) / m_wcs_ybin;
1711
1712 // Continue only if pixel is within IRF
1713 if (xpix > 0.0 && xpix < m_wcs_xpix_max &&
1714 ypix > 0.0 && ypix < m_wcs_ypix_max) {
1715
1716 // Get pointer to livetime array for all detectors
1717 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
1718
1719 // Compute IRF values for this pointing (photo peak only)
1720 irf_vector(cube, xpix, ypix, 0, livetimes, &wrk_irf);
1721
1722 // Copy IRF values into result vector
1723 for (int i = 0, idst = ipt * nirf; i < nirf; ++i, ++idst) {
1724 irfs[idst] = wrk_irf[i];
1725 }
1726
1727 } // endif: pixel was within IRF
1728
1729 } // endif: zenith angle was within validity range
1730
1731 } // endfor: looped over pointings
1732
1733 // Return IRF value
1734 return irfs;
1735}
1736
1737
1738/***********************************************************************//**
1739 * @brief Return instrument response to radial source sky model
1740 *
1741 * @param[in] model Sky model.
1742 * @param[in] obs Observation.
1743 * @param[in] gradients Gradients matrix.
1744 * @return Instrument response to radial source sky model.
1745 *
1746 * Returns the instrument response to a radial source sky model for all
1747 * events. The methods works on response vectors that are computed per
1748 * pointing. The method assumes that the spatial model component is of
1749 * type GModelSpatialRadial and that the events in the observation
1750 * are of type GSPIEventCube. No checking of these assumptions is performed
1751 * by the method, and not respecting the assumptions will results in a
1752 * segfault.
1753 *
1754 * The computation is performed by projecting for each pointing the IRF
1755 * pattern on the sky, and by evaluating the radial model at each IRF pixel.
1756 ***************************************************************************/
1757GVector GSPIResponse::irf_radial(const GModelSky& model,
1758 const GObservation& obs,
1759 GMatrix* gradients) const
1760{
1761 // Store entry time
1762 #if defined(G_DEBUG_IRF_TIMING)
1763 double cstart = gammalib::get_current_clock();
1764 #endif
1765
1766 #if defined(G_IRF_RADIAL_PROJECTION)
1767 #if defined(G_DEBUG_IRF_TIMING)
1768 std::cout << "GSPIResponse::irf_radial - projection" << std::endl;
1769 #endif
1770 // Set dummy energy and times
1771 static const GEnergy energy;
1772 static const GTime time;
1773
1774 // Get pointers to radial source model spatial component
1775 const GModelSpatialRadial* radial = static_cast<const GModelSpatialRadial*>(model.spatial());
1776 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
1777
1778 // Initialise response computation
1779 GSPIResponseIrf irf;
1780 irf_init(obs, &irf);
1781
1782 // Initialise result and working vector
1783 GVector irfs(irf.nevents);
1784
1785 // Loop over pointings
1786 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
1787
1788 // Compute distance between SPI pointing direction and centre
1789 // of radial model in radians
1790 double centre_distance = this->zenith(ipt, radial->dir());
1791
1792 // Compute IRF in case that radial model overlaps with SPI field
1793 // of view
1794 if (centre_distance < (m_max_zenith + radial->theta_max())) {
1795
1796 // Get pointer to livetime array for all detectors
1797 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
1798
1799 // Set rotation matrix for current pointing to transform from IRF
1800 // coordinate system to the celestial coordinate system
1801 irf_rot_matrix(ipt, &irf);
1802
1803 // Loop over IRF pixels
1804 for (int inx = 0; inx < irf.nxy; ++inx) {
1805
1806 // Skip zero IRF pixels
1807 if (irf.zero[inx]) {
1808 continue;
1809 }
1810
1811 // Set sky direction for IRF pixel
1812 irf_skydir(ipt, inx, &irf);
1813
1814 // Compute offset angle to centre of radial model
1815 double theta = radial->dir().dist(irf.dir);
1816
1817 // Continue only if radial offset angle value is valid
1818 if (theta < radial->theta_max()) {
1819
1820 // Compute radial model value
1821 double model = radial->eval(theta, energy, time);
1822
1823 // Continue only if model value is positive
1824 if (model > 0.0) {
1825
1826 // Multiply model value by solid angle of IRF pixel
1827 model *= irf.solidangle[inx];
1828
1829 // Loop over all detectors
1830 for (int idet = 0, idst = ipt * irf.nirf; idet < cube->m_num_det; ++idet) {
1831
1832 // Compute IRFs if livetime is positive
1833 if (livetimes[idet] > 0.0) {
1834
1835 // Loop over all energy bins
1836 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng, ++idst) {
1837
1838 // Get IRF map
1839 int map = irf_map(idet, ieng, &irf);
1840
1841 // Add IRF
1842 irfs[idst] += m_irfs(inx, map) * model;
1843
1844 } // endfor: looped over energy bins
1845
1846 } // endif: livetime was positive
1847
1848 // ... otherwise simply increment idst index
1849 else {
1850 idst += cube->m_num_ebin;
1851 }
1852
1853 } // endfor: looped over detectors
1854
1855 } // endif: model was positive
1856
1857 } // endif: theta was valid
1858
1859 } // endfor: looped over IRF pixels
1860
1861 } // endif: model was in field of view
1862
1863 } // endfor: looped over pointings
1864 #else
1865 #if defined(G_DEBUG_IRF_TIMING)
1866 std::cout << "GSPIResponse::irf_radial - numerical integration" << std::endl;
1867 #endif
1868 // Set constants
1869 static const int iter_rho = 6;
1870 static const int iter_omega = 6;
1871
1872 // Get pointers to radial spatial component
1873 const GModelSpatialRadial* radial = static_cast<const GModelSpatialRadial*>(model.spatial());
1874 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
1875
1876 // Get number of events and IRF working vector values
1877 int nevents = obs.events()->size();
1878 int nirf = cube->m_num_det * cube->m_num_ebin;
1879
1880 // Initialise result and working vector
1881 GVector irfs(nevents);
1882 GVector wrk_irfs(nirf);
1883
1884 // Set radial integration range (radians)
1885 double rho_min = 0.0;
1886 double rho_max = radial->theta_max() - 1.0e-12; // Kludge to stay inside
1887 // the boundary of a
1888 // sharp edged model
1889
1890 // Integrate only if radial integration interval is valid
1891 if (rho_min < rho_max) {
1892
1893 // Allocate Euler and rotation matrices
1894 GMatrix ry;
1895 GMatrix rz;
1896
1897 // Set up rotation matrix to rotate from native model coordinates to
1898 // celestial coordinates
1899 ry.eulery( radial->dir().dec_deg() - 90.0);
1900 rz.eulerz(-radial->dir().ra_deg());
1901 GMatrix rot = (ry * rz).transpose();
1902
1903 // Loop over pointings
1904 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
1905
1906 // Compute distance between SPI pointing direction and centre
1907 // of radial model in radians
1908 double centre_distance = this->zenith(ipt, radial->dir());
1909
1910 // Compute IRF in case that radial model overlaps with SPI field
1911 // of view
1912 if (centre_distance < (m_max_zenith + rho_max)) {
1913
1914 // Get pointer to livetime array for all detectors
1915 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
1916
1917
1918 // Setup integration kernel
1919 spi_radial_kerns_rho integrands(cube,
1920 this,
1921 radial,
1922 ipt,
1923 livetimes,
1924 rot,
1925 iter_omega,
1926 wrk_irfs);
1927
1928 // Integrate over model's radial coordinate
1929 GIntegrals integral(&integrands);
1930 integral.fixed_iter(iter_rho);
1931
1932 // Integrate over model
1933 wrk_irfs = integral.romberg(rho_min, rho_max, iter_rho);
1934
1935 // Copy IRF values into result vector
1936 for (int i = 0, idst = ipt * nirf; i < nirf; ++i, ++idst) {
1937 irfs[idst] = wrk_irfs[i];
1938 }
1939
1940 } // endif: Radial model overlapped with SPI field of view
1941
1942 } // endfor: looped over pointings
1943
1944 } // endif: integration range was valid
1945 #endif
1946
1947 // Dump CPU consumption
1948 #if defined(G_DEBUG_IRF_TIMING)
1949 double celapse = gammalib::get_current_clock() - cstart;
1950 std::cout << "GSPIResponse::irf_radial consumed " << celapse;
1951 std::cout << " seconds of CPU time." << std::endl;
1952 #endif
1953
1954 // Return IRF value
1955 return irfs;
1956}
1957
1958
1959/***********************************************************************//**
1960 * @brief Return instrument response to elliptical source sky model
1961 *
1962 * @param[in] model Sky model.
1963 * @param[in] obs Observation.
1964 * @param[in] gradients Gradients matrix.
1965 * @return Instrument response to ellipitical source sky model.
1966 *
1967 * Returns the instrument response to a elliptical source sky model for all
1968 * events. The methods works on response vectors that are computed per
1969 * pointing. The method assumes that the spatial model component is of
1970 * type GModelSpatialElliptical and that the events in the observation
1971 * are of type GSPIEventCube. No checking of these assumptions is performed
1972 * by the method, and not respecting the assumptions will results in a
1973 * segfault.
1974 *
1975 * The computation is performed by projecting for each pointing the IRF
1976 * pattern on the sky, and by evaluating the elliptical model at each IRF
1977 * pixel.
1978 ***************************************************************************/
1979GVector GSPIResponse::irf_elliptical(const GModelSky& model,
1980 const GObservation& obs,
1981 GMatrix* gradients) const
1982{
1983 // Store entry time
1984 #if defined(G_DEBUG_IRF_TIMING)
1985 double cstart = gammalib::get_current_clock();
1986 #endif
1987
1988 #if defined(G_IRF_ELLIPTICAL_PROJECTION)
1989 #if defined(G_DEBUG_IRF_TIMING)
1990 std::cout << "GSPIResponse::irf_elliptical - projection" << std::endl;
1991 #endif
1992 // Set dummy energy and times
1993 static const GEnergy energy;
1994 static const GTime time;
1995
1996 // Get pointers to elliptical source model spatial component
1997 const GModelSpatialElliptical* elliptical = static_cast<const GModelSpatialElliptical*>(model.spatial());
1998 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
1999
2000 // Initialise response computation
2001 GSPIResponseIrf irf;
2002 irf_init(obs, &irf);
2003
2004 // Initialise result and working vector
2005 GVector irfs(irf.nevents);
2006
2007 // Loop over pointings
2008 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
2009
2010 // Compute distance between SPI pointing direction and centre
2011 // of radial model in radians
2012 double centre_distance = this->zenith(ipt, elliptical->dir());
2013
2014 // Compute IRF in case that radial model overlaps with SPI field
2015 // of view
2016 if (centre_distance < (m_max_zenith + elliptical->theta_max())) {
2017
2018 // Get pointer to livetime array for all detectors
2019 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
2020
2021 // Set rotation matrix for current pointing to transform from IRF
2022 // coordinate system to the celestial coordinate system
2023 irf_rot_matrix(ipt, &irf);
2024
2025 // Loop over IRF pixels
2026 for (int inx = 0; inx < irf.nxy; ++inx) {
2027
2028 // Skip zero IRF pixels
2029 if (irf.zero[inx]) {
2030 continue;
2031 }
2032
2033 // Set sky direction for IRF pixel
2034 irf_skydir(ipt, inx, &irf);
2035
2036 // Compute offset angle to centre of elliptical model
2037 double theta = elliptical->dir().dist(irf.dir);
2038
2039 // Continue only if offset angle value is valid
2040 if (theta < elliptical->theta_max()) {
2041
2042 // Compute position angle of elliptical model in the
2043 // coordinate system of the elliptical model
2044 double posang = elliptical->dir().posang(irf.dir, elliptical->coordsys());
2045
2046 // Compute elliptical model value
2047 double model = elliptical->eval(theta, posang, energy, time);
2048
2049 // Continue only if model value is positive
2050 if (model > 0.0) {
2051
2052 // Multiply model value by solid angle of IRF pixel
2053 model *= irf.solidangle[inx];
2054
2055 // Loop over all detectors
2056 for (int idet = 0, idst = ipt * irf.nirf; idet < cube->m_num_det; ++idet) {
2057
2058 // Compute IRFs if livetime is positive
2059 if (livetimes[idet] > 0.0) {
2060
2061 // Loop over all energy bins
2062 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng, ++idst) {
2063
2064 // Get IRF map
2065 int map = irf_map(idet, ieng, &irf);
2066
2067 // Add IRF
2068 irfs[idst] += m_irfs(inx, map) * model;
2069
2070 } // endfor: looped over energy bins
2071
2072 } // endif: livetime was positive
2073
2074 // ... otherwise simply increment idst index
2075 else {
2076 idst += cube->m_num_ebin;
2077 }
2078
2079 } // endfor: looped over detectors
2080
2081 } // endif: model was positive
2082
2083 } // endif: theta was valid
2084
2085 } // endfor: looped over IRF pixels
2086
2087 } // endif: model was in field of view
2088
2089 } // endfor: looped over pointings
2090 #else
2091 #if defined(G_DEBUG_IRF_TIMING)
2092 std::cout << "GSPIResponse::irf_elliptical - numerical integration" << std::endl;
2093 #endif
2094 // Set constants
2095 static const int iter_rho = 6;
2096 static const int iter_omega = 6;
2097
2098 // Get pointers to radial spatial component
2099 const GModelSpatialElliptical* elliptical = static_cast<const GModelSpatialElliptical*>(model.spatial());
2100 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
2101
2102 // Get number of events and IRF working vector values
2103 int nevents = obs.events()->size();
2104 int nirf = cube->m_num_det * cube->m_num_ebin;
2105
2106 // Initialise result and working vector
2107 GVector irfs(nevents);
2108 GVector wrk_irfs(nirf);
2109
2110 // Set radial integration range (radians)
2111 double rho_min = 0.0;
2112 double rho_max = elliptical->theta_max() - 1.0e-6; // Kludge to stay inside
2113 // the boundary of a
2114 // sharp edged model
2115
2116 // Integrate only if radial integration interval is valid
2117 if (rho_min < rho_max) {
2118
2119 // Allocate Euler and rotation matrices
2120 GMatrix ry;
2121 GMatrix rz;
2122
2123 // Set up rotation matrix to rotate from native model coordinates to
2124 // celestial coordinates
2125 ry.eulery( elliptical->dir().dec_deg() - 90.0);
2126 rz.eulerz(-elliptical->dir().ra_deg());
2127 GMatrix rot = (ry * rz).transpose();
2128
2129 // Loop over pointings
2130 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
2131
2132 // Compute distance between SPI pointing direction and centre
2133 // of radial model in radians
2134 double centre_distance = this->zenith(ipt, elliptical->dir());
2135
2136 // Compute IRF in case that radial model overlaps with SPI field
2137 // of view
2138 if (centre_distance < (m_max_zenith + rho_max)) {
2139
2140 // Get pointer to livetime array for all detectors
2141 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
2142
2143 // Setup integration kernel
2144 spi_elliptical_kerns_rho integrands(cube,
2145 this,
2146 elliptical,
2147 ipt,
2148 livetimes,
2149 rot,
2150 iter_omega,
2151 wrk_irfs);
2152
2153 // Integrate over model's radial coordinate
2154 GIntegrals integral(&integrands);
2155 integral.fixed_iter(iter_rho);
2156
2157 // Integrate over model
2158 wrk_irfs = integral.romberg(rho_min, rho_max, iter_rho);
2159
2160 // Copy IRF values into result vector
2161 for (int i = 0, idst = ipt * nirf; i < nirf; ++i, ++idst) {
2162 irfs[idst] = wrk_irfs[i];
2163 }
2164
2165 } // endif: Radial model overlapped with SPI field of view
2166
2167 } // endfor: looped over pointings
2168
2169 } // endif: integration range was valid
2170 #endif
2171
2172 // Dump CPU consumption
2173 #if defined(G_DEBUG_IRF_TIMING)
2174 double celapse = gammalib::get_current_clock() - cstart;
2175 std::cout << "GSPIResponse::irf_elliptical consumed " << celapse;
2176 std::cout << " seconds of CPU time." << std::endl;
2177 #endif
2178
2179 // Return IRF value
2180 return irfs;
2181}
2182
2183
2184/***********************************************************************//**
2185 * @brief Return instrument response to diffuse source sky model
2186 *
2187 * @param[in] model Sky model.
2188 * @param[in] obs Observation.
2189 * @param[in] gradients Gradients matrix.
2190 * @return Instrument response to diffuse source sky model.
2191 *
2192 * Returns the instrument response to a diffuse source sky model for all
2193 * events. The methods works on response vectors that are computed per
2194 * pointing. The method assumes that the spatial model component is of
2195 * type GModelSpatialDiffuse and that the events in the observation
2196 * are of type GSPIEventCube. No checking of these assumptions is performed
2197 * by the method, and not respecting the assumptions will results in a
2198 * segfault.
2199 *
2200 * The computation is performed by projecting for each pointing the IRF
2201 * pattern on the sky, and by evaluating the diffuse model at each IRF pixel.
2202 ***************************************************************************/
2203GVector GSPIResponse::irf_diffuse(const GModelSky& model,
2204 const GObservation& obs,
2205 GMatrix* gradients) const
2206{
2207 // Store entry time
2208 #if defined(G_DEBUG_IRF_TIMING)
2209 double cstart = gammalib::get_current_clock();
2210 #endif
2211
2212 #if defined(G_IRF_DIFFUSE_PROJECTION)
2213 #if defined(G_DEBUG_IRF_TIMING)
2214 std::cout << "GSPIResponse::irf_diffuse - projection" << std::endl;
2215 #endif
2216 // Set dummy time
2217 static const GTime time;
2218
2219 // Get pointers to diffuse source model spatial component
2220 const GModelSpatialDiffuse* diffuse = static_cast<const GModelSpatialDiffuse*>(model.spatial());
2221 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
2222
2223 // Initialise response computation
2224 GSPIResponseIrf irf;
2225 irf_init(obs, &irf);
2226
2227 // Initialise result and working vector
2228 GVector irfs(irf.nevents);
2229
2230 // Loop over pointings
2231 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
2232
2233 // Get pointer to livetime array for all detectors
2234 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
2235
2236 // Set rotation matrix for current pointing to transform from IRF
2237 // coordinate system to the celestial coordinate system
2238 irf_rot_matrix(ipt, &irf);
2239
2240 // Loop over IRF pixels
2241 for (int inx = 0; inx < irf.nxy; ++inx) {
2242
2243 // Skip zero IRF pixels
2244 if (irf.zero[inx]) {
2245 continue;
2246 }
2247
2248 // Set sky direction for IRF pixel
2249 irf_skydir(ipt, inx, &irf);
2250
2251 // Get destination index offset
2252 int idst0 = ipt * irf.nirf;
2253
2254 // Loop over all energy bins
2255 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng) {
2256
2257 // Setup photon
2258 GPhoton photon(irf.dir, cube->m_energy[ieng], time);
2259
2260 // Compute model value
2261 double model = diffuse->eval(photon);
2262
2263 // Continue only if model is positive
2264 if (model > 0.0) {
2265
2266 // Multiply model value by solid angle of IRF pixel
2267 model *= irf.solidangle[inx];
2268
2269 // Loop over all detectors
2270 for (int idet = 0, idst = idst0 + ieng; idet < cube->m_num_det; ++idet) {
2271
2272 // Compute IRFs if livetime is positive
2273 if (livetimes[idet] > 0.0) {
2274
2275 // Get IRF map
2276 int map = irf_map(idet, ieng, &irf);
2277
2278 // Add IRF
2279 irfs[idst] += m_irfs(inx, map) * model;
2280
2281 } // endif: livetime was positive
2282
2283 // Increment destination index
2284 idst += cube->m_num_ebin;
2285
2286 } // endfor: looped over detectors
2287
2288 } // endif: model was positive
2289
2290 } // endfor: looped over energy bins
2291
2292 } // endfor: looped over IRF pixels
2293
2294 } // endfor: looped over pointings
2295 #else
2296 #if defined(G_DEBUG_IRF_TIMING)
2297 std::cout << "GSPIResponse::irf_diffuse - not projection" << std::endl;
2298 #endif
2299 // Get pointer to SPI event cube
2300 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
2301
2302 // Get number of events
2303 int nevents = obs.events()->size();
2304
2305 // Initialise result vector
2306 GVector irfs(nevents);
2307
2308 // If we have a diffuse map sky model then branch to dedicated method
2309 const GModelSpatialDiffuseMap* map = dynamic_cast<const GModelSpatialDiffuseMap*>(model.spatial());
2310 if (map != NULL) {
2311 irf_diffuse_map(map, obs, &irfs, gradients);
2312 }
2313
2314 // ... otherwise use bin-wise computation
2315 else {
2316
2317 // Loop over events
2318 for (int k = 0; k < nevents; ++k) {
2319
2320 // Get reference to event
2321 const GEvent& event = *((*obs.events())[k]);
2322
2323 // Get source energy and time (no dispersion)
2324 GEnergy srcEng = event.energy();
2325 GTime srcTime = event.time();
2326
2327 // Setup source
2328 GSource source(model.name(), model.spatial(), srcEng, srcTime);
2329
2330 // Get IRF value for event
2331 irfs[k] = GResponse::irf_diffuse(event, source, obs);
2332
2333 } // endfor: looped over events
2334
2335 } // endelse: used bin-wise computation
2336 #endif
2337
2338 // Dump CPU consumption
2339 #if defined(G_DEBUG_IRF_TIMING)
2340 double celapse = gammalib::get_current_clock() - cstart;
2341 std::cout << "GSPIResponse::irf_diffuse consumed " << celapse;
2342 std::cout << " seconds of CPU time." << std::endl;
2343 #endif
2344
2345 // Return IRF value
2346 return irfs;
2347}
2348
2349
2350/***********************************************************************//**
2351 * @brief Return instrument response to diffuse map sky model
2352 *
2353 * @param[in] diffuse Pointer to diffuse map spatial model.
2354 * @param[in] obs Observation.
2355 * @param[out] irfs Pointer to IRF vector.
2356 * @param[in] gradients Gradients matrix.
2357 *
2358 * Returns the instrument response to a diffuse source sky model for all
2359 * events. The methods works on response vectors that are computed per
2360 * pointing. The method assumes that the spatial model component is of
2361 * type GModelSpatialDiffuseMap and that the events in the observation
2362 * are of type GSPIEventCube. No checking of these assumptions is performed
2363 * by the method, and not respecting the assumptions will results in a
2364 * segfault.
2365 ***************************************************************************/
2366void GSPIResponse::irf_diffuse_map(const GModelSpatialDiffuseMap* diffuse,
2367 const GObservation& obs,
2368 GVector* irfs,
2369 GMatrix* gradients) const
2370{
2371 // Get pointer to SPI event cube
2372 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
2373
2374 // Get reference to diffuse map
2375 const GSkyMap& map = diffuse->map();
2376
2377 // Get number of events and IRF working vector values
2378 int nevents = obs.events()->size();
2379 int nirf = cube->m_num_det * cube->m_num_ebin;
2380 int nsky = map.npix();
2381
2382 // Initialise working vector
2383 GVector wrk_irf(nirf);
2384
2385 // Initialise vector of sky directions
2386 std::vector<GSkyDir> mapDirs(nsky);
2387 for (int isky = 0; isky < nsky; ++isky) {
2388 mapDirs[isky] = map.inx2dir(isky);
2389 }
2390
2391 // Loop over pointings
2392 for (int ipt = 0; ipt < cube->m_num_pt; ++ipt) {
2393
2394 // Loop over sky map pixels
2395 for (int isky = 0; isky < nsky; ++isky) {
2396
2397 // Get sky map value
2398 double value = map(isky);
2399
2400 // Go to next pixel if map value is zero
2401 if (value == 0.0) {
2402 continue;
2403 }
2404
2405 // Compute zenith angle for sky map pixel direction
2406 double zenith = this->zenith(ipt, mapDirs[isky]);
2407
2408 // Continue only if zenith angle is smaller than the maximum
2409 // zenith angle (checking whether sky map pixel is within
2410 // field of view)
2411 if (zenith < m_max_zenith) {
2412
2413 // Convert sky direction to azimuth angle
2414 double azimuth = this->azimuth(ipt, mapDirs[isky]);
2415
2416 // Compute pixel
2417 double xpix = (zenith * std::cos(azimuth) - m_wcs_xmin) / m_wcs_xbin;
2418 double ypix = (zenith * std::sin(azimuth) - m_wcs_ymin) / m_wcs_ybin;
2419
2420 // Continue only if pixel is within IRF
2421 if (xpix > 0.0 && xpix < m_wcs_xpix_max &&
2422 ypix > 0.0 && ypix < m_wcs_ypix_max) {
2423
2424 // Multiply sky map value with solid angle
2425 value *= map.solidangle(isky);
2426
2427 // Get pointer to livetime array for all detectors
2428 const double* livetimes = cube->m_livetime + ipt * cube->m_num_det;
2429
2430 // Compute IRF values for this pointing (photo peak only)
2431 irf_vector(cube, xpix, ypix, 0, livetimes, &wrk_irf);
2432
2433 // Add IRF values multiplied by sky map value into result vector
2434 for (int i = 0, idst = ipt * nirf; i < nirf; ++i, ++idst) {
2435 (*irfs)[idst] += wrk_irf[i] * value;
2436 }
2437
2438 } // endif: pixel was within IRF
2439
2440 } // endif: zenith angle was within validity range
2441
2442 } // endfor: looped over sky map pixels
2443
2444 } // endfor: looped over pointings
2445
2446 // Return
2447 return;
2448}
2449
2450
2451/***********************************************************************//**
2452 * @brief Initialise structure for projection computation of response
2453 *
2454 * @param[in] obs Observation.
2455 * @param[out] irf Stucture for projection computation of response.
2456 *
2457 * The method initialises the @p irf structure that is needed for the
2458 * projection computations of the response. The method initialises the
2459 * following members of the structure:
2460 *
2461 * nevents : Number of events in observation
2462 * nirf : Number of IRFs in observation
2463 * nx : Number of pixels in x direction of IRF
2464 * ny : Number of pixels in y direction of IRF
2465 * nxy : Number of pixels of IRF
2466 * ndet : Number of detectors in IRF
2467 * nreg : Number of regions in IRF
2468 * ireg : Used IRF region (set to 0, i.e. use photo peak response)
2469 * detids : Vector of detector indices
2470 * zero : Vector of booleans flagging IRF pixels that are zero
2471 * zeniths : Zenith angles for all IRF pixels (radians)
2472 * sin_zeniths : Sine of zenith angles for all IRF pixels
2473 * cos_zeniths : Cosine of zenith angles for all IRF pixels
2474 * azimuths : Azimuth angles for all IRF pixels (radians)
2475 * solidangle : Solid angles for all IRF pixels
2476 ***************************************************************************/
2477void GSPIResponse::irf_init(const GObservation& obs,
2478 GSPIResponseIrf* irf) const
2479{
2480 // Get pointers to event cube
2481 const GSPIEventCube* cube = static_cast<const GSPIEventCube*>(obs.events());
2482
2483 // Get number of events and IRF working vector values
2484 irf->nevents = obs.events()->size();
2485 irf->nirf = cube->m_num_det * cube->m_num_ebin;
2486 irf->nx = m_irfs.nx();
2487 irf->ny = m_irfs.ny();
2488 irf->nxy = irf->nx * irf->ny;
2489 irf->ndet = m_irfs.shape()[0];
2490 irf->nreg = m_irfs.shape()[1];
2491 irf->ireg = 0; // Photo peak
2492
2493 // Pre-compute detector indices
2494 irf->detids = std::vector<int>(cube->m_num_det);
2495 for (int idet = 0; idet < cube->m_num_det; ++idet) {
2496 irf->detids[idet] = irf_detid(cube->m_detid[idet]);
2497 }
2498
2499 // Pre-compute boolean vector that flags all zero IRF pixels
2500 irf->zero = std::vector<bool>(irf->nxy, true);
2501 for (int inx = 0; inx < irf->nxy; ++inx) {
2502 for (int idet = 0; idet < cube->m_num_det; ++idet) {
2503 for (int ieng = 0; ieng < cube->m_num_ebin; ++ieng) {
2504 int map = irf->detids[idet] + (irf->ireg + ieng * irf->nreg) * irf->ndet;
2505 if (m_irfs(inx, map) > 0.0) {
2506 irf->zero[inx] = false;
2507 break;
2508 }
2509 }
2510 if (!irf->zero[inx]) {
2511 break;
2512 }
2513 }
2514 }
2515
2516 // Initilise IRF related vectors for faster computation
2517 irf->zeniths = GVector(irf->nxy);
2518 irf->sin_zeniths = GVector(irf->nxy);
2519 irf->cos_zeniths = GVector(irf->nxy);
2520 irf->azimuths = GVector(irf->nxy);
2521 irf->solidangle = GVector(irf->nxy);
2522
2523 // Loop over y pixels
2524 for (int iy = 0, inx = 0; iy < irf->ny; ++iy) {
2525
2526 // Compute y pixel and its square
2527 double ypix = m_wcs_ymin + iy * m_wcs_ybin;
2528 double ypix2 = ypix * ypix;
2529
2530 // Loop over x pixels
2531 for (int ix = 0; ix < irf->nx; ++ix, ++inx) {
2532
2533 // Compute x pixel
2534 double xpix = m_wcs_xmin + ix * m_wcs_xbin;
2535
2536 // Pre-compute vector values
2537 irf->zeniths[inx] = std::sqrt(xpix*xpix + ypix2);
2538 irf->azimuths[inx] = std::atan2(ypix, xpix);
2539 irf->sin_zeniths[inx] = std::sin(irf->zeniths[inx]);
2540 irf->cos_zeniths[inx] = std::cos(irf->zeniths[inx]);
2541 irf->solidangle[inx] = m_irfs.solidangle(inx);
2542
2543 } // endfor: looped over x pixels
2544
2545 } // endfor: looped over y pixels
2546
2547 // Return
2548 return;
2549}
2550
2551
2552/***********************************************************************//**
2553 * @brief Set rotation matrix for pointing
2554 *
2555 * @param[in] ipt Pointing index.
2556 * @param[in,out] irf Stucture for projection computation of response.
2557 *
2558 * Sets the rotation matrix for the pointing with index @p ipt in the
2559 * @p irf structure that allows rotation from the IRF coordinate system
2560 * to the celestial coordinate system. The method sets the following
2561 * member of the @p irf structure:
2562 *
2563 * rot : Rotation matrix
2564 *
2565 * The client needs to make sure that the pointing index is comprised within
2566 * the valid range and that the GSPIResponse::irf_init() method was called
2567 * prior to invoking this method.
2568 ***************************************************************************/
2569void GSPIResponse::irf_rot_matrix(const int& ipt,
2570 GSPIResponseIrf* irf) const
2571{
2572 // Allocate Euler matrices
2573 GMatrix ry;
2574 GMatrix rz;
2575
2576 // Set Euler angles
2577 ry.eulery(m_spix[ipt].dec_deg() - 90.0);
2578 rz.eulerz(-m_spix[ipt].ra_deg());
2579
2580 // Compute rotation matrix
2581 irf->rot = (ry * rz).transpose();
2582
2583 // Return
2584 return;
2585}
2586
2587
2588/***********************************************************************//**
2589* @brief Set sky direction for pointing and IRF pixel index
2590*
2591* @param[in] ipt Pointing index.
2592* @param[in] inx IRF pixel index.
2593* @param[in,out] irf Stucture for projection computation of response.
2594*
2595* Set the sky direction of the IRF pixel index @p inx for a given pointing
2596* index @p ipt. The method makes use of the rotation matrix that needs
2597* to be computed before by the GSPIResponse::irf_rot_matrix() method. The
2598* method sets the following member of the @p irf structure:
2599*
2600* dir : Sky direction
2601*
2602* The client needs to make sure that the pointing index and thr IRF pixel
2603* index is comprised within the valid ranges and that the
2604* GSPIResponse::irf_init() and GSPIResponse::irf_rot_matrix() methods were
2605* called prior to invoking this method.
2606***************************************************************************/
2607void GSPIResponse::irf_skydir(const int& ipt,
2608 const int& inx,
2609 GSPIResponseIrf* irf) const
2610{
2611 // Compute azimuth angle (radians)
2612 double azimuth = m_posang[ipt] - irf->azimuths[inx];
2613
2614 // Compute cosine and sine
2615 double cos_phi = std::cos(azimuth);
2616 double sin_phi = std::sin(azimuth);
2617
2618 // Get vector or IRF pixel
2619 GVector native(-cos_phi * irf->sin_zeniths[inx],
2620 sin_phi * irf->sin_zeniths[inx],
2621 irf->cos_zeniths[inx]);
2622
2623 // Rotate IRF vector into sky system
2624 GVector vector = irf->rot * native;
2625
2626 // Transform IRF vector into sky direction
2627 irf->dir.celvector(vector);
2628
2629 // Debug: check coordinate transformations
2630 #if defined(G_DEBUG_IRF_SKYDIR)
2631 double zenith_check = this->zenith(ipt, irf->dir);
2632 double azimuth_check = this->azimuth(ipt, irf->dir);
2633 double xpix_check = (zenith_check * std::cos(azimuth_check) - m_wcs_xmin) / m_wcs_xbin;
2634 double ypix_check = (zenith_check * std::sin(azimuth_check) - m_wcs_ymin) / m_wcs_ybin;
2635 int ix = inx % irf->nx;
2636 int iy = inx / irf->nx;
2637 if ((std::abs(double(ix)-xpix_check) > 1.0e-3) ||
2638 (std::abs(double(iy)-ypix_check) > 1.0e-3)) {
2639 std::cout << "*** Bad coordinate transformation: (";
2640 std::cout << ix << "," << iy << ") differs from (";
2641 std::cout << xpix_check << "," << ypix_check << ")" << std::endl;
2642 }
2643 #endif
2644
2645 // Return
2646 return;
2647}
G_IRF
#define G_IRF
Definition GCOMResponse.cpp:63
G_EBOUNDS
#define G_EBOUNDS
Definition GCOMResponse.cpp:68
G_NROI
#define G_NROI
Definition GCOMResponse.cpp:66
GEbounds.hpp
Energy boundaries class interface definition.
GEnergies.hpp
Energy container class definition.
GException.hpp
Exception handler interface definition.
GFitsBinTable.hpp
FITS binary table class definition.
GFitsTableDoubleCol.hpp
FITS table double column class interface definition.
GFitsTableShortCol.hpp
FITS table short integer column class interface definition.
GFits.hpp
FITS file class interface definition.
GIntegrals.hpp
Integration class for set of functions interface definition.
GMatrix.hpp
Generic matrix class definition.
GModelSky.hpp
Sky model class interface definition.
GModelSpatialDiffuseMap.hpp
Spatial map model class interface definition.
GModelSpatialDiffuse.hpp
Abstract diffuse spatial model base class interface definition.
GModelSpatialElliptical.hpp
Abstract elliptical spatial model base class interface definition.
GModelSpatialPointSource.hpp
Point source spatial model class interface definition.
GModelSpatialRadial.hpp
Abstract radial spatial model base class interface definition.
GNodeArray.hpp
Node array class interface definition.
GSPIEventBin.hpp
INTEGRAL/SPI event bin class definition.
GSPIEventCube.hpp
INTEGRAL/SPI event bin container class definition.
GSPIObservation.hpp
INTEGRAL/SPI observation class definition.
G_LOAD_IRF
#define G_LOAD_IRF
Definition GSPIResponse.cpp:67
G_LOAD_IRFS
#define G_LOAD_IRFS
Definition GSPIResponse.cpp:68
GSPIResponse.hpp
INTEGRAL/SPI instrument response function class definition.
GSPITools.hpp
Definition of INTEGRAL/SPI tools.
GSource.hpp
Source class definition.
GTools.hpp
Gammalib tools definition.
GChatter
GChatter
Definition GTypemaps.hpp:33
SILENT
@ SILENT
Definition GTypemaps.hpp:34
GVector.hpp
Vector class interface definition.
GEbounds
Energy boundaries container class.
Definition GEbounds.hpp:60
GEbounds::read
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition GEbounds.cpp:761
GEbounds::append
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition GEbounds.cpp:301
GEbounds::clear
void clear(void)
Clear energy boundaries.
Definition GEbounds.cpp:268
GEbounds::is_empty
bool is_empty(void) const
Signal if there are no energy boundaries.
Definition GEbounds.hpp:175
GEbounds::write
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
GEnergy
Class that handles energies in a unit independent way.
Definition GEnergy.hpp:48
GEnergy::keV
double keV(void) const
Return energy in keV.
Definition GEnergy.cpp:327
GEvent
Abstract interface for the event classes.
Definition GEvent.hpp:71
GEvents::ebounds
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition GEvents.cpp:136
GEvents::size
virtual int size(void) const =0
GException::feature_not_implemented
Definition GException.hpp:144
GException::invalid_argument
Definition GException.hpp:88
GException::invalid_value
Definition GException.hpp:81
GFilename
Filename class.
Definition GFilename.hpp:62
GFilename::path
std::string path(void) const
Return access path.
Definition GFilename.hpp:217
GFilename::url
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition GFilename.hpp:189
GFilename::clear
void clear(void)
Clear file name.
Definition GFilename.cpp:188
GFitsBinTable
FITS binary table class.
Definition GFitsBinTable.hpp:39
GFitsHDU::has_card
bool has_card(const int &cardno) const
Check existence of header card.
Definition GFitsHDU.hpp:233
GFitsHDU::extname
const std::string & extname(void) const
Return extension name.
Definition GFitsHDU.hpp:162
GFitsHDU::real
double real(const std::string &keyname) const
Return card value as double precision.
Definition GFitsHDU.hpp:423
GFitsHDU::string
std::string string(const std::string &keyname) const
Return card value as string.
Definition GFitsHDU.hpp:410
GFitsHDU::card
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition GFitsHDU.hpp:259
GFitsHDU::integer
int integer(const std::string &keyname) const
Return card value as integer.
Definition GFitsHDU.hpp:436
GFitsImage
Abstract FITS image base class.
Definition GFitsImage.hpp:43
GFitsImage::pixel
virtual double pixel(const int &ix) const =0
GFitsTableCol
Abstract interface for FITS table column.
Definition GFitsTableCol.hpp:44
GFitsTableCol::integer
virtual int integer(const int &row, const int &inx=0) const =0
GFitsTableCol::real
virtual double real(const int &row, const int &inx=0) const =0
GFitsTableCol::unit
void unit(const std::string &unit)
Set column unit.
Definition GFitsTableCol.hpp:189
GFitsTableDoubleCol
FITS table double column.
Definition GFitsTableDoubleCol.hpp:42
GFitsTableShortCol
FITS table short integer column.
Definition GFitsTableShortCol.hpp:42
GFitsTable
Abstract interface for FITS table.
Definition GFitsTable.hpp:44
GFitsTable::append
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition GFitsTable.hpp:147
GFitsTable::nrows
const int & nrows(void) const
Return number of rows in table.
Definition GFitsTable.hpp:119
GFits
FITS file class.
Definition GFits.hpp:63
GFits::saveto
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition GFits.cpp:1293
GFits::contains
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition GFits.hpp:282
GFits::append
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition GFits.cpp:678
GFits::image
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition GFits.cpp:368
GFits::close
void close(void)
Close FITS file.
Definition GFits.cpp:1342
GFits::table
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition GFits.cpp:482
GIntegrals
Integration class for set of functions.
Definition GIntegrals.hpp:47
GIntegrals::romberg
GVector romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition GIntegrals.cpp:210
GIntegrals::fixed_iter
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition GIntegrals.hpp:179
GMatrix
Generic matrix class definition.
Definition GMatrix.hpp:79
GMatrix::eulerz
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition GMatrix.cpp:1194
GMatrix::eulery
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition GMatrix.cpp:1157
GModelSky
Sky model class.
Definition GModelSky.hpp:122
GModelSky::spatial
GModelSpatial * spatial(void) const
Return spatial model component.
Definition GModelSky.hpp:296
GModelSpatialDiffuseMap::map
const GSkyMap & map(void) const
Get map.
Definition GModelSpatialDiffuseMap.hpp:188
GModelSpatialDiffuse
Abstract diffuse spatial model base class.
Definition GModelSpatialDiffuse.hpp:52
GModelSpatialDiffuse::eval
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
GModelSpatialElliptical
Abstract elliptical spatial model base class.
Definition GModelSpatialElliptical.hpp:55
GModelSpatialElliptical::eval
virtual double eval(const double &theta, const double &posangle, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
GModelSpatialElliptical::dir
const GSkyDir & dir(void) const
Return position of elliptical spatial model.
Definition GModelSpatialElliptical.cpp:306
GModelSpatialElliptical::coordsys
std::string coordsys(void) const
Return coordinate system.
Definition GModelSpatialElliptical.hpp:258
GModelSpatialElliptical::theta_max
virtual double theta_max(void) const =0
GModelSpatialPointSource
Point source spatial model.
Definition GModelSpatialPointSource.hpp:57
GModelSpatialPointSource::dir
const GSkyDir & dir(void) const
Return position of point source.
Definition GModelSpatialPointSource.cpp:565
GModelSpatialRadial
Abstract radial spatial model base class.
Definition GModelSpatialRadial.hpp:55
GModelSpatialRadial::theta_max
virtual double theta_max(void) const =0
GModelSpatialRadial::eval
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
GModelSpatialRadial::dir
const GSkyDir & dir(void) const
Return position of radial spatial model.
Definition GModelSpatialRadial.cpp:297
GModel::name
const std::string & name(void) const
Return parameter name.
Definition GModel.hpp:265
GNodeArray::set_value
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition GNodeArray.cpp:587
GNodeArray::inx_right
const int & inx_right(void) const
Returns right node index.
Definition GNodeArray.hpp:249
GNodeArray::inx_left
const int & inx_left(void) const
Returns left node index.
Definition GNodeArray.hpp:235
GNodeArray::reserve
void reserve(const int &num)
Reserves space for nodes in node array.
Definition GNodeArray.hpp:220
GNodeArray::wgt_right
const double & wgt_right(void) const
Returns right node weight.
Definition GNodeArray.hpp:279
GNodeArray::wgt_left
const double & wgt_left(void) const
Returns left node weight.
Definition GNodeArray.hpp:264
GNodeArray::clear
void clear(void)
Clear node array.
Definition GNodeArray.cpp:235
GNodeArray::size
int size(void) const
Return number of nodes in node array.
Definition GNodeArray.hpp:192
GNodeArray::append
void append(const double &node)
Append one node to array.
Definition GNodeArray.cpp:351
GObservation
Abstract observation base class.
Definition GObservation.hpp:70
GObservation::events
virtual GEvents * events(void)
Return events.
Definition GObservation.cpp:176
GPhoton
Class that handles photons.
Definition GPhoton.hpp:47
GPhoton::dir
const GSkyDir & dir(void) const
Return photon sky direction.
Definition GPhoton.hpp:110
GPhoton::energy
const GEnergy & energy(void) const
Return photon energy.
Definition GPhoton.hpp:122
GResponse
Abstract instrument response base class.
Definition GResponse.hpp:77
GResponse::irf_diffuse
virtual double irf_diffuse(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to diffuse source.
Definition GResponse.cpp:1125
GResponse::init_members
void init_members(void)
Initialise class members.
Definition GResponse.cpp:794
GResponse::free_members
void free_members(void)
Delete class members.
Definition GResponse.cpp:842
GResponse::operator=
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition GResponse.cpp:139
GSPIEventBin
INTEGRAL/SPI event bin class.
Definition GSPIEventBin.hpp:62
GSPIEventBin::ipt
const int & ipt(void) const
Return event bin pointing index.
Definition GSPIEventBin.hpp:268
GSPIEventBin::iebin
const int & iebin(void) const
Return event bin energy index.
Definition GSPIEventBin.hpp:296
GSPIEventBin::dir
virtual const GSPIInstDir & dir(void) const
Return instrument direction.
Definition GSPIEventBin.hpp:142
GSPIEventBin::livetime
const double & livetime(void) const
Return livetime of event bin.
Definition GSPIEventBin.hpp:240
GSPIEventCube
INTEGRAL/SPI event bin container class.
Definition GSPIEventCube.hpp:56
GSPIEventCube::dir
const GSPIInstDir & dir(const int &ipt, const int &idet) const
Return instrument direction.
Definition GSPIEventCube.cpp:715
GSPIEventCube::spi_z
const GSkyDir & spi_z(const int &ipt) const
Return SPI Z direction.
Definition GSPIEventCube.cpp:772
GSPIEventCube::m_livetime
double * m_livetime
Livetime array.
Definition GSPIEventCube.hpp:132
GSPIEventCube::naxis
virtual int naxis(const int &axis) const
Return number of bins in axis.
Definition GSPIEventCube.cpp:263
GSPIEventCube::m_energy
GEnergy * m_energy
Energy array.
Definition GSPIEventCube.hpp:139
GSPIEventCube::m_num_pt
int m_num_pt
Number of pointings.
Definition GSPIEventCube.hpp:122
GSPIEventCube::m_detid
int * m_detid
Detector ID of event cube ADD!!!!
Definition GSPIEventCube.hpp:130
GSPIEventCube::m_num_ebin
int m_num_ebin
Number of energy bins.
Definition GSPIEventCube.hpp:124
GSPIEventCube::m_num_det
int m_num_det
Number of detectors.
Definition GSPIEventCube.hpp:123
GSPIEventCube::spi_x
const GSkyDir & spi_x(const int &ipt) const
Return SPI X direction (pointing direction)
Definition GSPIEventCube.cpp:746
GSPIInstDir::detid
void detid(const int &detid)
Set detector identifier.
Definition GSPIInstDir.hpp:138
GSPIObservation
INTEGRAL/SPI observation class.
Definition GSPIObservation.hpp:54
GSPIResponse
INTEGRAL/SPI instrument response function class.
Definition GSPIResponse.hpp:89
GSPIResponse::irf_init
void irf_init(const GObservation &obs, GSPIResponseIrf *irf) const
Initialise structure for projection computation of response.
Definition GSPIResponse.cpp:2477
GSPIResponse::m_wcs_xpix_max
double m_wcs_xpix_max
Maximum X pixel index.
Definition GSPIResponse.hpp:216
GSPIResponse::set_detids
void set_detids(const GSPIEventCube *cube)
Set vector of detector identifiers used by the observation.
Definition GSPIResponse.cpp:1426
GSPIResponse::operator=
virtual GSPIResponse & operator=(const GSPIResponse &rsp)
Assignment operator.
Definition GSPIResponse.cpp:176
GSPIResponse::write
void write(GFits &fits) const
Write SPI response into FITS object.
Definition GSPIResponse.cpp:605
GSPIResponse::irf_map
int irf_map(const int &idet, const int &ieng, const GSPIResponseIrf *irf) const
Return map index for detector and energy index.
Definition GSPIResponse.hpp:400
GSPIResponse::irf_skydir
void irf_skydir(const int &ipt, const int &inx, GSPIResponseIrf *irf) const
Set sky direction for pointing and IRF pixel index.
Definition GSPIResponse.cpp:2607
GSPIResponse::irf_radial
virtual GVector irf_radial(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to radial source sky model.
Definition GSPIResponse.cpp:1757
GSPIResponse::irf_ptsrc
virtual GVector irf_ptsrc(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to point source sky model.
Definition GSPIResponse.cpp:1678
GSPIResponse::irf_elliptical
virtual GVector irf_elliptical(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to elliptical source sky model.
Definition GSPIResponse.cpp:1979
GSPIResponse::m_wcs_xbin
double m_wcs_xbin
X value bin size (radians)
Definition GSPIResponse.hpp:214
GSPIResponse::save
void save(const GFilename &filename, const bool &clobber=false) const
Save SPI response into file.
Definition GSPIResponse.cpp:652
GSPIResponse::m_irfs
GSkyMap m_irfs
IRFs stored as sky map.
Definition GSPIResponse.hpp:201
GSPIResponse::m_energies
GNodeArray m_energies
Node array of IRF energies.
Definition GSPIResponse.hpp:199
GSPIResponse::nroi
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.
Definition GSPIResponse.cpp:431
GSPIResponse::irf_weight
double irf_weight(const double &beta, const double &emin, const double &emax) const
Compute weight of logarithmic energy bin.
Definition GSPIResponse.cpp:1641
GSPIResponse::m_wcs_ymin
double m_wcs_ymin
Minimum Y value (radians)
Definition GSPIResponse.hpp:211
GSPIResponse::print
virtual std::string print(const GChatter &chatter=NORMAL) const
Print INTEGRAL/SPI response information.
Definition GSPIResponse.cpp:677
GSPIResponse::m_wcs_ybin
double m_wcs_ybin
Y value bin size (radians)
Definition GSPIResponse.hpp:215
GSPIResponse::set_wcs
void set_wcs(const GFitsImage *image)
Set IRF image limits.
Definition GSPIResponse.cpp:1380
GSPIResponse::irf_value
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.
Definition GSPIResponse.cpp:340
GSPIResponse::read
void read(const GFits &fits)
Read SPI response from FITS object.
Definition GSPIResponse.cpp:576
GSPIResponse::compute_irf
GSkyMap compute_irf(const double &emin, const double &emax) const
Compute as sky map.
Definition GSPIResponse.cpp:1234
GSPIResponse::irf
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return value of INTEGRAL/SPI instrument response for a photon.
Definition GSPIResponse.cpp:260
GSPIResponse::load_irf
GSkyMap load_irf(const GFilename &rspname) const
Load IRF as sky map.
Definition GSPIResponse.cpp:1021
GSPIResponse::m_spix
std::vector< GSkyDir > m_spix
SPI pointing direction.
Definition GSPIResponse.hpp:207
GSPIResponse::m_dlogE
double m_dlogE
Logarithmic energy step for IRF band.
Definition GSPIResponse.hpp:203
GSPIResponse::write_detids
void write_detids(GFits &fits) const
Write detector identifiers into FITS object.
Definition GSPIResponse.cpp:929
GSPIResponse::init_members
void init_members(void)
Initialise class members.
Definition GSPIResponse.cpp:743
GSPIResponse::~GSPIResponse
virtual ~GSPIResponse(void)
Destructor.
Definition GSPIResponse.cpp:149
GSPIResponse::azimuth
double azimuth(const int &ipt, const GSkyDir &dir) const
Return azimuth angle of sky direction for pointing in radians.
Definition GSPIResponse.hpp:374
GSPIResponse::copy_members
void copy_members(const GSPIResponse &rsp)
Copy class members.
Definition GSPIResponse.cpp:779
GSPIResponse::m_detids
std::vector< int > m_detids
Vector of detector IDs.
Definition GSPIResponse.hpp:198
GSPIResponse::irf_diffuse
virtual GVector irf_diffuse(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to diffuse source sky model.
Definition GSPIResponse.cpp:2203
GSPIResponse::free_members
void free_members(void)
Delete class members.
Definition GSPIResponse.cpp:813
GSPIResponse::m_rspname
GFilename m_rspname
File name of response group.
Definition GSPIResponse.hpp:197
GSPIResponse::write_energies
void write_energies(GFits &fits) const
Write energies into FITS object.
Definition GSPIResponse.cpp:963
GSPIResponse::load
void load(const GFilename &filename)
Load SPI response from file.
Definition GSPIResponse.cpp:628
GSPIResponse::m_wcs_ypix_max
double m_wcs_ypix_max
Maximum Y pixel index.
Definition GSPIResponse.hpp:217
GSPIResponse::m_gamma
double m_gamma
Power-law spectral index for IRF band.
Definition GSPIResponse.hpp:204
GSPIResponse::read_detids
void read_detids(const GFits &fits)
Read detector identifiers from FITS object.
Definition GSPIResponse.cpp:827
GSPIResponse::set_cache
void set_cache(const GSPIEventCube *cube)
Set computation cache.
Definition GSPIResponse.cpp:1466
GSPIResponse::irf_detid
int irf_detid(const int &detid) const
Convert detector identifier into IRF detector identifier.
Definition GSPIResponse.cpp:1613
GSPIResponse::zenith
double zenith(const int &ipt, const GSkyDir &dir) const
Return zenith angle of sky direction for pointing in radians.
Definition GSPIResponse.hpp:358
GSPIResponse::m_wcs_xmax
double m_wcs_xmax
Maximum X value (radians)
Definition GSPIResponse.hpp:212
GSPIResponse::m_wcs_ymax
double m_wcs_ymax
Maximum Y value (radians)
Definition GSPIResponse.hpp:213
GSPIResponse::irf_rot_matrix
void irf_rot_matrix(const int &ipt, GSPIResponseIrf *irf) const
Set rotation matrix for pointing.
Definition GSPIResponse.cpp:2569
GSPIResponse::irf_vector
void irf_vector(const GSPIEventCube *cube, const double &xpix, const double &ypix, const int &ireg, const double *livetimes, GVector *irf) const
Fill vector of INTEGRAL/SPI instrument response for IRF pixel.
Definition GSPIResponse.cpp:1521
GSPIResponse::clone
virtual GSPIResponse * clone(void) const
Clone instance.
Definition GSPIResponse.cpp:232
GSPIResponse::m_wcs_xmin
double m_wcs_xmin
Minimum X value (radians)
Definition GSPIResponse.hpp:210
GSPIResponse::m_energy_keV
double m_energy_keV
IRF line energy (optional)
Definition GSPIResponse.hpp:202
GSPIResponse::m_ebounds
GEbounds m_ebounds
Energy bounaries of IRF.
Definition GSPIResponse.hpp:200
GSPIResponse::m_max_zenith
double m_max_zenith
Maximum zenith angle (radians)
Definition GSPIResponse.hpp:218
GSPIResponse::read_energies
void read_energies(const GFits &fits)
Read energies from FITS object.
Definition GSPIResponse.cpp:860
GSPIResponse::GSPIResponse
GSPIResponse(void)
Void constructor.
Definition GSPIResponse.cpp:96
GSPIResponse::clear
virtual void clear(void)
Clear instance.
Definition GSPIResponse.cpp:212
GSPIResponse::rspname
const GFilename & rspname(void) const
Get response group file name.
Definition GSPIResponse.hpp:285
GSPIResponse::m_has_wcs
bool m_has_wcs
Has WCS information.
Definition GSPIResponse.hpp:209
GSPIResponse::ebounds
virtual GEbounds ebounds(const GEnergy &obsEnergy) const
Return true energy boundaries for a specific observed energy.
Definition GSPIResponse.cpp:457
GSPIResponse::set
void set(const GSPIObservation &obs, const GEnergy &energy=GEnergy())
Set response for a specific observation.
Definition GSPIResponse.cpp:485
GSPIResponse::m_posang
std::vector< double > m_posang
Position angle of Y axis (CEL, radians)
Definition GSPIResponse.hpp:208
GSPIResponse::load_irfs
void load_irfs(const int &region)
Load Instrument Response Functions.
Definition GSPIResponse.cpp:1128
GSPIResponse::irf_diffuse_map
void irf_diffuse_map(const GModelSpatialDiffuseMap *diffuse, const GObservation &obs, GVector *irfs, GMatrix *gradients=NULL) const
Return instrument response to diffuse map sky model.
Definition GSPIResponse.cpp:2366
GSkyDir
Sky direction class.
Definition GSkyDir.hpp:62
GSkyDir::dec_deg
double dec_deg(void) const
Returns Declination in degrees.
Definition GSkyDir.hpp:258
GSkyDir::ra_deg
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition GSkyDir.hpp:231
GSkyDir::dist
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition GSkyDir.hpp:273
GSkyDir::posang
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition GSkyDir.cpp:1300
GSkyMap
Sky map class.
Definition GSkyMap.hpp:89
GSkyMap::solidangle
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition GSkyMap.cpp:1858
GSkyMap::nmaps
const int & nmaps(void) const
Returns number of maps.
Definition GSkyMap.hpp:427
GSkyMap::clear
void clear(void)
Clear instance.
Definition GSkyMap.cpp:1100
GSkyMap::inx2dir
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition GSkyMap.cpp:1331
GSkyMap::nx
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition GSkyMap.hpp:399
GSkyMap::read
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition GSkyMap.cpp:2664
GSkyMap::npix
const int & npix(void) const
Returns number of pixels.
Definition GSkyMap.hpp:383
GSkyMap::ny
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition GSkyMap.hpp:415
GSkyMap::write
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition GSkyMap.cpp:2697
GSkyMap::shape
const std::vector< int > & shape(void) const
Returns shape of maps.
Definition GSkyMap.hpp:439
GSource
Class that handles gamma-ray sources.
Definition GSource.hpp:53
GTime
Time class.
Definition GTime.hpp:55
GVector
Vector class.
Definition GVector.hpp:46
spi_elliptical_kerns_rho
Kernel for rho angle integration of elliptical models.
Definition spi_helpers_response_vector.hpp:129
spi_radial_kerns_rho
Kernel for rho angle integration of radial models.
Definition spi_helpers_response_vector.hpp:49
gammalib::parformat
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition GTools.cpp:1162
gammalib::is_infinite
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition GTools.hpp:186
gammalib::is_notanumber
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition GTools.hpp:203
gammalib::ln10
const double ln10
Definition GMath.hpp:46
gammalib::extname_ebounds
const std::string extname_ebounds
Definition GEbounds.hpp:44
gammalib::str
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition GTools.cpp:508
gammalib::pihalf
const double pihalf
Definition GMath.hpp:38
gammalib::spi_num_hdus
int spi_num_hdus(const GFits &fits, const std::string &extname)
Return number of HDU versions.
Definition GSPITools.cpp:166
gammalib::deg2rad
const double deg2rad
Definition GMath.hpp:43
gammalib::get_current_clock
double get_current_clock(void)
Get current clock in seconds.
Definition GTools.cpp:2587
gammalib::spi_hdu
const GFitsTable * spi_hdu(const GFits &fits, const std::string &extname, const int &extver=1)
Return FITS table.
Definition GSPITools.cpp:57
gammalib::rad2deg
const double rad2deg
Definition GMath.hpp:44
gammalib::extname_energies
const std::string extname_energies
Definition GEnergies.hpp:44
spi_helpers_response_vector.hpp
Defintion of SPI helper classes for vector response.