Rolf Buehler / gammalib

Commit 8cd863ba2dafcf6f80c100e5f8a1a362949b81e7

Authored by Jürgen Knödlseder
1 parent 6265eb80

Add GCOMIaq class to implement IAQ computation (#2186) 

There is a first version that performs an IAQ computation, without however the correction coefficients.
Inline Side-by-side
Showing 13 changed files with 1077 additions and 73 deletions
inst/com/Makefile.am
  View file @8cd863b
... ... @@ -42,6 +42,7 @@ sources = src/GCOMSupport.cpp \
42 42 src/GCOMD1Response.cpp \
43 43 src/GCOMD2Response.cpp \
44 44 src/GCOMInstChars.cpp \
  45 + src/GCOMIaq.cpp \
45 46 src/GCOMObservation.cpp \
46 47 src/GCOMModelDRBFitting.cpp
47 48  
... ... @@ -53,6 +54,7 @@ pkginclude_HEADERS = include/GCOMEventCube.hpp \
53 54 include/GCOMD1Response.hpp \
54 55 include/GCOMD2Response.hpp \
55 56 include/GCOMInstChars.hpp \
  57 + include/GCOMIaq.hpp \
56 58 include/GCOMObservation.hpp \
57 59 include/GCOMModelDRBFitting.hpp \
58 60 include/GCOMLib.hpp
... ...
inst/com/include/GCOMD1Response.hpp
  View file @8cd863b
... ... @@ -31,7 +31,6 @@
31 31 #include <vector>
32 32 #include <string>
33 33 #include "GCaldb.hpp"
34   -#include "GEnergy.hpp"
35 34 #include "GNodeArray.hpp"
36 35  
37 36 /* __ Type definitions ___________________________________________________ */
... ... @@ -55,8 +54,8 @@ public:
55 54 ~GCOMD1Response(void);
56 55  
57 56 // Operators
58   - GCOMD1Response& operator=(const GCOMD1Response & rsp);
59   - double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  57 + GCOMD1Response& operator=(const GCOMD1Response& rsp);
  58 + double operator()(const double& etrue, const double& ereco) const;
60 59  
61 60 // Methods
62 61 void clear(void);
... ... @@ -66,6 +65,12 @@ public:
66 65 const GCaldb& caldb(void) const;
67 66 void load(const std::string& sdaname);
68 67 void read(const GFitsTable& table);
  68 + double position(const double& etrue) const;
  69 + double sigma(const double& etrue) const;
  70 + double amplitude(const double& etrue) const;
  71 + double emin(const double& etrue) const;
  72 + double ewidth(const double& etrue) const;
  73 + double emax(const double& etrue) const;
69 74 std::string print(const GChatter& chatter = NORMAL) const;
70 75  
71 76 private:
... ... @@ -73,7 +78,7 @@ private:
73 78 void init_members(void);
74 79 void copy_members(const GCOMD1Response& rsp);
75 80 void free_members(void);
76   - void update_cache(const GEnergy& etrue) const;
  81 + void update_cache(const double& etrue) const;
77 82  
78 83 // Private data members
79 84 GCaldb m_caldb; //!< Calibration database
... ... @@ -86,13 +91,13 @@ private:
86 91 std::vector<double> m_emaxs; //!< Upper energy limit of D1
87 92  
88 93 // Pre-computation cache
89   - mutable GEnergy m_energy;
90   - mutable double m_position;
91   - mutable double m_sigma;
92   - mutable double m_amplitude;
93   - mutable double m_emin;
94   - mutable double m_ewidth;
95   - mutable double m_emax;
  94 + mutable double m_energy;
  95 + mutable double m_position;
  96 + mutable double m_sigma;
  97 + mutable double m_amplitude;
  98 + mutable double m_emin;
  99 + mutable double m_ewidth;
  100 + mutable double m_emax;
96 101 };
97 102  
98 103  
... ... @@ -116,7 +121,7 @@ std::string GCOMD1Response::classname(void) const
116 121 inline
117 122 const GCaldb& GCOMD1Response::caldb(void) const
118 123 {
119   - return m_caldb;
  124 + return (m_caldb);
120 125 }
121 126  
122 127  
... ... @@ -134,4 +139,88 @@ void GCOMD1Response::caldb(const GCaldb&amp; caldb)
134 139 return;
135 140 }
136 141  
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Return photo peak position
  145 + *
  146 + * @param[in] etrue True energy (MeV).
  147 + * @return Photo peak position (MeV).
  148 + ***************************************************************************/
  149 +inline
  150 +double GCOMD1Response::position(const double& etrue) const
  151 +{
  152 + update_cache(etrue);
  153 + return (m_position);
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief Return photo peak standard deviation
  159 + *
  160 + * @param[in] etrue True energy (MeV).
  161 + * @return Photo peak standard deviation (MeV).
  162 + ***************************************************************************/
  163 +inline
  164 +double GCOMD1Response::sigma(const double& etrue) const
  165 +{
  166 + update_cache(etrue);
  167 + return (m_sigma);
  168 +}
  169 +
  170 +
  171 +/***********************************************************************//**
  172 + * @brief Return photo peak amplitude
  173 + *
  174 + * @param[in] etrue True energy (MeV).
  175 + * @return Photo peak amplitude.
  176 + ***************************************************************************/
  177 +inline
  178 +double GCOMD1Response::amplitude(const double& etrue) const
  179 +{
  180 + update_cache(etrue);
  181 + return (m_amplitude);
  182 +}
  183 +
  184 +
  185 +/***********************************************************************//**
  186 + * @brief Return minimum energy
  187 + *
  188 + * @param[in] etrue True energy (MeV).
  189 + * @return Minimum energy (MeV).
  190 + ***************************************************************************/
  191 +inline
  192 +double GCOMD1Response::emin(const double& etrue) const
  193 +{
  194 + update_cache(etrue);
  195 + return (m_emin);
  196 +}
  197 +
  198 +
  199 +/***********************************************************************//**
  200 + * @brief Return energy threshold width
  201 + *
  202 + * @param[in] etrue True energy (MeV).
  203 + * @return Energy threshold width (MeV).
  204 + ***************************************************************************/
  205 +inline
  206 +double GCOMD1Response::ewidth(const double& etrue) const
  207 +{
  208 + update_cache(etrue);
  209 + return (m_ewidth);
  210 +}
  211 +
  212 +
  213 +/***********************************************************************//**
  214 + * @brief Return maximum energy
  215 + *
  216 + * @param[in] etrue True energy (MeV).
  217 + * @return Maximum energy (MeV).
  218 + ***************************************************************************/
  219 +inline
  220 +double GCOMD1Response::emax(const double& etrue) const
  221 +{
  222 + update_cache(etrue);
  223 + return (m_emax);
  224 +}
  225 +
137 226 #endif /* GCOMD1RESPONSE_HPP */
... ...
inst/com/include/GCOMD2Response.hpp
  View file @8cd863b
... ... @@ -31,7 +31,6 @@
31 31 #include <vector>
32 32 #include <string>
33 33 #include "GCaldb.hpp"
34   -#include "GEnergy.hpp"
35 34 #include "GNodeArray.hpp"
36 35  
37 36 /* __ Type definitions ___________________________________________________ */
... ... @@ -55,8 +54,8 @@ public:
55 54 ~GCOMD2Response(void);
56 55  
57 56 // Operators
58   - GCOMD2Response& operator=(const GCOMD2Response & rsp);
59   - double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  57 + GCOMD2Response& operator=(const GCOMD2Response& rsp);
  58 + double operator()(const double& etrue, const double& ereco) const;
60 59  
61 60 // Methods
62 61 void clear(void);
... ... @@ -73,7 +72,7 @@ private:
73 72 void init_members(void);
74 73 void copy_members(const GCOMD2Response& rsp);
75 74 void free_members(void);
76   - void update_cache(const GEnergy& etrue) const;
  75 + void update_cache(const double& etrue) const;
77 76  
78 77 // Private data members
79 78 GCaldb m_caldb; //!< Calibration database
... ... @@ -90,17 +89,17 @@ private:
90 89 std::vector<double> m_emaxs; //!< Upper energy limit of D2
91 90  
92 91 // Pre-computation cache
93   - mutable GEnergy m_energy;
94   - mutable double m_position;
95   - mutable double m_sigma;
96   - mutable double m_amplitude;
97   - mutable double m_escape1;
98   - mutable double m_escape2;
99   - mutable double m_tail;
100   - mutable double m_background;
101   - mutable double m_emin;
102   - mutable double m_ewidth;
103   - mutable double m_emax;
  92 + mutable double 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;
104 103 };
105 104  
106 105  
... ...
inst/com/include/GCOMIaq.hpp 0 → 100644
  View file @8cd863b
  1 +/***************************************************************************
  2 + * GCOMIaq.hpp - COMPTEL instrument response representation *
  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 GCOMIaq.hpp
  23 + * @brief COMPTEL instrument response representation class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCOMIAQ_HPP
  28 +#define GCOMIAQ_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GMath.hpp"
  33 +#include "GEbounds.hpp"
  34 +#include "GFunction.hpp"
  35 +#include "GFitsImageFloat.hpp"
  36 +#include "GCOMD1Response.hpp"
  37 +#include "GCOMD2Response.hpp"
  38 +#include "GCOMInstChars.hpp"
  39 +
  40 +/* __ Type definitions ___________________________________________________ */
  41 +
  42 +/* __ Forward declarations _______________________________________________ */
  43 +class GFilename;
  44 +
  45 +
  46 +/***********************************************************************//**
  47 + * @class GCOMIaq
  48 + *
  49 + * @brief Interface for the COMPTEL instrument response representation class
  50 + ***************************************************************************/
  51 +class GCOMIaq : public GBase {
  52 +
  53 +public:
  54 + // Constructors and destructors
  55 + GCOMIaq(void);
  56 + GCOMIaq(const GCOMIaq& iaq);
  57 + GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
  58 + const double& phibar_max, const double& phibar_bin_size,
  59 + const GEbounds& ebounds);
  60 + ~GCOMIaq(void);
  61 +
  62 + // Operators
  63 + GCOMIaq& operator=(const GCOMIaq & iaq);
  64 +
  65 + // Methods
  66 + void clear(void);
  67 + GCOMIaq* clone(void) const;
  68 + std::string classname(void) const;
  69 + void save(const GFilename& filename, const bool& clobber) const;
  70 + std::string print(const GChatter& chatter = NORMAL) const;
  71 + void iaqwei(const GEnergy& energy, const double& weight);
  72 +
  73 +private:
  74 + // Private methods
  75 + void init_members(void);
  76 + void copy_members(const GCOMIaq& iaq);
  77 + void free_members(void);
  78 +
  79 + // RESPSI methods
  80 + //void iaqwei(const GEnergy& energy, const double& weight);
  81 + double respsc(const double& etrue1, const double& etrue2, const double& phibar);
  82 +
  83 + // Reponse integration kernel
  84 + class response_kernel : public GFunction {
  85 + public:
  86 + response_kernel(const GCOMD1Response& response_d1,
  87 + const GCOMD2Response& response_d2,
  88 + const double& etrue1,
  89 + const double& etrue2,
  90 + const double& phibar,
  91 + const double& etmin,
  92 + const double& etmax,
  93 + const double& e1min,
  94 + const double& e1max,
  95 + const double& e2min,
  96 + const double& e2max) :
  97 + m_response_d1(response_d1),
  98 + m_response_d2(response_d2),
  99 + m_etrue1(etrue1),
  100 + m_etrue2(etrue2),
  101 + m_cos_phibar(std::cos(phibar*gammalib::deg2rad)),
  102 + m_sin_phibar(std::sin(phibar*gammalib::deg2rad)),
  103 + m_etmin(etmin),
  104 + m_etmax(etmax),
  105 + m_e1min(e1min),
  106 + m_e1max(e1max),
  107 + m_e2min(e2min),
  108 + m_e2max(e2max) {}
  109 + double eval(const double& energy1);
  110 + protected:
  111 + const GCOMD1Response& m_response_d1; //!< Reference to D1 module response
  112 + const GCOMD2Response& m_response_d2; //!< Reference to D2 module response
  113 + double m_etrue1; //!< True D1 energy (MeV)
  114 + double m_etrue2; //!< True D2 energy (MeV)
  115 + double m_cos_phibar; //!< cos(phibar)
  116 + double m_sin_phibar; //!< sin(phibar)
  117 + double m_etmin; //!< Minimum total energy (MeV)
  118 + double m_etmax; //!< Maximum total energy (MeV)
  119 + double m_e1min; //!< Minimum D1 energy (MeV)
  120 + double m_e1max; //!< Maximum D1 energy (MeV)
  121 + double m_e2min; //!< Minimum D2 energy (MeV)
  122 + double m_e2max; //!< Maximum D2 energy (MeV)
  123 + };
  124 +
  125 + // Private data members
  126 + GFitsImageFloat m_iaq; //!< Response
  127 + GEbounds m_ebounds; //!< Energy boundaries
  128 + GCOMD1Response m_response_d1; //!< D1 module response
  129 + GCOMD2Response m_response_d2; //!< D2 module response
  130 + GCOMInstChars m_ict; //!< Instrument characteristics
  131 + double m_phigeo_max; //!< Maximum geometrical scatter angle (deg)
  132 + double m_phibar_max; //!< Maximum Compton scatter angle (deg)
  133 + double m_phigeo_bin_size; //!< Bin size in geometrical scatter angle (deg)
  134 + double m_phibar_bin_size; //!< Bin size in Compton scatter angle (deg)
  135 + double m_e1min; //!< Minimum D1 energy (MeV)
  136 + double m_e1max; //!< Maximum D1 energy (MeV)
  137 + double m_e2min; //!< Minimum D2 energy (MeV)
  138 + double m_e2max; //!< Maximum D2 energy (MeV)
  139 +};
  140 +
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Return class name
  144 + *
  145 + * @return String containing the class name ("GCOMIaq").
  146 + ***************************************************************************/
  147 +inline
  148 +std::string GCOMIaq::classname(void) const
  149 +{
  150 + return ("GCOMIaq");
  151 +}
  152 +
  153 +#endif /* GCOMIAQ_HPP */
