Rolf Buehler / gammalib

Commit b6f589051644c560941df406910b51ab7d9e76b3

Authored by Jürgen Knödlseder
1 parent dd822618

Add initial versions of GCOMD1Response and GCOMD2Response classes (#2186)

Inline Side-by-side
Showing 13 changed files with 1471 additions and 14 deletions
... ... @@ -12,7 +12,8 @@ made since the last release of GammaLib.
12 12 2. Public interface modifications
13 13 ---------------------------------
14 14 The following classes have been added:
15   -- n.a.
  15 +- GCOMD1Response
  16 +- GCOMD2Response
16 17  
17 18 The following classes have been removed:
18 19 - n.a.
... ... @@ -49,8 +50,8 @@ None
49 50  
50 51 4. Application module
51 52 ---------------------
52   -A time parameter with type "t" has been added to allow specifying time values in
53   -various formats to time systems to an application. Methods have been added to set
  53 +A time parameter with type "t" has was to allow specifying time values in
  54 +various formats to time systems to an application. Methods were added to set
54 55 and retrieve the time value (#1864).
55 56  
56 57  
... ... @@ -82,11 +83,11 @@ None
82 83  
83 84 10. Observation module
84 85 ----------------------
85   -A time string constructor and set method has been added to GTime to allow
  86 +A time string constructor and set method was added to GTime to allow
86 87 setting various time formats from a string. Variants exist to set the time in
87 88 the native time reference or in a specific time reference (#1864).
88 89  
89   -A rounding error has been fixed in GTime::utc() so that times such as
  90 +A rounding error was fixed in GTime::utc() so that times such as
90 91 "2019-12-31T23:59:60" are now "2020-01-01T00:00:00" as they should be (#2177).
91 92  
92 93  
... ... @@ -145,3 +146,7 @@ None
145 146 COMPTEL SDA, SDB and ICT response files were added so that the information to
146 147 compute the instrument response for any arbitrary energy band is available
147 148 (#2186).
  149 +
  150 +The classes GCOMD1Response and GCOMD2Response were added to handle the energy
  151 +response of the D1 and D2 modules, respectively. These classes handle reading
  152 +of the SDA and SDB response files (#2186).
... ...
include/GMath.hpp
  View file @b6f5890
... ... @@ -49,8 +49,9 @@ namespace gammalib {
49 49 const double twothird = 2.0/3.0;
50 50 const double fourthird = 4.0/3.0;
51 51 const double sqrt_onehalf = std::sqrt(1.0/2.0);
52   - const double sqrt_pihalf = std::sqrt(pihalf);
53 52 const double sqrt_two = std::sqrt(2.0);
  53 + const double sqrt_pihalf = std::sqrt(pihalf);
  54 + const double sqrt_twopi = std::sqrt(twopi);
54 55 const double inv_ln10 = 1.0/ln10;
55 56 }
56 57  
... ...
inst/com/Makefile.am
  View file @b6f5890
1 1 #############################################################################
2 2 # Makefile.am for com module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2017 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -39,6 +39,8 @@ sources = src/GCOMSupport.cpp \
39 39 src/GCOMEventBin.cpp \
40 40 src/GCOMInstDir.cpp \
41 41 src/GCOMResponse.cpp \
  42 + src/GCOMD1Response.cpp \
  43 + src/GCOMD2Response.cpp \
42 44 src/GCOMObservation.cpp \
43 45 src/GCOMModelDRBFitting.cpp
44 46  
... ... @@ -47,6 +49,8 @@ pkginclude_HEADERS = include/GCOMEventCube.hpp \
47 49 include/GCOMEventBin.hpp \
48 50 include/GCOMInstDir.hpp \
49 51 include/GCOMResponse.hpp \
  52 + include/GCOMD1Response.hpp \
  53 + include/GCOMD2Response.hpp \
50 54 include/GCOMObservation.hpp \
51 55 include/GCOMModelDRBFitting.hpp \
52 56 include/GCOMLib.hpp
... ...
inst/com/include/GCOMD1Response.hpp 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD1Response.hpp - COMPTEL D1 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD1Response.hpp
  23 + * @brief COMPTEL D1 module response class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCOMD1RESPONSE_HPP
  28 +#define GCOMD1RESPONSE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GCaldb.hpp"
  34 +#include "GEnergy.hpp"
  35 +#include "GNodeArray.hpp"
  36 +
  37 +/* __ Type definitions ___________________________________________________ */
  38 +
  39 +/* __ Forward declarations _______________________________________________ */
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCOMD1Response
  44 + *
  45 + * @brief Interface for the COMPTEL D1 module response class
  46 + ***************************************************************************/
  47 +class GCOMD1Response : public GBase {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GCOMD1Response(void);
  52 + GCOMD1Response(const GCOMD1Response& rsp);
  53 + GCOMD1Response(const std::string& sdaname, const GCaldb& caldb);
  54 + ~GCOMD1Response(void);
  55 +
  56 + // Operators
  57 + GCOMD1Response& operator=(const GCOMD1Response & rsp);
  58 + double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCOMD1Response* clone(void) const;
  63 + std::string classname(void) const;
  64 + void caldb(const GCaldb& caldb);
  65 + const GCaldb& caldb(void) const;
  66 + void load(const std::string& sdaname);
  67 + void read(const GFitsTable& table);
  68 + std::string print(const GChatter& chatter = NORMAL) const;
  69 +
  70 +private:
  71 + // Private methods
  72 + void init_members(void);
  73 + void copy_members(const GCOMD1Response& rsp);
  74 + void free_members(void);
  75 + void update_cache(const GEnergy& etrue) const;
  76 +
  77 + // Private data members
  78 + GCaldb m_caldb; //!< Calibration database
  79 + GNodeArray m_energies; //!< Input energies
  80 + std::vector<double> m_positions; //!< Photo peak position in MeV
  81 + std::vector<double> m_sigmas; //!< Photo peak width in MeV
  82 + std::vector<double> m_amplitudes; //!< Photo peak amplitude
  83 + std::vector<double> m_emins; //!< Lower energy threshold of D1
  84 + std::vector<double> m_ewidths; //!< Lower energy threshold width of D1
  85 + std::vector<double> m_emaxs; //!< Upper energy limit of D1
  86 +
  87 + // Pre-computation cache
  88 + mutable GEnergy m_energy;
  89 + mutable double m_position;
  90 + mutable double m_sigma;
  91 + mutable double m_amplitude;
  92 + mutable double m_emin;
  93 + mutable double m_ewidth;
  94 + mutable double m_emax;
  95 +};
  96 +
  97 +
  98 +/***********************************************************************//**
  99 + * @brief Return class name
  100 + *
  101 + * @return String containing the class name ("GCOMD1Response").
  102 + ***************************************************************************/
  103 +inline
  104 +std::string GCOMD1Response::classname(void) const
  105 +{
  106 + return ("GCOMD1Response");
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Return calibration database
  112 + *
  113 + * @return Calibration database.
  114 + ***************************************************************************/
  115 +inline
  116 +const GCaldb& GCOMD1Response::caldb(void) const
  117 +{
  118 + return m_caldb;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Set calibration database
  124 + *
  125 + * @param[in] caldb Calibration database.
  126 + *
  127 + * Sets the calibration database containing the COMPTEL D1 module response.
  128 + ***************************************************************************/
  129 +inline
  130 +void GCOMD1Response::caldb(const GCaldb& caldb)
  131 +{
  132 + m_caldb = caldb;
  133 + return;
  134 +}
  135 +
  136 +#endif /* GCOMD1RESPONSE_HPP */
... ...
inst/com/include/GCOMD2Response.hpp 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD2Response.hpp - COMPTEL D2 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD2Response.hpp
  23 + * @brief COMPTEL D2 module response class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCOMD2RESPONSE_HPP
  28 +#define GCOMD2RESPONSE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GCaldb.hpp"
  34 +#include "GEnergy.hpp"
  35 +#include "GNodeArray.hpp"
  36 +
  37 +/* __ Type definitions ___________________________________________________ */
  38 +
  39 +/* __ Forward declarations _______________________________________________ */
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCOMD2Response
  44 + *
  45 + * @brief Interface for the COMPTEL D2 module response class
  46 + ***************************************************************************/
  47 +class GCOMD2Response : public GBase {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GCOMD2Response(void);
  52 + GCOMD2Response(const GCOMD2Response& rsp);
  53 + GCOMD2Response(const std::string& sdbname, const GCaldb& caldb);
  54 + ~GCOMD2Response(void);
  55 +
  56 + // Operators
  57 + GCOMD2Response& operator=(const GCOMD2Response & rsp);
  58 + double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCOMD2Response* clone(void) const;
  63 + std::string classname(void) const;
  64 + void caldb(const GCaldb& caldb);
  65 + const GCaldb& caldb(void) const;
  66 + void load(const std::string& sdbname);
  67 + void read(const GFitsTable& table);
  68 + std::string print(const GChatter& chatter = NORMAL) const;
  69 +
  70 +private:
  71 + // Private methods
  72 + void init_members(void);
  73 + void copy_members(const GCOMD2Response& rsp);
  74 + void free_members(void);
  75 + void update_cache(const GEnergy& etrue) const;
  76 +
  77 + // Private data members
  78 + GCaldb m_caldb; //!< Calibration database
  79 + GNodeArray m_energies; //!< Input energies
  80 + std::vector<double> m_positions; //!< Photo peak position in MeV
  81 + std::vector<double> m_sigmas; //!< Photo peak width in MeV
  82 + std::vector<double> m_amplitudes; //!< Photo peak amplitude
  83 + std::vector<double> m_escapes1; //!< Amplitude of first escape peak
  84 + std::vector<double> m_escapes2; //!< Amplitude of second escape peak
  85 + std::vector<double> m_tails; //!< Amplitude of Compton tail
  86 + std::vector<double> m_backgrounds; //!< Amplitude of Compton background
  87 + std::vector<double> m_emins; //!< Lower energy threshold of D2
  88 + std::vector<double> m_ewidths; //!< Lower energy threshold width of D2
  89 + std::vector<double> m_emaxs; //!< Upper energy limit of D2
  90 +
  91 + // Pre-computation cache
  92 + mutable GEnergy m_energy;
  93 + mutable double m_position;
  94 + mutable double m_sigma;
  95 + mutable double m_amplitude;
  96 + mutable double m_escape1;
  97 + mutable double m_escape2;
  98 + mutable double m_tail;
  99 + mutable double m_background;
  100 + mutable double m_emin;
  101 + mutable double m_ewidth;
  102 + mutable double m_emax;
  103 +};
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Return class name
  108 + *
  109 + * @return String containing the class name ("GCOMD2Response").
  110 + ***************************************************************************/
  111 +inline
  112 +std::string GCOMD2Response::classname(void) const
  113 +{
  114 + return ("GCOMD2Response");
  115 +}
  116 +
  117 +
  118 +/***********************************************************************//**
  119 + * @brief Return calibration database
  120 + *
  121 + * @return Calibration database.
  122 + ***************************************************************************/
  123 +inline
  124 +const GCaldb& GCOMD2Response::caldb(void) const
  125 +{
  126 + return m_caldb;
  127 +}
  128 +
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Set calibration database
  132 + *
  133 + * @param[in] caldb Calibration database.
  134 + *
  135 + * Sets the calibration database containing the COMPTEL D2 module response.
  136 + ***************************************************************************/
  137 +inline
  138 +void GCOMD2Response::caldb(const GCaldb& caldb)
  139 +{
  140 + m_caldb = caldb;
  141 + return;
  142 +}
  143 +
  144 +#endif /* GCOMD2RESPONSE_HPP */
... ...
inst/com/include/GCOMLib.hpp
  View file @b6f5890
1 1 /***************************************************************************
2 2 * GCOMLib.hpp - COMPTEL Support Header files *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,6 +35,8 @@
35 35 #include "GCOMEventBin.hpp"
36 36 #include "GCOMInstDir.hpp"
37 37 #include "GCOMResponse.hpp"
  38 +#include "GCOMD1Response.hpp"
  39 +#include "GCOMD2Response.hpp"
38 40 #include "GCOMObservation.hpp"
39 41 #include "GCOMModelDRBFitting.hpp"
40 42  
... ...
inst/com/include/GCOMResponse.hpp
  View file @b6f5890
... ... @@ -94,7 +94,7 @@ private:
94 94 void free_members(void);
95 95  
96 96 // Private data members
97   - GCaldb m_caldb; //!< Name of or path to the calibration database
  97 + GCaldb m_caldb; //!< Calibration database
98 98 std::string m_rspname; //!< Response name
99 99 std::vector<double> m_iaq; //!< IAQ array
100 100 int m_phigeo_bins; //!< Number of Phigeo bins
... ... @@ -159,7 +159,7 @@ const GCaldb&amp; GCOMResponse::caldb(void) const
159 159  
160 160  
161 161 /***********************************************************************//**
162   - * @brief Set path to the calibration database
  162 + * @brief Set calibration database
163 163 *
164 164 * @param[in] caldb Calibration database.
165 165 *
... ...
inst/com/pyext/GCOMD1Response.i 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD1Response.i - COMPTEL D1 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD1Response.i
  23 + * @brief COMPTEL D1 module response class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCOMD1Response.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCOMD1Response
  34 + *
  35 + * @brief Interface for the COMPTEL D1 module response class
  36 + ***************************************************************************/
  37 +class GCOMD1Response : public GBase {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GCOMD1Response(void);
  42 + GCOMD1Response(const GCOMD1Response& rsp);
  43 + GCOMD1Response(const std::string& sdaname, const GCaldb& caldb);
  44 + ~GCOMD1Response(void);
  45 +
  46 + // Operators
  47 + double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  48 +
  49 + // Methods
  50 + void clear(void);
  51 + GCOMD1Response* clone(void) const;
  52 + std::string classname(void) const;
  53 + void caldb(const GCaldb& caldb);
  54 + const GCaldb& caldb(void) const;
  55 + void load(const std::string& sdaname);
  56 + void read(const GFitsTable& hdu);
  57 +};
  58 +
  59 +
  60 +/***********************************************************************//**
  61 + * @brief GCOMD1Response class extension
  62 + ***************************************************************************/
  63 +%extend GCOMD1Response {
  64 + GCOMD1Response copy() {
  65 + return (*self);
  66 + }
  67 +};
... ...
inst/com/pyext/GCOMD2Response.i 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD2Response.i - COMPTEL D2 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD2Response.i
  23 + * @brief COMPTEL D2 module response class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCOMD2Response.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCOMD2Response
  34 + *
  35 + * @brief Interface for the COMPTEL D2 module response class
  36 + ***************************************************************************/
  37 +class GCOMD2Response : public GBase {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GCOMD2Response(void);
  42 + GCOMD2Response(const GCOMD2Response& rsp);
  43 + GCOMD2Response(const std::string& sdaname, const GCaldb& caldb);
  44 + ~GCOMD2Response(void);
  45 +
  46 + // Operators
  47 + double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  48 +
  49 + // Methods
  50 + void clear(void);
  51 + GCOMD2Response* clone(void) const;
  52 + std::string classname(void) const;
  53 + void caldb(const GCaldb& caldb);
  54 + const GCaldb& caldb(void) const;
  55 + void load(const std::string& sdaname);
  56 + void read(const GFitsTable& hdu);
  57 +};
  58 +
  59 +
  60 +/***********************************************************************//**
  61 + * @brief GCOMD2Response class extension
  62 + ***************************************************************************/
  63 +%extend GCOMD2Response {
  64 + GCOMD2Response copy() {
  65 + return (*self);
  66 + }
  67 +};
... ...
inst/com/pyext/com.i
  View file @b6f5890
1 1 /***************************************************************************
2 2 * com.i - COMPTEL module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -67,4 +67,5 @@
67 67 %include "GCOMModelDRBFitting.i"
68 68 %include "GCOMObservation.i"
69 69 %include "GCOMResponse.i"
70   -
  70 +%include "GCOMD1Response.i"
  71 +%include "GCOMD2Response.i"
... ...
inst/com/src/GCOMD1Response.cpp 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD1Response.cpp - COMPTEL D1 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD1Response.cpp
  23 + * @brief COMPTEL D1 module response class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GMath.hpp"
  32 +#include "GCOMD1Response.hpp"
  33 +
  34 +/* __ Method name definitions ____________________________________________ */
  35 +
  36 +/* __ Macros _____________________________________________________________ */
  37 +
  38 +/* __ Coding definitions _________________________________________________ */
  39 +
  40 +/* __ Debug definitions __________________________________________________ */
  41 +
  42 +/* __ Constants __________________________________________________________ */
  43 +
  44 +
  45 +/*==========================================================================
  46 + = =
  47 + = Constructors/destructors =
  48 + = =
  49 + ==========================================================================*/
  50 +
  51 +/***********************************************************************//**
  52 + * @brief Void constructor
  53 + *
  54 + * Creates an empty COMPTEL D1 module response.
  55 + ***************************************************************************/
  56 +GCOMD1Response::GCOMD1Response(void)
  57 +{
  58 + // Initialise members
  59 + init_members();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Copy constructor
  68 + *
  69 + * @param[in] rsp COMPTEL D1 module response.
  70 + **************************************************************************/
  71 +GCOMD1Response::GCOMD1Response(const GCOMD1Response& rsp)
  72 +{
  73 + // Initialise members
  74 + init_members();
  75 +
  76 + // Copy members
  77 + copy_members(rsp);
  78 +
  79 + // Return
  80 + return;
  81 +}
  82 +
  83 +
  84 +/***********************************************************************//**
  85 + * @brief Response constructor
  86 + *
  87 + * @param[in] sdaname SDA response name.
  88 + * @param[in] caldb Calibration database.
  89 + *
  90 + * Create COMPTEL D1 module response response by loading an SDA file from a
  91 + * calibration database.
  92 + ***************************************************************************/
  93 +GCOMD1Response::GCOMD1Response(const std::string& sdaname, const GCaldb& caldb)
  94 +{
  95 + // Initialise members
  96 + init_members();
  97 +
  98 + // Set calibration database
  99 + this->caldb(caldb);
  100 +
  101 + // Load D1 module response
  102 + this->load(sdaname);
  103 +
  104 + // Return
  105 + return;
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Destructor
  111 + *
  112 + * Destroys instance of COMPTEL response object.
  113 + ***************************************************************************/
  114 +GCOMD1Response::~GCOMD1Response(void)
  115 +{
  116 + // Free members
  117 + free_members();
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/*==========================================================================
  125 + = =
  126 + = Operators =
  127 + = =
  128 + ==========================================================================*/
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Assignment operator
  132 + *
  133 + * @param[in] rsp COMPTEL D1 module response.
  134 + * @return COMPTEL D1 module response.
  135 + ***************************************************************************/
  136 +GCOMD1Response& GCOMD1Response::operator=(const GCOMD1Response& rsp)
  137 +{
  138 + // Execute only if object is not identical
  139 + if (this != &rsp) {
  140 +
  141 + // Free members
  142 + free_members();
  143 +
  144 + // Initialise members
  145 + init_members();
  146 +
  147 + // Copy members
  148 + copy_members(rsp);
  149 +
  150 + } // endif: object was not identical
  151 +
  152 + // Return this object
  153 + return *this;
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief D1 module response evaluation operator
  159 + *
  160 + * @param[in] etrue True energy.
  161 + * @param[in] ereco Reconstructed energy.
  162 + * @return COMPTEL D1 module response.
  163 + ***************************************************************************/
  164 +double GCOMD1Response::operator()(const GEnergy& etrue, const GEnergy& ereco) const
  165 +{
  166 + // Initialise response with zero
  167 + double response = 0.0;
  168 +
  169 + // Continue only if a response was loaded
  170 + if (!m_energies.is_empty()) {
  171 +
  172 + // Update response evaluation cache
  173 + update_cache(etrue);
  174 +
  175 + // Continue only if amplitude is positive
  176 + if (m_amplitude > 0.0) {
  177 +
  178 + // Compute D1 module response (here is where the real magic
  179 + // happens)
  180 + double arg = (m_energy.MeV()-ereco.MeV()) / m_sigma;
  181 + response = m_amplitude * std::exp(-0.5 * arg * arg);
  182 +
  183 + }
  184 +
  185 + } // endif: response was loaded
  186 +
  187 + // Return response
  188 + return response;
  189 +}
  190 +
  191 +
  192 +/*==========================================================================
  193 + = =
  194 + = Public methods =
  195 + = =
  196 + ==========================================================================*/
  197 +
  198 +/***********************************************************************//**
  199 + * @brief Clear instance
  200 + *
  201 + * Clears COMPTEL D1 module response object by resetting all members to an
  202 + * initial state. Any information that was present in the object before will
  203 + * be lost.
  204 + ***************************************************************************/
  205 +void GCOMD1Response::clear(void)
  206 +{
  207 + // Free class members
  208 + free_members();
  209 +
  210 + // Initialise members
  211 + init_members();
  212 +
  213 + // Return
  214 + return;
  215 +}
  216 +
  217 +
  218 +/***********************************************************************//**
  219 + * @brief Clone instance
  220 + *
  221 + * @return Pointer to deep copy of COMPTEL D1 module response.
  222 + ***************************************************************************/
  223 +GCOMD1Response* GCOMD1Response::clone(void) const
  224 +{
  225 + return new GCOMD1Response(*this);
  226 +}
  227 +
  228 +
  229 +/***********************************************************************//**
  230 + * @brief Load COMPTEL D1 module response.
  231 + *
  232 + * @param[in] sdaname COMPTEL D1 module response name.
  233 + *
  234 + * Loads the COMPTEL D1 module response with specified name @p sdaname. The
  235 + * method first searchs for an appropriate response in the calibration
  236 + * database. If no appropriate response is found, the method takes the
  237 + * database root path and response name to build the full path to the
  238 + * response file, and tries to load the response from these paths.
  239 + ***************************************************************************/
  240 +void GCOMD1Response::load(const std::string& sdaname)
  241 +{
  242 + // Clear instance but conserve calibration database
  243 + GCaldb caldb = m_caldb;
  244 + clear();
  245 + m_caldb = caldb;
  246 +
  247 + // First attempt reading the response using the GCaldb interface
  248 + GFilename filename = m_caldb.filename("","","SDA","","",sdaname);
  249 +
  250 + // If filename is empty then build filename from CALDB root path and
  251 + // response name
  252 + if (filename.is_empty()) {
  253 + filename = gammalib::filepath(m_caldb.rootdir(), sdaname);
  254 + if (!filename.exists()) {
  255 + GFilename testname = filename + ".fits";
  256 + if (testname.exists()) {
  257 + filename = testname;
  258 + }
  259 + }
  260 + }
  261 +
  262 + // Open FITS file
  263 + GFits fits(filename);
  264 +
  265 + // Get SDA table
  266 + const GFitsTable& sda = *fits.table(1);
  267 +
  268 + // Read SDA
  269 + read(sda);
  270 +
  271 + // Close SDA FITS file
  272 + fits.close();
  273 +
  274 + // Return
  275 + return;
  276 +}
  277 +
  278 +
  279 +/***********************************************************************//**
  280 + * @brief Read COMPTEL D1 module response.
  281 + *
  282 + * @param[in] table FITS table.
  283 + *
  284 + * Read the COMPTEL D1 module response from a SDA FITS table.
  285 + ***************************************************************************/
  286 +void GCOMD1Response::read(const GFitsTable& table)
  287 +{
  288 + // Initialise COMPTEL D1 module response vectors
  289 + m_energies.clear();
  290 + m_positions.clear();
  291 + m_sigmas.clear();
  292 + m_amplitudes.clear();
  293 + m_emins.clear();
  294 + m_ewidths.clear();
  295 + m_emaxs.clear();
  296 +
  297 + // Extract number of entries in table
  298 + int num = table.nrows();
  299 +
  300 + // If there are entries then read them
  301 + if (num > 0) {
  302 +
  303 + // Get column pointers
  304 + const GFitsTableCol* ptr_energy = table["ENERGY"];
  305 + const GFitsTableCol* ptr_position = table["POSITION"];
  306 + const GFitsTableCol* ptr_sigma = table["WIDTH"];
  307 + const GFitsTableCol* ptr_amplitude = table["AMPLITUDE"];
  308 + const GFitsTableCol* ptr_emin = table["EMIN"];
  309 + const GFitsTableCol* ptr_ewidth = table["EWIDTH"];
  310 + const GFitsTableCol* ptr_emax = table["EMAX"];
  311 +
  312 + // Copy data from table into vectors
  313 + for (int i = 0; i < num; ++i) {
  314 + m_energies.append(ptr_energy->real(i));
  315 + m_positions.push_back(ptr_position->real(i));
  316 + m_sigmas.push_back(ptr_sigma->real(i));
  317 + m_amplitudes.push_back(ptr_amplitude->real(i));
  318 + m_emins.push_back(ptr_emin->real(i));
  319 + m_ewidths.push_back(ptr_ewidth->real(i));
  320 + m_emaxs.push_back(ptr_emax->real(i));
  321 + }
  322 +
  323 + } // endif: there were entries
  324 +
  325 + // Return
  326 + return;
  327 +}
  328 +
  329 +
  330 +/***********************************************************************//**
  331 + * @brief Print COMPTEL D1 module response information
  332 + *
  333 + * @param[in] chatter Chattiness.
  334 + * @return String containing COMPTEL D1 module response information.
  335 + ***************************************************************************/
  336 +std::string GCOMD1Response::print(const GChatter& chatter) const
  337 +{
  338 + // Initialise result string
  339 + std::string result;
  340 +
  341 + // Continue only if chatter is not silent
  342 + if (chatter != SILENT) {
  343 +
  344 + // Append header
  345 + result.append("=== GCOMD1Response ===");
  346 +
  347 + // Append D1 module response information
  348 + result.append("\n"+gammalib::parformat("Energy range"));
  349 + if (m_energies.size() > 0) {
  350 + result.append(gammalib::str(m_energies[0])+ " - ");
  351 + result.append(gammalib::str(m_energies[m_energies.size()-1])+ " MeV");
  352 + }
  353 + else {
  354 + result.append("not defined");
  355 + }
  356 + result.append("\n"+gammalib::parformat("Entries"));
  357 + result.append(gammalib::str(m_energies.size()));
  358 +
  359 + // Append calibration database
  360 + result.append("\n"+m_caldb.print(chatter));
  361 +
  362 + // Append information
  363 +
  364 + } // endif: chatter was not silent
  365 +
  366 + // Return result
  367 + return result;
  368 +}
  369 +
  370 +
  371 +/*==========================================================================
  372 + = =
  373 + = Private methods =
  374 + = =
  375 + ==========================================================================*/
  376 +
  377 +/***********************************************************************//**
  378 + * @brief Initialise class members
  379 + ***************************************************************************/
  380 +void GCOMD1Response::init_members(void)
  381 +{
  382 + // Initialise members
  383 + m_caldb.clear();
  384 + m_energies.clear();
  385 + m_positions.clear();
  386 + m_sigmas.clear();
  387 + m_amplitudes.clear();
  388 + m_emins.clear();
  389 + m_ewidths.clear();
  390 + m_emaxs.clear();
  391 +
  392 + // Initialise pre-computation cache
  393 + m_energy = GEnergy(-1.0e30, "MeV"); // To assure initialisation
  394 + m_position = 0.0;
  395 + m_sigma = 0.0;
  396 + m_amplitude = 0.0;
  397 + m_emin = 0.0;
  398 + m_ewidth = 0.0;
  399 + m_emax = 0.0;
  400 +
  401 + // Return
  402 + return;
  403 +}
  404 +
  405 +
  406 +/***********************************************************************//**
  407 + * @brief Copy class members
  408 + *
  409 + * @param[in] rsp COMPTEL response.
  410 + ***************************************************************************/
  411 +void GCOMD1Response::copy_members(const GCOMD1Response& rsp)
  412 +{
  413 + // Copy attributes
  414 + m_caldb = rsp.m_caldb;
  415 + m_energies = rsp.m_energies;
  416 + m_positions = rsp.m_positions;
  417 + m_sigmas = rsp.m_sigmas;
  418 + m_amplitudes = rsp.m_amplitudes;
  419 + m_emins = rsp.m_emins;
  420 + m_ewidths = rsp.m_ewidths;
  421 + m_emaxs = rsp.m_emaxs;
  422 +
  423 + // Copy pre-computation cache
  424 + m_energy = rsp.m_energy;
  425 + m_position = rsp.m_position;
  426 + m_sigma = rsp.m_sigma;
  427 + m_amplitude = rsp.m_amplitude;
  428 + m_emin = rsp.m_emin;
  429 + m_ewidth = rsp.m_ewidth;
  430 + m_emax = rsp.m_emax;
  431 +
  432 + // Return
  433 + return;
  434 +}
  435 +
  436 +
  437 +/***********************************************************************//**
  438 + * @brief Delete class members
  439 + ***************************************************************************/
  440 +void GCOMD1Response::free_members(void)
  441 +{
  442 + // Return
  443 + return;
  444 +}
  445 +
  446 +
  447 +/***********************************************************************//**
  448 + * @brief Update computation cache
  449 + *
  450 + * @param[in] etrue True energy.
  451 + *
  452 + * The method assumes that there is a valid D1 module response.
  453 + ***************************************************************************/
  454 +void GCOMD1Response::update_cache(const GEnergy& etrue) const
  455 +{
  456 + // Update only if the true energy has changed
  457 + if (etrue != m_energy) {
  458 +
  459 + // Set true energy
  460 + m_energy = etrue;
  461 +
  462 + // Get true energy in MeV
  463 + double etrue_MeV = etrue.MeV();
  464 +
  465 + // If true energy is below lowest energy or above largest energy
  466 + // then set response to zero
  467 + if ((etrue_MeV < m_energies[0]) ||
  468 + (etrue_MeV > m_energies[m_energies.size()-1])) {
  469 + m_position = 0.0;
  470 + m_sigma = 0.0;
  471 + m_amplitude = 0.0;
  472 + m_emin = 0.0;
  473 + m_ewidth = 0.0;
  474 + m_emax = 0.0;
  475 + }
  476 +
  477 + // ... otherwise interpolate response parameters
  478 + else {
  479 +
  480 + // Interpolate response parameters
  481 + m_position = m_energies.interpolate(etrue_MeV, m_positions);
  482 + m_sigma = m_energies.interpolate(etrue_MeV, m_sigmas);
  483 + m_amplitude = m_energies.interpolate(etrue_MeV, m_amplitudes);
  484 + m_emin = m_energies.interpolate(etrue_MeV, m_emins);
  485 + m_ewidth = m_energies.interpolate(etrue_MeV, m_ewidths);
  486 + m_emax = m_energies.interpolate(etrue_MeV, m_emaxs);
  487 +
  488 + // Re-compute Gaussian amplitude to assure normalization
  489 + m_amplitude = 1.0 / (m_sigma * gammalib::sqrt_twopi);
  490 +
  491 + }
  492 +
  493 + } // endif: true energy has changed
  494 +
  495 + // Return
  496 + return;
  497 +}
... ...
inst/com/src/GCOMD2Response.cpp 0 → 100644
  View file @b6f5890
  1 +/***************************************************************************
  2 + * GCOMD2Response.cpp - COMPTEL D2 module response class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GCOMD2Response.cpp
  23 + * @brief COMPTEL D2 module response class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GMath.hpp"
  32 +#include "GCOMD2Response.hpp"
  33 +
  34 +/* __ Method name definitions ____________________________________________ */
  35 +
  36 +/* __ Macros _____________________________________________________________ */
  37 +
  38 +/* __ Coding definitions _________________________________________________ */
  39 +
  40 +/* __ Debug definitions __________________________________________________ */
  41 +
  42 +/* __ Constants __________________________________________________________ */
  43 +
  44 +
  45 +/*==========================================================================
  46 + = =
  47 + = Constructors/destructors =
  48 + = =
  49 + ==========================================================================*/
  50 +
  51 +/***********************************************************************//**
  52 + * @brief Void constructor
  53 + *
  54 + * Creates an empty COMPTEL D2 module response.
  55 + ***************************************************************************/
  56 +GCOMD2Response::GCOMD2Response(void)
  57 +{
  58 + // Initialise members
  59 + init_members();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Copy constructor
  68 + *
  69 + * @param[in] rsp COMPTEL D2 module response.
  70 + **************************************************************************/
  71 +GCOMD2Response::GCOMD2Response(const GCOMD2Response& rsp)
  72 +{
  73 + // Initialise members
  74 + init_members();
  75 +
  76 + // Copy members
  77 + copy_members(rsp);
  78 +
  79 + // Return
  80 + return;
  81 +}
  82 +
  83 +
  84 +/***********************************************************************//**
  85 + * @brief Response constructor
  86 + *
  87 + * @param[in] sdbname SDA response name.
  88 + * @param[in] caldb Calibration database.
  89 + *
  90 + * Create COMPTEL D2 module response response by loading an SDB file from a
  91 + * calibration database.
  92 + ***************************************************************************/
  93 +GCOMD2Response::GCOMD2Response(const std::string& sdbname, const GCaldb& caldb)
  94 +{
  95 + // Initialise members
  96 + init_members();
  97 +
  98 + // Set calibration database
  99 + this->caldb(caldb);
  100 +
  101 + // Load D2 module response
  102 + this->load(sdbname);
  103 +
  104 + // Return
  105 + return;
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Destructor
  111 + *
  112 + * Destroys instance of COMPTEL response object.
  113 + ***************************************************************************/
  114 +GCOMD2Response::~GCOMD2Response(void)
  115 +{
  116 + // Free members
  117 + free_members();
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/*==========================================================================
  125 + = =
  126 + = Operators =
  127 + = =
  128 + ==========================================================================*/
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Assignment operator
  132 + *
  133 + * @param[in] rsp COMPTEL D2 module response.
  134 + * @return COMPTEL D2 module response.
  135 + ***************************************************************************/
  136 +GCOMD2Response& GCOMD2Response::operator=(const GCOMD2Response& rsp)
  137 +{
  138 + // Execute only if object is not identical
  139 + if (this != &rsp) {
  140 +
  141 + // Free members
  142 + free_members();
  143 +
  144 + // Initialise members
  145 + init_members();
  146 +
  147 + // Copy members
  148 + copy_members(rsp);
  149 +
  150 + } // endif: object was not identical
  151 +
  152 + // Return this object
  153 + return *this;
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief D2 module response evaluation operator
  159 + *
  160 + * @param[in] etrue True energy.
  161 + * @param[in] ereco Reconstructed energy.
  162 + * @return COMPTEL D2 module response.
  163 + ***************************************************************************/
  164 +double GCOMD2Response::operator()(const GEnergy& etrue, const GEnergy& ereco) const
  165 +{
  166 + // Initialise response with zero
  167 + double response = 0.0;
  168 +
  169 + // Continue only if a response was loaded
  170 + if (!m_energies.is_empty()) {
  171 +
  172 + // Update response evaluation cache
  173 + update_cache(etrue);
  174 +
  175 + // Continue only if amplitude is positive
  176 + if (m_amplitude > 0.0) {
  177 +
  178 + // Compute D2 module response (here is where the real magic
  179 + // happens)
  180 + //TODO: Implement the correct formula
  181 + double arg = (m_energy.MeV()-ereco.MeV()) / m_sigma;
  182 + response = m_amplitude * std::exp(-0.5 * arg * arg);
  183 +
  184 + }
  185 +
  186 + } // endif: response was loaded
  187 +
  188 + // Return response
  189 + return response;
  190 +}
  191 +
  192 +
  193 +/*==========================================================================
  194 + = =
  195 + = Public methods =
  196 + = =
  197 + ==========================================================================*/
  198 +
  199 +/***********************************************************************//**
  200 + * @brief Clear instance
  201 + *
  202 + * Clears COMPTEL D2 module response object by resetting all members to an
  203 + * initial state. Any information that was present in the object before will
  204 + * be lost.
  205 + ***************************************************************************/
  206 +void GCOMD2Response::clear(void)
  207 +{
  208 + // Free class members
  209 + free_members();
  210 +
  211 + // Initialise members
  212 + init_members();
  213 +
  214 + // Return
  215 + return;
  216 +}
  217 +
  218 +
  219 +/***********************************************************************//**
  220 + * @brief Clone instance
  221 + *
  222 + * @return Pointer to deep copy of COMPTEL D2 module response.
  223 + ***************************************************************************/
  224 +GCOMD2Response* GCOMD2Response::clone(void) const
  225 +{
  226 + return new GCOMD2Response(*this);
  227 +}
  228 +
  229 +
  230 +/***********************************************************************//**
  231 + * @brief Load COMPTEL D2 module response.
  232 + *
  233 + * @param[in] sdbname COMPTEL D2 module response name.
  234 + *
  235 + * Loads the COMPTEL D2 module response with specified name @p sdbname. The
  236 + * method first searchs for an appropriate response in the calibration
  237 + * database. If no appropriate response is found, the method takes the
  238 + * database root path and response name to build the full path to the
  239 + * response file, and tries to load the response from these paths.
  240 + ***************************************************************************/
  241 +void GCOMD2Response::load(const std::string& sdbname)
  242 +{
  243 + // Clear instance but conserve calibration database
  244 + GCaldb caldb = m_caldb;
  245 + clear();
  246 + m_caldb = caldb;
  247 +
  248 + // First attempt reading the response using the GCaldb interface
  249 + GFilename filename = m_caldb.filename("","","SDB","","",sdbname);
  250 +
  251 + // If filename is empty then build filename from CALDB root path and
  252 + // response name
  253 + if (filename.is_empty()) {
  254 + filename = gammalib::filepath(m_caldb.rootdir(), sdbname);
  255 + if (!filename.exists()) {
  256 + GFilename testname = filename + ".fits";
  257 + if (testname.exists()) {
  258 + filename = testname;
  259 + }
  260 + }
  261 + }
  262 +
  263 + // Open FITS file
  264 + GFits fits(filename);
  265 +
  266 + // Get SDB table
  267 + const GFitsTable& sdb = *fits.table(1);
  268 +
  269 + // Read SDB
  270 + read(sdb);
  271 +
  272 + // Close SDB FITS file
  273 + fits.close();
  274 +
  275 + // Return
  276 + return;
  277 +}
  278 +
  279 +
  280 +/***********************************************************************//**
  281 + * @brief Read COMPTEL D2 module response.
  282 + *
  283 + * @param[in] table FITS table.
  284 + *
  285 + * Read the COMPTEL D2 module response from a SDB FITS table.
  286 + ***************************************************************************/
  287 +void GCOMD2Response::read(const GFitsTable& table)
  288 +{
  289 + // Initialise COMPTEL D2 module response vectors
  290 + m_energies.clear();
  291 + m_positions.clear();
  292 + m_sigmas.clear();
  293 + m_amplitudes.clear();
  294 + m_escapes1.clear();
  295 + m_escapes2.clear();
  296 + m_tails.clear();
  297 + m_backgrounds.clear();
  298 + m_emins.clear();
  299 + m_ewidths.clear();
  300 + m_emaxs.clear();
  301 +
  302 + // Extract number of entries in table
  303 + int num = table.nrows();
  304 +
  305 + // If there are entries then read them
  306 + if (num > 0) {
  307 +
  308 + // Get column pointers
  309 + const GFitsTableCol* ptr_energy = table["ENERGY"];
  310 + const GFitsTableCol* ptr_position = table["POSITION"];
  311 + const GFitsTableCol* ptr_sigma = table["WIDTH"];
  312 + const GFitsTableCol* ptr_amplitude = table["AMPLITUDE"];
  313 + const GFitsTableCol* ptr_escape1 = table["ESCAPE1"];
  314 + const GFitsTableCol* ptr_escape2 = table["ESCAPE2"];
  315 + const GFitsTableCol* ptr_tail = table["TAIL"];
  316 + const GFitsTableCol* ptr_background = table["BACKGROUND"];
  317 + const GFitsTableCol* ptr_emin = table["EMIN"];
  318 + const GFitsTableCol* ptr_ewidth = table["EWIDTH"];
  319 + const GFitsTableCol* ptr_emax = table["EMAX"];
  320 +
  321 + // Copy data from table into vectors
  322 + for (int i = 0; i < num; ++i) {
  323 + m_energies.append(ptr_energy->real(i));
  324 + m_positions.push_back(ptr_position->real(i));
  325 + m_sigmas.push_back(ptr_sigma->real(i));
  326 + m_amplitudes.push_back(ptr_amplitude->real(i));
  327 + m_escapes1.push_back(ptr_escape1->real(i));
  328 + m_escapes2.push_back(ptr_escape2->real(i));
  329 + m_tails.push_back(ptr_tail->real(i));
  330 + m_backgrounds.push_back(ptr_background->real(i));
  331 + m_emins.push_back(ptr_emin->real(i));
  332 + m_ewidths.push_back(ptr_ewidth->real(i));
  333 + m_emaxs.push_back(ptr_emax->real(i));
  334 + }
  335 +
  336 + } // endif: there were entries
  337 +
  338 + // Return
  339 + return;
  340 +}
  341 +
  342 +
  343 +/***********************************************************************//**
  344 + * @brief Print COMPTEL D2 module response information
  345 + *
  346 + * @param[in] chatter Chattiness.
  347 + * @return String containing COMPTEL D2 module response information.
  348 + ***************************************************************************/
  349 +std::string GCOMD2Response::print(const GChatter& chatter) const
  350 +{
  351 + // Initialise result string
  352 + std::string result;
  353 +
  354 + // Continue only if chatter is not silent
  355 + if (chatter != SILENT) {
  356 +
  357 + // Append header
  358 + result.append("=== GCOMD2Response ===");
  359 +
  360 + // Append D2 module response information
  361 + result.append("\n"+gammalib::parformat("Energy range"));
  362 + if (m_energies.size() > 0) {
  363 + result.append(gammalib::str(m_energies[0])+ " - ");
  364 + result.append(gammalib::str(m_energies[m_energies.size()-1])+ " MeV");
  365 + }
  366 + else {
  367 + result.append("not defined");
  368 + }
  369 + result.append("\n"+gammalib::parformat("Entries"));
  370 + result.append(gammalib::str(m_energies.size()));
  371 +
  372 + // Append calibration database
  373 + result.append("\n"+m_caldb.print(chatter));
  374 +
  375 + // Append information
  376 +
  377 + } // endif: chatter was not silent
  378 +
  379 + // Return result
  380 + return result;
  381 +}
  382 +
  383 +
  384 +/*==========================================================================
  385 + = =
  386 + = Private methods =
  387 + = =
  388 + ==========================================================================*/
  389 +
  390 +/***********************************************************************//**
  391 + * @brief Initialise class members
  392 + ***************************************************************************/
  393 +void GCOMD2Response::init_members(void)
  394 +{
  395 + // Initialise members
  396 + m_caldb.clear();
  397 + m_energies.clear();
  398 + m_positions.clear();
  399 + m_sigmas.clear();
  400 + m_amplitudes.clear();
  401 + m_escapes1.clear();
  402 + m_escapes2.clear();
  403 + m_tails.clear();
  404 + m_backgrounds.clear();
  405 + m_emins.clear();
  406 + m_ewidths.clear();
  407 + m_emaxs.clear();
  408 +
  409 + // Initialise pre-computation cache
  410 + m_energy = GEnergy(-1.0e30, "MeV"); // To assure initialisation
  411 + m_position = 0.0;
  412 + m_sigma = 0.0;
  413 + m_amplitude = 0.0;
  414 + m_escape1 = 0.0;
  415 + m_escape2 = 0.0;
  416 + m_tail = 0.0;
  417 + m_background = 0.0;
  418 + m_emin = 0.0;
  419 + m_ewidth = 0.0;
  420 + m_emax = 0.0;
  421 +
  422 + // Return
  423 + return;
  424 +}
  425 +
  426 +
  427 +/***********************************************************************//**
  428 + * @brief Copy class members
  429 + *
  430 + * @param[in] rsp COMPTEL response.
  431 + ***************************************************************************/
  432 +void GCOMD2Response::copy_members(const GCOMD2Response& rsp)
  433 +{
  434 + // Copy attributes
  435 + m_caldb = rsp.m_caldb;
  436 + m_energies = rsp.m_energies;
  437 + m_positions = rsp.m_positions;
  438 + m_sigmas = rsp.m_sigmas;
  439 + m_amplitudes = rsp.m_amplitudes;
  440 + m_escapes1 = rsp.m_escapes1;
  441 + m_escapes2 = rsp.m_escapes2;
  442 + m_tails = rsp.m_tails;
  443 + m_backgrounds = rsp.m_backgrounds;
  444 + m_emins = rsp.m_emins;
  445 + m_ewidths = rsp.m_ewidths;
  446 + m_emaxs = rsp.m_emaxs;
  447 +
  448 + // Copy pre-computation cache
  449 + m_energy = rsp.m_energy;
  450 + m_position = rsp.m_position;
  451 + m_sigma = rsp.m_sigma;
  452 + m_amplitude = rsp.m_amplitude;
  453 + m_escape1 = rsp.m_escape1;
  454 + m_escape2 = rsp.m_escape2;
  455 + m_tail = rsp.m_tail;
  456 + m_background = rsp.m_background;
  457 + m_emin = rsp.m_emin;
  458 + m_ewidth = rsp.m_ewidth;
  459 + m_emax = rsp.m_emax;
  460 +
  461 + // Return
  462 + return;
  463 +}
  464 +
  465 +
  466 +/***********************************************************************//**
  467 + * @brief Delete class members
  468 + ***************************************************************************/
  469 +void GCOMD2Response::free_members(void)
  470 +{
  471 + // Return
  472 + return;
  473 +}
  474 +
  475 +
  476 +/***********************************************************************//**
  477 + * @brief Update computation cache
  478 + *
  479 + * @param[in] etrue True energy.
  480 + *
  481 + * The method assumes that there is a valid D2 module response.
  482 + ***************************************************************************/
  483 +void GCOMD2Response::update_cache(const GEnergy& etrue) const
  484 +{
  485 + // Update only if the true energy has changed
  486 + if (etrue != m_energy) {
  487 +
  488 + // Set true energy
  489 + m_energy = etrue;
  490 +
  491 + // Get true energy in MeV
  492 + double etrue_MeV = etrue.MeV();
  493 +
  494 + // If true energy is below lowest energy or above largest energy
  495 + // then set response to zero
  496 + if ((etrue_MeV < m_energies[0]) ||
  497 + (etrue_MeV > m_energies[m_energies.size()-1])) {
  498 + m_position = 0.0;
  499 + m_sigma = 0.0;
  500 + m_amplitude = 0.0;
  501 + m_escape1 = 0.0;
  502 + m_escape2 = 0.0;
  503 + m_tail = 0.0;
  504 + m_background = 0.0;
  505 + m_emin = 0.0;
  506 + m_ewidth = 0.0;
  507 + m_emax = 0.0;
  508 + }
  509 +
  510 + // ... otherwise interpolate response parameters
  511 + else {
  512 +
  513 + // Interpolate response parameters
  514 + m_position = m_energies.interpolate(etrue_MeV, m_positions);
  515 + m_sigma = m_energies.interpolate(etrue_MeV, m_sigmas);
  516 + m_amplitude = m_energies.interpolate(etrue_MeV, m_amplitudes);
  517 + m_escape1 = m_energies.interpolate(etrue_MeV, m_escapes1);
  518 + m_escape2 = m_energies.interpolate(etrue_MeV, m_escapes2);
  519 + m_tail = m_energies.interpolate(etrue_MeV, m_tails);
  520 + m_background = m_energies.interpolate(etrue_MeV, m_backgrounds);
  521 + m_emin = m_energies.interpolate(etrue_MeV, m_emins);
  522 + m_ewidth = m_energies.interpolate(etrue_MeV, m_ewidths);
  523 + m_emax = m_energies.interpolate(etrue_MeV, m_emaxs);
  524 +
  525 + //TODO: Assure normalization
  526 +
  527 + }
  528 +
  529 + } // endif: true energy has changed
  530 +
  531 + // Return
  532 + return;
  533 +}
... ...
inst/com/src/GCOMResponse.cpp
  View file @b6f5890
... ... @@ -488,7 +488,7 @@ void GCOMResponse::load(const std::string&amp; rspname)
488 488 // Read IAQ
489 489 read(iaq);
490 490  
491   - // Close ARF FITS file
  491 + // Close IAQ FITS file
492 492 fits.close();
493 493  
494 494 // Return
... ... @@ -556,7 +556,7 @@ void GCOMResponse::read(const GFitsImage&amp; image)
556 556 /***********************************************************************//**
557 557 * @brief Print COMPTEL response information
558 558 *
559   - * @param[in] chatter Chattiness (defaults to NORMAL).
  559 + * @param[in] chatter Chattiness.
560 560 * @return String containing COMPTEL response information.
561 561 ***************************************************************************/
562 562 std::string GCOMResponse::print(const GChatter& chatter) const
... ...