... ...
inst/com/include/GCOMResponse.hpp
  View file @8cd863b
... ... @@ -56,9 +56,6 @@ public:
56 56 GCOMResponse(void);
57 57 GCOMResponse(const GCOMResponse& rsp);
58 58 GCOMResponse(const GCaldb& caldb, const std::string& rspname);
59   - //GCOMResponse(const int& phigeo_num, const int& phibar_num,
60   - // const double& phigeo_max, const double& phibar_max,
61   - // const double& phigeo_bin, const double& phibar_bin);
62 59 virtual ~GCOMResponse(void);
63 60  
64 61 // Operators
... ...
inst/com/pyext/GCOMD1Response.i
  View file @8cd863b
... ... @@ -44,7 +44,7 @@ public:
44 44 ~GCOMD1Response(void);
45 45  
46 46 // Operators
47   - double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  47 + double operator()(const double& etrue, const double& ereco) const;
48 48  
49 49 // Methods
50 50 void clear(void);
... ... @@ -54,6 +54,12 @@ public:
54 54 const GCaldb& caldb(void) const;
55 55 void load(const std::string& sdaname);
56 56 void read(const GFitsTable& hdu);
  57 + double position(const double& etrue) const;
  58 + double sigma(const double& etrue) const;
  59 + double amplitude(const double& etrue) const;
  60 + double emin(const double& etrue) const;
  61 + double ewidth(const double& etrue) const;
  62 + double emax(const double& etrue) const;
57 63 };
58 64  
59 65  
... ...
inst/com/pyext/GCOMD2Response.i
  View file @8cd863b
... ... @@ -44,7 +44,7 @@ public:
44 44 ~GCOMD2Response(void);
45 45  
46 46 // Operators
47   - double operator()(const GEnergy& etrue, const GEnergy& ereco) const;
  47 + double operator()(const double& etrue, const double& ereco) const;
48 48  
49 49 // Methods
50 50 void clear(void);
... ...
inst/com/pyext/GCOMIaq.i 0 → 100644
  View file @8cd863b
  1 +/***************************************************************************
  2 + * GCOMIaq.i - COMPTEL instrument response representation *
  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 GCOMIaq.i
  23 + * @brief COMPTEL instrument response representation class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCOMIaq.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCOMIaq
  34 + *
  35 + * @brief Interface for the COMPTEL instrument response representation class
  36 + ***************************************************************************/
  37 +class GCOMIaq : public GBase {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GCOMIaq(void);
  42 + GCOMIaq(const GCOMIaq& iaq);
  43 + GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
  44 + const double& phibar_max, const double& phibar_bin_size,
  45 + const GEbounds& ebounds);
  46 + ~GCOMIaq(void);
  47 +
  48 + // Methods
  49 + void clear(void);
  50 + GCOMIaq* clone(void) const;
  51 + std::string classname(void) const;
  52 + void save(const GFilename& filename, const bool& clobber) const;
  53 + void iaqwei(const GEnergy& energy, const double& weight);
  54 +};
  55 +
  56 +
  57 +/***********************************************************************//**
  58 + * @brief GCOMIaq class extension
  59 + ***************************************************************************/
  60 +%extend GCOMIaq {
  61 + GCOMIaq copy() {
  62 + return (*self);
  63 + }
  64 +};
... ...
inst/com/pyext/com.i
  View file @8cd863b
... ... @@ -70,3 +70,4 @@
70 70 %include "GCOMD1Response.i"
71 71 %include "GCOMD2Response.i"
72 72 %include "GCOMInstChar.i"
  73 +%include "GCOMIaq.i"
... ...
inst/com/src/GCOMD1Response.cpp
  View file @8cd863b
... ... @@ -157,11 +157,11 @@ GCOMD1Response&amp; GCOMD1Response::operator=(const GCOMD1Response&amp; rsp)
157 157 /***********************************************************************//**
158 158 * @brief D1 module response evaluation operator
159 159 *
160   - * @param[in] etrue True energy.
161   - * @param[in] ereco Reconstructed energy.
  160 + * @param[in] etrue True energy (MeV).
  161 + * @param[in] ereco Reconstructed energy (MeV).
162 162 * @return COMPTEL D1 module response.
163 163 ***************************************************************************/
164   -double GCOMD1Response::operator()(const GEnergy& etrue, const GEnergy& ereco) const
  164 +double GCOMD1Response::operator()(const double& etrue, const double& ereco) const
165 165 {
166 166 // Initialise response with zero
167 167 double response = 0.0;
... ... @@ -177,7 +177,7 @@ double GCOMD1Response::operator()(const GEnergy&amp; etrue, const GEnergy&amp; ereco) co
177 177  
178 178 // Compute D1 module response (here is where the real magic
179 179 // happens)
180   - double arg = (m_energy.MeV()-ereco.MeV()) / m_sigma;
  180 + double arg = (m_energy-ereco) / m_sigma;
181 181 response = m_amplitude * std::exp(-0.5 * arg * arg);
182 182  
183 183 }
... ... @@ -390,7 +390,7 @@ void GCOMD1Response::init_members(void)
390 390 m_emaxs.clear();
391 391  
392 392 // Initialise pre-computation cache
393   - m_energy = GEnergy(-1.0e30, "MeV"); // To assure initialisation
  393 + m_energy = 0.0;
394 394 m_position = 0.0;
395 395 m_sigma = 0.0;
396 396 m_amplitude = 0.0;
... ... @@ -447,11 +447,11 @@ void GCOMD1Response::free_members(void)
447 447 /***********************************************************************//**
448 448 * @brief Update computation cache
449 449 *
450   - * @param[in] etrue True energy.
  450 + * @param[in] etrue True energy (MeV).
451 451 *
452 452 * The method assumes that there is a valid D1 module response.
453 453 ***************************************************************************/
454   -void GCOMD1Response::update_cache(const GEnergy& etrue) const
  454 +void GCOMD1Response::update_cache(const double& etrue) const
455 455 {
456 456 // Update only if the true energy has changed
457 457 if (etrue != m_energy) {
... ... @@ -459,13 +459,10 @@ void GCOMD1Response::update_cache(const GEnergy&amp; etrue) const
459 459 // Set true energy
460 460 m_energy = etrue;
461 461  
462   - // Get true energy in MeV
463   - double etrue_MeV = etrue.MeV();
464   -
465 462 // If true energy is below lowest energy or above largest energy
466 463 // then set response to zero
467   - if ((etrue_MeV < m_energies[0]) ||
468   - (etrue_MeV > m_energies[m_energies.size()-1])) {
  464 + if ((etrue < m_energies[0]) ||
  465 + (etrue > m_energies[m_energies.size()-1])) {
469 466 m_position = 0.0;
470 467 m_sigma = 0.0;
471 468 m_amplitude = 0.0;
... ... @@ -478,12 +475,12 @@ void GCOMD1Response::update_cache(const GEnergy&amp; etrue) const
478 475 else {
479 476  
480 477 // 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);
  478 + m_position = m_energies.interpolate(etrue, m_positions);
  479 + m_sigma = m_energies.interpolate(etrue, m_sigmas);
  480 + m_amplitude = m_energies.interpolate(etrue, m_amplitudes);
  481 + m_emin = m_energies.interpolate(etrue, m_emins);
  482 + m_ewidth = m_energies.interpolate(etrue, m_ewidths);
  483 + m_emax = m_energies.interpolate(etrue, m_emaxs);
487 484  
488 485 // Re-compute Gaussian amplitude to assure normalization
489 486 m_amplitude = 1.0 / (m_sigma * gammalib::sqrt_twopi);
... ...
inst/com/src/GCOMD2Response.cpp
  View file @8cd863b
... ... @@ -157,11 +157,11 @@ GCOMD2Response&amp; GCOMD2Response::operator=(const GCOMD2Response&amp; rsp)
157 157 /***********************************************************************//**
158 158 * @brief D2 module response evaluation operator
159 159 *
160   - * @param[in] etrue True energy.
161   - * @param[in] ereco Reconstructed energy.
  160 + * @param[in] etrue True energy (MeV).
  161 + * @param[in] ereco Reconstructed energy (MeV).
162 162 * @return COMPTEL D2 module response.
163 163 ***************************************************************************/
164   -double GCOMD2Response::operator()(const GEnergy& etrue, const GEnergy& ereco) const
  164 +double GCOMD2Response::operator()(const double& etrue, const double& ereco) const
165 165 {
166 166 // Initialise response with zero
167 167 double response = 0.0;
... ... @@ -178,7 +178,7 @@ double GCOMD2Response::operator()(const GEnergy&amp; etrue, const GEnergy&amp; ereco) co
178 178 // Compute D2 module response (here is where the real magic
179 179 // happens)
180 180 //TODO: Implement the correct formula
181   - double arg = (m_energy.MeV()-ereco.MeV()) / m_sigma;
  181 + double arg = (m_energy-ereco) / m_sigma;
182 182 response = m_amplitude * std::exp(-0.5 * arg * arg);
183 183  
184 184 }
... ... @@ -407,7 +407,7 @@ void GCOMD2Response::init_members(void)
407 407 m_emaxs.clear();
408 408  
409 409 // Initialise pre-computation cache
410   - m_energy = GEnergy(-1.0e30, "MeV"); // To assure initialisation
  410 + m_energy = 0.0;
411 411 m_position = 0.0;
412 412 m_sigma = 0.0;
413 413 m_amplitude = 0.0;
... ... @@ -476,11 +476,11 @@ void GCOMD2Response::free_members(void)
476 476 /***********************************************************************//**
477 477 * @brief Update computation cache
478 478 *
479   - * @param[in] etrue True energy.
  479 + * @param[in] etrue True energy (MeV).
480 480 *
481 481 * The method assumes that there is a valid D2 module response.
482 482 ***************************************************************************/
483   -void GCOMD2Response::update_cache(const GEnergy& etrue) const
  483 +void GCOMD2Response::update_cache(const double& etrue) const
484 484 {
485 485 // Update only if the true energy has changed
486 486 if (etrue != m_energy) {
... ... @@ -488,13 +488,10 @@ void GCOMD2Response::update_cache(const GEnergy&amp; etrue) const
488 488 // Set true energy
489 489 m_energy = etrue;
490 490  
491   - // Get true energy in MeV
492   - double etrue_MeV = etrue.MeV();
493   -
494 491 // If true energy is below lowest energy or above largest energy
495 492 // then set response to zero
496   - if ((etrue_MeV < m_energies[0]) ||
497   - (etrue_MeV > m_energies[m_energies.size()-1])) {
  493 + if ((etrue < m_energies[0]) ||
  494 + (etrue > m_energies[m_energies.size()-1])) {
498 495 m_position = 0.0;
499 496 m_sigma = 0.0;
500 497 m_amplitude = 0.0;
... ... @@ -511,16 +508,16 @@ void GCOMD2Response::update_cache(const GEnergy&amp; etrue) const
511 508 else {
512 509  
513 510 // 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);
  511 + m_position = m_energies.interpolate(etrue, m_positions);
  512 + m_sigma = m_energies.interpolate(etrue, m_sigmas);
  513 + m_amplitude = m_energies.interpolate(etrue, m_amplitudes);
  514 + m_escape1 = m_energies.interpolate(etrue, m_escapes1);
  515 + m_escape2 = m_energies.interpolate(etrue, m_escapes2);
  516 + m_tail = m_energies.interpolate(etrue, m_tails);
  517 + m_background = m_energies.interpolate(etrue, m_backgrounds);
  518 + m_emin = m_energies.interpolate(etrue, m_emins);
  519 + m_ewidth = m_energies.interpolate(etrue, m_ewidths);
  520 + m_emax = m_energies.interpolate(etrue, m_emaxs);
524 521  
525 522 //TODO: Assure normalization
526 523  
... ...
inst/com/src/GCOMIaq.cpp 0 → 100644
  View file @8cd863b
  1 +/***************************************************************************
  2 + * GCOMIaq.cpp - COMPTEL instrument response representation *
  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 GCOMIaq.cpp
  23 + * @brief COMPTEL instrument response representation class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GFits.hpp"
  33 +#include "GFilename.hpp"
  34 +#include "GIntegral.hpp"
  35 +#include "GCOMIaq.hpp"
  36 +#include "GCOMSupport.hpp"
  37 +
  38 +/* __ Method name definitions ____________________________________________ */
  39 +
  40 +/* __ Macros _____________________________________________________________ */
  41 +
  42 +/* __ Coding definitions _________________________________________________ */
  43 +//#define G_RESPONSE_KERNEL_LIMITS
  44 +//#define G_RESPSC_NO_WARNINGS
  45 +
  46 +/* __ Debug definitions __________________________________________________ */
  47 +//#define G_DEBUG_IAQWEI
  48 +//#define G_DEBUG_RESPSC
  49 +//#define G_DEBUG_RESPONSE_KERNEL
  50 +
  51 +/* __ Constants __________________________________________________________ */
  52 +
  53 +
  54 +/*==========================================================================
  55 + = =
  56 + = Constructors/destructors =
  57 + = =
  58 + ==========================================================================*/
  59 +
  60 +/***********************************************************************//**
  61 + * @brief Void constructor
  62 + *
  63 + * Creates an empty COMPTEL instrument response representation.
  64 + ***************************************************************************/
  65 +GCOMIaq::GCOMIaq(void)
  66 +{
  67 + // Initialise members
  68 + init_members();
  69 +
  70 + // Return
  71 + return;
  72 +}
  73 +
  74 +
  75 +/***********************************************************************//**
  76 + * @brief Copy constructor
  77 + *
  78 + * @param[in] iaq COMPTEL instrument response representation.
  79 + **************************************************************************/
  80 +GCOMIaq::GCOMIaq(const GCOMIaq& iaq)
  81 +{
  82 + // Initialise members
  83 + init_members();
  84 +
  85 + // Copy members
  86 + copy_members(iaq);
  87 +
  88 + // Return
  89 + return;
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief COMPTEL instrument response representation constructor
  95 + *
  96 + * @param[in] phigeo_max Maximum geometrical scatter angle (deg).
  97 + * @param[in] phigeo_bin_size Bin size in geometrical scatter angle (deg).
  98 + * @param[in] phibar_max Maximum Compton scatter angle (deg).
  99 + * @param[in] phibar_bin_size Bin size in Compton scatter angle (deg).
  100 + * @param[in] ebounds Boundaries of observed energy.
  101 + *
  102 + * @todo Test input argument validity.
  103 + **************************************************************************/
  104 +GCOMIaq::GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
  105 + const double& phibar_max, const double& phibar_bin_size,
  106 + const GEbounds& ebounds)
  107 +{
  108 + // Initialise members
  109 + init_members();
  110 +
  111 + // Compute the number of bins in phigeo and phibar
  112 + int nphigeo = int(phigeo_max / phigeo_bin_size + 0.5);
  113 + int nphibar = int(phibar_max / phibar_bin_size + 0.5);
  114 +
  115 + // Store the bin definition. Use the number of bins to re-derive the
  116 + // maximum angles
  117 + m_phigeo_bin_size = phigeo_bin_size;
  118 + m_phibar_bin_size = phibar_bin_size;
  119 + m_phigeo_max = phigeo_bin_size * double(nphigeo);
  120 + m_phibar_max = phibar_bin_size * double(nphibar);
  121 +
  122 + // Store the energy boundaries
  123 + m_ebounds = ebounds;
  124 +
  125 + // Allocate FITS image
  126 + m_iaq = GFitsImageFloat(nphigeo, nphibar);
  127 +
  128 + // Set FITS image keywords
  129 + m_iaq.card("CTYPE1", "nphig", "Number of geometrical scatter angle pixels");
  130 + m_iaq.card("CRVAL1", 0.5*m_phigeo_bin_size, "[deg] Geometrical scatter angle for reference bin");
  131 + m_iaq.card("CRPIX1", 1.0, "Reference bin of geometrical scatter angle");
  132 + m_iaq.card("CDELT1", m_phigeo_bin_size, "[deg] Geometrical scatter angle bin size");
  133 + m_iaq.card("CTYPE2", "nphib", "Number of Compton scatter angle bins");
  134 + m_iaq.card("CRVAL2", 0.5*m_phibar_bin_size, "[deg] Compton scatter angle for reference bin");
  135 + m_iaq.card("CRPIX2", 1.0, "Reference bin of Compton scatter angle");
  136 + m_iaq.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
  137 + m_iaq.card("ENERGIE", 0.0, "[MeV] Source photon energy");
  138 + m_iaq.card("EMIN", m_ebounds.emin().MeV(), "[MeV] Minimum measured photon energy");
  139 + m_iaq.card("EMAX", m_ebounds.emax().MeV(), "[MeV] Maximum measured photon energy");
  140 + m_iaq.card("TELESCOP", "GRO", "Name of telescope");
  141 + m_iaq.card("INSTRUME", "COMPTEL", "Name of instrument");
  142 + m_iaq.card("ORIGIN", "GammaLib", "Origin of FITS file");
  143 + m_iaq.card("OBSERVER", "Unknown", "Observer that created FITS file");
  144 +
  145 + // Return
  146 + return;
  147 +}
  148 +
  149 +
  150 +/***********************************************************************//**
  151 + * @brief Destructor
  152 + *
  153 + * Destroys instance of COMPTEL response object.
  154 + ***************************************************************************/
  155 +GCOMIaq::~GCOMIaq(void)
  156 +{
  157 + // Free members
  158 + free_members();
  159 +
  160 + // Return
  161 + return;
  162 +}
  163 +
  164 +
  165 +/*==========================================================================
  166 + = =
  167 + = Operators =
  168 + = =
  169 + ==========================================================================*/
  170 +
  171 +/***********************************************************************//**
  172 + * @brief Assignment operator
  173 + *
  174 + * @param[in] iaq COMPTEL instrument response representation.
  175 + * @return COMPTEL instrument response representation.
  176 + ***************************************************************************/
  177 +GCOMIaq& GCOMIaq::operator=(const GCOMIaq& iaq)
  178 +{
  179 + // Execute only if object is not identical
  180 + if (this != &iaq) {
  181 +
  182 + // Free members
  183 + free_members();
  184 +
  185 + // Initialise members
  186 + init_members();
  187 +
  188 + // Copy members
  189 + copy_members(iaq);
  190 +
  191 + } // endif: object was not identical
  192 +
  193 + // Return this object
  194 + return *this;
  195 +}
  196 +
  197 +
  198 +/*==========================================================================
  199 + = =
  200 + = Public methods =
  201 + = =
  202 + ==========================================================================*/
  203 +
  204 +/***********************************************************************//**
  205 + * @brief Clear instance
  206 + *
  207 + * Clears COMPTEL instrument response representation by resetting all members
  208 + * to an initial state. Any information that was present in the object before
  209 + * will be lost.
  210 + ***************************************************************************/
  211 +void GCOMIaq::clear(void)
  212 +{
  213 + // Free class members
  214 + free_members();
  215 +
  216 + // Initialise members
  217 + init_members();
  218 +
  219 + // Return
  220 + return;
  221 +}
  222 +
  223 +
  224 +/***********************************************************************//**
  225 + * @brief Clone instance
  226 + *
  227 + * @return Pointer to deep copy of COMPTEL instrument response
  228 + * representation.
  229 + ***************************************************************************/
  230 +GCOMIaq* GCOMIaq::clone(void) const
  231 +{
  232 + return new GCOMIaq(*this);
  233 +}
  234 +
  235 +
  236 +/***********************************************************************//**
  237 + * @brief Save COMPTEL instrument response representation into FITS file
  238 + *
  239 + * @param[in] filename FITS file name.
  240 + * @param[in] clobber Overwrite existing file?
  241 + *
  242 + * Saves the COMPTEL instrument response representation into a FITS file. If
  243 + * the file exists already and the @p clobber parameter is true, the method
  244 + * will overwrite the content of the existing file. Otherwise, an exception
  245 + * is thrown.
  246 + ***************************************************************************/
  247 +void GCOMIaq::save(const GFilename& filename, const bool& clobber) const
  248 +{
  249 + // Work on a copy of the image to set the extension name
  250 + GFitsImageFloat iaq = m_iaq;
  251 +
  252 + // Set extension name
  253 + if (filename.has_extname()) {
  254 + iaq.extname(filename.extname());
  255 + }
  256 +
  257 + // Initialise empty FITS file
  258 + GFits fits;
  259 +
  260 + // Append IAQ to FITS file
  261 + fits.append(iaq);
  262 +
  263 + // Save FITS file (without extension name which was extracted earlier
  264 + // and set in the HDU)
  265 + fits.saveto(filename.url(), clobber);
  266 +
  267 + // Return
  268 + return;
  269 +}
  270 +
  271 +
  272 +/***********************************************************************//**
  273 + * @brief Print COMPTEL instrument response representation information
  274 + *
  275 + * @param[in] chatter Chattiness.
  276 + * @return String containing COMPTEL instrument response representation
  277 + * information.
  278 + ***************************************************************************/
  279 +std::string GCOMIaq::print(const GChatter& chatter) const
  280 +{
  281 + // Initialise result string
  282 + std::string result;
  283 +
  284 + // Continue only if chatter is not silent
  285 + if (chatter != SILENT) {
  286 +
  287 + // Append header
  288 + result.append("=== GCOMIaq ===");
  289 +
  290 + // Append information
  291 + result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
  292 + if (m_iaq.naxis() > 0) {
  293 + result.append(gammalib::str(m_iaq.naxes(0)));
  294 + }
  295 + else {
  296 + result.append("0");
  297 + }
  298 + result.append("\n"+gammalib::parformat("Number of Phibar bins"));
  299 + if (m_iaq.naxis() > 1) {
  300 + result.append(gammalib::str(m_iaq.naxes(1)));
  301 + }
  302 + else {
  303 + result.append("0");
  304 + }
  305 + result.append("\n"+gammalib::parformat("Maximum Phigeo value"));
  306 + result.append(gammalib::str(m_phigeo_max)+" deg");
  307 + result.append("\n"+gammalib::parformat("Maximum Phibar value"));
  308 + result.append(gammalib::str(m_phibar_max)+" deg");
  309 + result.append("\n"+gammalib::parformat("Phigeo bin size"));
  310 + result.append(gammalib::str(m_phigeo_bin_size)+" deg");
  311 + result.append("\n"+gammalib::parformat("Phibar bin size"));
  312 + result.append(gammalib::str(m_phibar_bin_size)+" deg");
  313 + result.append("\n"+gammalib::parformat("Energy range"));
  314 + result.append(m_ebounds.emin().print()+" - ");
  315 + result.append(m_ebounds.emax().print());
  316 +
  317 + } // endif: chatter was not silent
  318 +
  319 + // Return result
  320 + return result;
  321 +}
  322 +
  323 +
  324 +/*==========================================================================
  325 + = =
  326 + = Private methods =
  327 + = =
  328 + ==========================================================================*/
  329 +
  330 +/***********************************************************************//**
  331 + * @brief Initialise class members
  332 + ***************************************************************************/
  333 +void GCOMIaq::init_members(void)
  334 +{
  335 + // Initialise members
  336 + m_iaq.clear();
  337 + m_ebounds.clear();
  338 + m_response_d1.clear();
  339 + m_response_d2.clear();
  340 + m_ict.clear();
  341 + m_phigeo_max = 0.0;
  342 + m_phibar_max = 0.0;
  343 + m_phigeo_bin_size = 0.0;
  344 + m_phibar_bin_size = 0.0;
  345 + m_e1min = 0.070; //!< Default: 70 keV
  346 + m_e1max = 20.0; //!< Default: 20 MeV
  347 + m_e2min = 0.650; //!< Default: 650 keV
  348 + m_e2max = 30.0; //!< Default: 30 MeV
  349 +
  350 + // Return
  351 + return;
  352 +}
  353 +
  354 +
  355 +/***********************************************************************//**
  356 + * @brief Copy class members
  357 + *
  358 + * @param[in] iaq COMPTEL instrument response representation.
  359 + ***************************************************************************/
  360 +void GCOMIaq::copy_members(const GCOMIaq& iaq)
  361 +{
  362 + // Copy attributes
  363 + m_iaq = iaq.m_iaq;
  364 + m_ebounds = iaq.m_ebounds;
  365 + m_response_d1 = iaq.m_response_d1;
  366 + m_response_d2 = iaq.m_response_d2;
  367 + m_ict = iaq.m_ict;
  368 + m_phigeo_max = iaq.m_phigeo_max;
  369 + m_phibar_max = iaq.m_phibar_max;
  370 + m_phigeo_bin_size = iaq.m_phigeo_bin_size;
  371 + m_phibar_bin_size = iaq.m_phibar_bin_size;
  372 + m_e1min = iaq.m_e1min;
  373 + m_e1max = iaq.m_e1max;
  374 + m_e2min = iaq.m_e2min;
  375 + m_e2max = iaq.m_e2max;
  376 +
  377 + // Return
  378 + return;
  379 +}
  380 +
  381 +
  382 +/***********************************************************************//**
  383 + * @brief Delete class members
  384 + ***************************************************************************/
  385 +void GCOMIaq::free_members(void)
  386 +{
  387 + // Return
  388 + return;
  389 +}
  390 +
  391 +
  392 +/***********************************************************************//**
  393 + * @brief Add IAQ for a given photon energy
  394 + *
  395 + * @param[in] energy Input photon energy.
  396 + * @param[in] weight Spectral weight.
  397 + *
  398 + * The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F
  399 + * (release 1.0, 24-FEB-93).
  400 + ***************************************************************************/
  401 +void GCOMIaq::iaqwei(const GEnergy& energy, const double& weight)
  402 +{
  403 + // Get response information
  404 + GCaldb caldb("cgro","comptel");
  405 + m_response_d1 = GCOMD1Response(caldb, "DEFAULT");
  406 + m_response_d2 = GCOMD2Response(caldb, "DEFAULT");
  407 + m_ict = GCOMInstChars(caldb, "DEFAULT");
  408 +
  409 + // Get IAQ dimensions
  410 + int n_phigeo = m_iaq.naxes(0);
  411 + int n_phibar = m_iaq.naxes(1);
  412 +
  413 + // Set fine phibar step size to improve the accuracy
  414 + double finstp = 0.25;
  415 + int nfin = int(m_phibar_bin_size / finstp + 0.5);
  416 + finstp = m_phibar_bin_size / double(nfin);
  417 +
  418 + // Loop over phigeo
  419 + for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
  420 +
  421 + // Get geometrical scatter angle (deg)
  422 + double phigeo = (double(i_phigeo) + 0.5) * m_phigeo_bin_size;
  423 +
  424 + // Compute the true D1 and D2 energy deposits based on the
  425 + // geometrical scatter angle
  426 + double etrue2 = com_energy2(energy.MeV(), phigeo);
  427 + double etrue1 = energy.MeV() - etrue2;
  428 +
  429 + // Debug
  430 + #if defined(G_DEBUG_IAQWEI)
  431 + std::cout << "phigeo=" << phigeo;
  432 + std::cout << " etrue1=" << etrue1;
  433 + std::cout << " etrue2=" << etrue2 << std::endl;
  434 + double sum = 0.0;
  435 + #endif
  436 +
  437 + // Loop over phibar
  438 + for (int i_phibar = 0; i_phibar < n_phibar; ++i_phibar) {
  439 +
  440 + // Initialise start of phibar for this layer
  441 + double phibar = double(i_phibar) * m_phibar_bin_size + 0.5 * finstp;
  442 +
  443 + // Initialise response for this phibar layer
  444 + double response = 0.0;
  445 +
  446 + // Loop fine phibar layers
  447 + for (int i = 0; i < nfin; ++i, phibar += finstp) {
  448 +
  449 + // Compute response
  450 + response += respsc(etrue1, etrue2, phibar);
  451 +//std::cout << " phibar=" << phibar << " response=" << response << std::endl;
  452 +
  453 + } // endfor: looper over fine phibar layers
  454 +
  455 + // Multiply response by size of fine phibar layers
  456 + response *= finstp;
  457 +
  458 + // Store result
  459 + m_iaq(i_phigeo, i_phibar) = response;
  460 +
  461 + // Debug
  462 + #if defined(G_DEBUG_IAQWEI)
  463 + sum += response;
  464 + #endif
  465 +
  466 + } // endfor: looped over phibar
  467 +
  468 + // Debug
  469 + #if defined(G_DEBUG_IAQWEI)
  470 + std::cout << " Sum=" << sum << std::endl;
  471 + #endif
  472 +
  473 + } // endfor: looped over phigeo
  474 +
  475 + // Return
  476 + return;
  477 +}
  478 +
  479 +
  480 +/***********************************************************************//**
  481 + * @brief Integrates response over energy interval
  482 + *
  483 + * @param[in] etrue1 True D1 energy deposit (MeV).
  484 + * @param[in] etrue2 True D2 energy deposit (MeV).
  485 + * @param[in] phibar Compton scatter angle (deg).
  486 + * @param[in] e1min Minimum measured D1 energy deposit for integration (MeV).
  487 + * @param[in] e1max Maximum measured D1 energy deposit for integration (MeV).
  488 + *
  489 + * Compute the integral
  490 + *
  491 + * \f[
  492 + * R(\bar{\varphi}|\hat{E_1},\hat{E_2}) =
  493 + * \int_{E_1^{\rm min}}^{E_1^{\rm max}}
  494 + * R(E_1,E_2|\hat{E_1},\hat{E_2}) dE_1
  495 + * \f]
  496 + *
  497 + * where
  498 + * \f$R(E_1,E_2|\hat{E_1},\hat{E_2})\f$ is computed using
  499 + * GCOMIaq::response_kernel::eval,
  500 + *
  501 + * \f$E_1\f$ and \f$E_2\f$ are the measured D1 and D2 energy deposits,
  502 + * respectively,
  503 + * \f$\hat{E_1}\f$ and \f$\hat{E_2}\f$ are the true D1 and D2 energy
  504 + * deposits, respectively, and \f$\bar{\varphi}\f$ the Compton scatter
  505 + * angle.
  506 + *
  507 + * The code implementation is based on the COMPASS RESPSIT2 functions
  508 + * RESPSC.F (release 2.0, 14-Mar-91) and RESINT.F (release 3.0, 21-Oct-91).
  509 + * Since RESPSC.F only defines the threshold while RESINT.F does the real
  510 + * job, both functions were merged.
  511 + ***************************************************************************/
  512 +double GCOMIaq::respsc(const double& etrue1,
  513 + const double& etrue2,
  514 + const double& phibar)
  515 +{
  516 + // Initialise response
  517 + double response = 0.0;
  518 +
  519 + // Set integration limits
  520 + double e1min = m_response_d1.emin(etrue1) - 0.5 * m_response_d1.ewidth(etrue1);
  521 + double e1max = m_response_d1.emax(etrue1);
  522 +
  523 + // Set integration interval. Do not integrate more than 3 standard
  524 + // deviations away from photo peak.
  525 + double position = m_response_d1.position(etrue1);
  526 + double sigma = m_response_d1.sigma(etrue1);
  527 + double e_low = position - 3.0 * sigma;
  528 + double e_high = position + 3.0 * sigma;
  529 + double e_min = (e1min < e_low) ? e_low : e1min;
  530 + double e_max = (e1max > e_high) ? e_high : e1max;
  531 +
  532 + // Continue only if the integration interval is positive
  533 + if (e_min < e_max) {
  534 +
  535 + // Setup integration kernel
  536 + response_kernel integrand(m_response_d1,
  537 + m_response_d2,
  538 + etrue1,
  539 + etrue2,
  540 + phibar,
  541 + m_ebounds.emin().MeV(),
  542 + m_ebounds.emax().MeV(),
  543 + m_e1min,
  544 + m_e1max,
  545 + m_e2min,
  546 + m_e2max);
  547 +
  548 + // Setup integral
  549 + GIntegral integral(&integrand);
  550 +
  551 + // No warnings
  552 + #if defined(G_RESPSC_NO_WARNINGS)
  553 + integral.silent(true);
  554 + #endif
  555 +
  556 + // Perform integration
  557 + response = integral.romberg(e_min, e_max);
  558 +
  559 + } // endif: integration interval was positive
  560 +
  561 + // Debug
  562 + #if defined(G_DEBUG_RESPSC)
  563 + std::cout << " phibar=" << phibar;
  564 + std::cout << " e1min=" << e1min;
  565 + std::cout << " e1max=" << e1max;
  566 + std::cout << " e_min=" << e_min;
  567 + std::cout << " e_max=" << e_max;
  568 + std::cout << " response=" << response << std::endl;
  569 + #endif
  570 +
  571 + // Return response
  572 + return response;
  573 +}
  574 +
  575 +
  576 +/***********************************************************************//**
  577 + * @brief Computes product of D1 and D2 responses at energy E1 and phibar
  578 + *
  579 + * @param[in] energy1 D1 energy deposit (MeV).
  580 + *
  581 + * The method computes
  582 + *
  583 + * \f[
  584 + * R(E_1,E_2|\hat{E_1},\hat{E_2}) =
  585 + R_1(E_1|\hat{E_1}) \times R_2(E_2|\hat{E_2}) \times J
  586 + * \f]
  587 + *
  588 + * where
  589 + * \f$R_1(E_1|\hat{E_1})\f$ is the D1 module response function,
  590 + * \f$R_2(E_2|\hat{E_2})\f$ is the D2 module response function,
  591 + * \f$E_1\f$ and \f$E_2\f$ are the measured D1 and D2 energy deposits,
  592 + * respectively,
  593 + * \f$\hat{E_1}\f$ and \f$\hat{E_2}\f$ are the true D1 and D2 energy
  594 + * deposits, respectively, and \f$J\f$ is the Jacobian.
  595 + *
  596 + * The measured D2 energy deposit is computed using
  597 + *
  598 + * \f[
  599 + * E_2 = \frac{1}{2} \times \left(
  600 + * \sqrt{\frac{4 m_e c^2 E_1}{1 - \cos \bar{\varphi}} + E_1^2} - E_1
  601 + * \right)
  602 + * \f]
  603 + *
  604 + * where
  605 + * \f$m_e\,c^2\f$ the electron rest mass, and
  606 + * \f$\bar{\varphi}\f$ the Compton scatter angle.
  607 + *
  608 + * The Jacobian is given by
  609 + *
  610 + * \f[
  611 + * J = \frac{m_e c^2 \sqrt{E_1} \sin \bar{\varphi}}
  612 + * {(1 - \cos \bar{\varphi})^2
  613 + * \sqrt{\frac{4 m_e c^2}{1 - \cos \bar{\varphi}} + E_1}}
  614 + * \f]
  615 + *
  616 + * The code implementation is based on the COMPASS RESPSIT2 function FUNC.F
  617 + * (release 4.0, 21-Oct-91).
  618 + ***************************************************************************/
  619 +double GCOMIaq::response_kernel::eval(const double& energy1)
  620 +{
  621 + // Initialise result
  622 + double value = 0.0;
  623 +
  624 + // Set E1 and 1-cos(phibar). If the values are too small, limit to small
  625 + // positive numbers.
  626 + #if defined(G_RESPONSE_KERNEL_LIMITS)
  627 + double e1 = (energy1 < 1.0e-20) ? 1.0e-20 : energy1;
  628 + double cp = 1.0 - m_cos_phibar;
  629 + if (std::abs(cp) < 1.0e-20) {
  630 + cp = 1.0e-20;
  631 + }
  632 + #else
  633 + double e1 = energy1;
  634 + double cp = 1.0 - m_cos_phibar;
  635 + #endif
  636 +
  637 + // Compute E2 and Jacobian
  638 + double term1 = 4.0 * gammalib::mec2 / cp;
  639 + double e2 = 0.5 * (std::sqrt(term1 * e1 + e1*e1) - e1);
  640 + double jc = gammalib::mec2 * std::sqrt(e1) * m_sin_phibar /
  641 + (cp * cp * std::sqrt(term1 + e1));
  642 +
  643 + // Debug
  644 + #if defined(G_DEBUG_RESPONSE_KERNEL)
  645 + std::cout << " e1=" << e1;
  646 + std::cout << " e2=" << e2;
  647 + std::cout << " jc=" << jc << std::endl;
  648 + #endif
  649 +
  650 + // Allow to break now at any time
  651 + do {
  652 +
  653 + // Break if the D1 energy deposit is outside the energy range
  654 + if ((e1 < m_e1min) || (e1 > m_e1max)) {
  655 + break;
  656 + }
  657 +
  658 + // Break if the D2 energy deposit is outside the energy range
  659 + if ((e2 < m_e2min) || (e2 > m_e2max)) {
  660 + break;
  661 + }
  662 +
  663 + // Break if the total energy deposit is outside the energy range
  664 + double etot = e1 + e2;
  665 + if ((etot < m_etmin) || (etot > m_etmax)) {
  666 + break;
  667 + }
  668 +
  669 + // Break if the Jacobian precision is too poor
  670 + #if defined(G_RESPONSE_KERNEL_LIMITS)
  671 + if (jc < 1.0e-20) {
  672 + break;
  673 + }
  674 + #endif
  675 +
  676 + // Get D1 response. Break if it is too small.
  677 + double d1 = m_response_d1(m_etrue1, e1);
  678 + #if defined(G_RESPONSE_KERNEL_LIMITS)
  679 + if (d1 < 1.0e-20) {
  680 + break;
  681 + }
  682 + #endif
  683 +
  684 + // Get D2 response. Break if it is too small.
  685 + double d2 = m_response_d2(m_etrue2, e2);
  686 + #if defined(G_RESPONSE_KERNEL_LIMITS)
  687 + if (d2 < 1.0e-20) {
  688 + break;
  689 + }
  690 + #endif
  691 +
  692 + // Compute value
  693 + value = jc * d1 * d2;
  694 +
  695 + } while(false);
  696 +
  697 + // Return value
  698 + return value;
  699 +}
... ...
inst/com/src/GCOMInstChars.cpp
  View file @8cd863b
... ... @@ -503,7 +503,7 @@ double GCOMInstChars::multi_scatter(const GEnergy&amp; energy, const double&amp; phigeo)
503 503 * deposit in MeV. Coefficients are taken from Boron calibration data
504 504 * (ZA=1.14) to remain consistent with Rob van Dijk's SIMPSF corrections.
505 505 *
506   - * The code was implemented from the COMPASS RESPSIT2 function PSDACP.F
  506 + * The code implemented is based on the COMPASS RESPSIT2 function PSDACP.F
507 507 * (release 1.0, 11-DEC-92).
508 508 ***************************************************************************/
509 509 double GCOMInstChars::psd_correction(const GEnergy& energy, const double& phigeo) const
... ...