Leonardo Di Venere / gammalib

Commit bd7183cbf34e5de8b4ef1d687bddbe42a2d4b779

Authored by Jürgen Knödlseder
1 parent 9d9b7459

Remove GCTACubeSourcePoint and GCTACubeSourceDiffuse classes

Inline Side-by-side
Showing 5 changed files with 0 additions and 1239 deletions
inst/cta/include/GCTACubeSourceDiffuse.hpp deleted
View file @9d9b745
1 -/***************************************************************************  
2 - * GCTACubeSourceDiffuse.hpp - CTA cube analysis diffuse source class *  
3 - * ----------------------------------------------------------------------- *  
4 - * copyright (C) 2014-2016 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 GCTACubeSourceDiffuse.hpp  
23 - * @brief CTA cube analysis diffuse source class definition  
24 - * @author Juergen Knoedlseder  
25 - */  
26 -  
27 -#ifndef GCTACUBESOURCEDIFFUSE_HPP  
28 -#define GCTACUBESOURCEDIFFUSE_HPP  
29 -  
30 -/* __ Includes ___________________________________________________________ */  
31 -#include <string>  
32 -#include <vector>  
33 -#include "GSkyMap.hpp"  
34 -#include "GCTACubeSource.hpp"  
35 -#include "GCTAResponseCube.hpp"  
36 -  
37 -/* __ Type definitions ___________________________________________________ */  
38 -  
39 -/* __ Forward declarations _______________________________________________ */  
40 -class GModelSpatial;  
41 -class GObservation;  
42 -class GModelSpatialDiffuse;  
43 -class GSkyDir;  
44 -class GEnergy;  
45 -class GTime;  
46 -  
47 -  
48 -/***********************************************************************//**  
49 - * @class GCTACubeSourceDiffuse  
50 - *  
51 - * @brief CTA diffuse source cube class  
52 - *  
53 - * This class handles pre-computed response information for a diffuse source  
54 - * in a stacked cube analysis. It derives from the abstract GCTACubeSource  
55 - * class.  
56 - ***************************************************************************/  
57 -class GCTACubeSourceDiffuse : public GCTACubeSource {  
58 -  
59 -public:  
60 - // Constructors and destructors  
61 - GCTACubeSourceDiffuse(void);  
62 - GCTACubeSourceDiffuse(const GCTACubeSourceDiffuse& source);  
63 - virtual ~GCTACubeSourceDiffuse(void);  
64 -  
65 - // Operators  
66 - GCTACubeSourceDiffuse& operator=(const GCTACubeSourceDiffuse & source);  
67 -  
68 - // Implemented pure virtual methods  
69 - void clear(void);  
70 - GCTACubeSourceDiffuse* clone(void) const;  
71 - std::string classname(void) const;  
72 - GCTAClassCode code(void) const;  
73 - void set(const std::string& name,  
74 - const GModelSpatial& model,  
75 - const GObservation& obs);  
76 - std::string print(const GChatter& chatter = NORMAL) const;  
77 -  
78 - // Other methods  
79 - double par(const int& index) const;  
80 - double irf(const int& pixel, const int& iebin) const;  
81 - double irf(const int& pixel, const GEnergy& srcEng) const;  
82 - double psf(const GCTAResponseCube* rsp,  
83 - const GModelSpatial* model,  
84 - const GSkyDir& srcDir,  
85 - const GEnergy& srcEng,  
86 - const GTime& srcTime) const;  
87 -  
88 -protected:  
89 - // Protected methods  
90 - void init_members(void);  
91 - void copy_members(const GCTACubeSourceDiffuse& source);  
92 - void free_members(void);  
93 -  
94 - // Data members  
95 - GSkyMap m_cube; //!< Diffuse map convolved with IRF  
96 - GNodeArray m_logE; //!< Node array of log10 energy values in TeV  
97 -};  
98 -  
99 -  
100 -/***********************************************************************//**  
101 - * @brief Return class name  
102 - *  
103 - * @return String containing the class name ("GCTACubeSourceDiffuse").  
104 - ***************************************************************************/  
105 -inline  
106 -std::string GCTACubeSourceDiffuse::classname(void) const  
107 -{  
108 - return ("GCTACubeSourceDiffuse");  
109 -}  
110 -  
111 -  
112 -/***********************************************************************//**  
113 - * @brief Return class type code  
114 - *  
115 - * @return GCTA_SOURCE_CUBE_DIFFUSE.  
116 - *  
117 - * Returns the class type code GCTA_CUBE_SOURCE_DIFFUSE.  
118 - ***************************************************************************/  
119 -inline  
120 -GCTAClassCode GCTACubeSourceDiffuse::code(void) const  
121 -{  
122 - return (GCTA_CUBE_SOURCE_DIFFUSE);  
123 -}  
124 -  
125 -  
126 -/***********************************************************************//**  
127 - * @brief Return instrument response function  
128 - *  
129 - * @param[in] pixel Spatial pixel index [0,...,???]  
130 - * @param[in] iebin Energy index [0,...,???]  
131 - * @return Instrument response function.  
132 - *  
133 - * Returns the instrument response function.  
134 - ***************************************************************************/  
135 -inline  
136 -double GCTACubeSourceDiffuse::irf(const int& pixel, const int& iebin) const  
137 -{  
138 - return (m_cube(pixel, iebin));  
139 -}  
140 -  
141 -#endif /* GCTACUBESOURCEDIFFUSE_HPP */  
inst/cta/include/GCTACubeSourcePoint.hpp deleted
View file @9d9b745
1 -/***************************************************************************  
2 - * GCTACubeSourcePoint.hpp - CTA cube analysis point source class *  
3 - * ----------------------------------------------------------------------- *  
4 - * copyright (C) 2014-2015 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 GCTACubeSourcePoint.hpp  
23 - * @brief CTA cube analysis point source class definition  
24 - * @author Juergen Knoedlseder  
25 - */  
26 -  
27 -#ifndef GCTACUBESOURCEPOINT_HPP  
28 -#define GCTACUBESOURCEPOINT_HPP  
29 -  
30 -/* __ Includes ___________________________________________________________ */  
31 -#include <string>  
32 -#include <vector>  
33 -#include "GSkyDir.hpp"  
34 -#include "GCTACubeSource.hpp"  
35 -#include "GNodeArray.hpp"  
36 -  
37 -/* __ Type definitions ___________________________________________________ */  
38 -  
39 -/* __ Forward declarations _______________________________________________ */  
40 -class GModelSpatial;  
41 -class GObservation;  
42 -  
43 -  
44 -/***********************************************************************//**  
45 - * @class GCTACubeSourcePoint  
46 - *  
47 - * @brief CTA source cube base class  
48 - *  
49 - * This class handles pre-computed response information for a point source  
50 - * in a cube-style analysis. It derives from the abstract GCTACubeSource  
51 - * class.  
52 - ***************************************************************************/  
53 -class GCTACubeSourcePoint : public GCTACubeSource {  
54 -  
55 -public:  
56 - // Constructors and destructors  
57 - GCTACubeSourcePoint(void);  
58 - GCTACubeSourcePoint(const GCTACubeSourcePoint& source);  
59 - virtual ~GCTACubeSourcePoint(void);  
60 -  
61 - // Operators  
62 - GCTACubeSourcePoint& operator=(const GCTACubeSourcePoint & source);  
63 -  
64 - // Implemented pure virtual methods  
65 - void clear(void);  
66 - GCTACubeSourcePoint* clone(void) const;  
67 - std::string classname(void) const;  
68 - GCTAClassCode code(void) const;  
69 - void set(const std::string& name,  
70 - const GModelSpatial& model,  
71 - const GObservation& obs);  
72 - std::string print(const GChatter& chatter = NORMAL) const;  
73 -  
74 - // Other methods  
75 - double aeff(const int& index) const;  
76 - double delta(const int& index) const;  
77 - double psf(const int& ieng, const double& delta) const;  
78 - const GSkyDir& dir(void) const;  
79 -  
80 -protected:  
81 - // Protected methods  
82 - void init_members(void);  
83 - void copy_members(const GCTACubeSourcePoint& source);  
84 - void free_members(void);  
85 -  
86 - // Data members  
87 - GSkyDir m_dir; //!< Source direction  
88 - std::vector<double> m_aeff; //!< Deadtime corrected effective area (cm2)  
89 - std::vector<double> m_delta_map; //!< Distance of bin from source (radians)  
90 - std::vector<double> m_psf; //!< Point spread function  
91 - GNodeArray m_deltas; //!< Delta node array for PSF  
92 -};  
93 -  
94 -  
95 -/***********************************************************************//**  
96 - * @brief Return class name  
97 - *  
98 - * @return String containing the class name ("GCTACubeSourcePoint").  
99 - ***************************************************************************/  
100 -inline  
101 -std::string GCTACubeSourcePoint::classname(void) const  
102 -{  
103 - return ("GCTACubeSourcePoint");  
104 -}  
105 -  
106 -  
107 -/***********************************************************************//**  
108 - * @brief Return class type code  
109 - *  
110 - * @return GCTA_SOURCE_CUBE_POINT_SOURCE.  
111 - *  
112 - * Returns the class type code GCTA_CUBE_SOURCE_POINT.  
113 - ***************************************************************************/  
114 -inline  
115 -GCTAClassCode GCTACubeSourcePoint::code(void) const  
116 -{  
117 - return (GCTA_CUBE_SOURCE_POINT);  
118 -}  
119 -  
120 -  
121 -/***********************************************************************//**  
122 - * @brief Return deadtime corrected effective area in cm2  
123 - *  
124 - * @param[in] index Energy index [0,...,m_aeff.size()-1]  
125 - * @return Deadtime corrected effective area (cm2).  
126 - *  
127 - * Returns the deadtime corrected effective area.  
128 - ***************************************************************************/  
129 -inline  
130 -double GCTACubeSourcePoint::aeff(const int& index) const  
131 -{  
132 - double aeff = (index >= 0 && index < (int)m_aeff.size()) ? m_aeff[index] : 0.0;  
133 - return (aeff);  
134 -}  
135 -  
136 -  
137 -/***********************************************************************//**  
138 - * @brief Return angular distance to source in radians  
139 - *  
140 - * @param[in] index Spatial index [0,...,m_delta_map.size()-1]  
141 - * @return Angular distance to source (radians).  
142 - *  
143 - * Returns the angular seperation between a event cube pixel with @p index  
144 - * and the source position. If the index is not within the valid range, 99.9  
145 - * is returned.  
146 - ***************************************************************************/  
147 -inline  
148 -double GCTACubeSourcePoint::delta(const int& index) const  
149 -{  
150 - double delta = (index >= 0 && index < (int)m_delta_map.size())  
151 - ? m_delta_map[index] : 99.9;  
152 - return (delta);  
153 -}  
154 -  
155 -  
156 -/***********************************************************************//**  
157 - * @brief Return point source sky direction  
158 - *  
159 - * @return Point source sky direction.  
160 - ***************************************************************************/  
161 -inline  
162 -const GSkyDir& GCTACubeSourcePoint::dir(void) const  
163 -{  
164 - return (m_dir);  
165 -}  
166 -  
167 -#endif /* GCTACUBESOURCEPOINT_HPP */  
inst/cta/include/GCTALib.hpp
  View file @bd7183c
@@ -71,8 +71,6 @@ @@ -71,8 +71,6 @@
71 #include "GCTACubePsf.hpp" 71 #include "GCTACubePsf.hpp"
72 #include "GCTACubeEdisp.hpp" 72 #include "GCTACubeEdisp.hpp"
73 #include "GCTACubeSource.hpp" 73 #include "GCTACubeSource.hpp"
74 -#include "GCTACubeSourcePoint.hpp"  
75 -#include "GCTACubeSourceDiffuse.hpp"  
76 #include "GCTAModelBackground.hpp" 74 #include "GCTAModelBackground.hpp"
77 #include "GCTAModelSpatial.hpp" 75 #include "GCTAModelSpatial.hpp"
78 #include "GCTAModelSpatialGradient.hpp" 76 #include "GCTAModelSpatialGradient.hpp"
inst/cta/src/GCTACubeSourceDiffuse.cpp deleted
View file @9d9b745
1 -/***************************************************************************  
2 - * GCTACubeSourceDiffuse.cpp - CTA cube analysis diffuse source class *  
3 - * ----------------------------------------------------------------------- *  
4 - * copyright (C) 2014-2018 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 GCTACubeSourceDiffuse.cpp  
23 - * @brief CTA cube analysis diffuse source 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 "GCTACubeSourceDiffuse.hpp"  
33 -#include "GModelSpatialDiffuse.hpp"  
34 -#include "GObservation.hpp"  
35 -#include "GSkyDir.hpp"  
36 -#include "GEnergy.hpp"  
37 -#include "GTime.hpp"  
38 -#include "GIntegral.hpp"  
39 -#include "GCTAEventCube.hpp"  
40 -#include "GCTAResponseCube.hpp"  
41 -#include "GCTAResponse_helpers.hpp"  
42 -  
43 -/* __ OpenMP section _____________________________________________________ */  
44 -#ifdef _OPENMP  
45 -#include <omp.h>  
46 -#endif  
47 -  
48 -/* __ Method name definitions ____________________________________________ */  
49 -#define G_SET "GCTACubeSourceDiffuse::set(GModelSpatial&, GObservation&)"  
50 -  
51 -/* __ Macros _____________________________________________________________ */  
52 -  
53 -/* __ Coding definitions _________________________________________________ */  
54 -#define G_PSF_INTEGRATE //!< Integrate map over PSF (default)  
55 -  
56 -/* __ Debug definitions __________________________________________________ */  
57 -//#define G_DEBUG_SET  
58 -  
59 -/* __ Constants __________________________________________________________ */  
60 -  
61 -  
62 -/*==========================================================================  
63 - = =  
64 - = Constructors/destructors =  
65 - = =  
66 - ==========================================================================*/  
67 -  
68 -/***********************************************************************//**  
69 - * @brief Void constructor  
70 - ***************************************************************************/  
71 -GCTACubeSourceDiffuse::GCTACubeSourceDiffuse(void) : GCTACubeSource()  
72 -{  
73 - // Initialise class members  
74 - init_members();  
75 -  
76 - // Return  
77 - return;  
78 -}  
79 -  
80 -  
81 -/***********************************************************************//**  
82 - * @brief Copy constructor  
83 - *  
84 - * @param[in] source Diffuse source cube.  
85 - ***************************************************************************/  
86 -GCTACubeSourceDiffuse::GCTACubeSourceDiffuse(const GCTACubeSourceDiffuse& source) :  
87 - GCTACubeSource(source)  
88 -{  
89 - // Initialise class members  
90 - init_members();  
91 -  
92 - // Copy members  
93 - copy_members(source);  
94 -  
95 - // Return  
96 - return;  
97 -}  
98 -  
99 -  
100 -/***********************************************************************//**  
101 - * @brief Destructor  
102 - ***************************************************************************/  
103 -GCTACubeSourceDiffuse::~GCTACubeSourceDiffuse(void)  
104 -{  
105 - // Free members  
106 - free_members();  
107 -  
108 - // Return  
109 - return;  
110 -}  
111 -  
112 -  
113 -/*==========================================================================  
114 - = =  
115 - = Operators =  
116 - = =  
117 - ==========================================================================*/  
118 -  
119 -/***********************************************************************//**  
120 - * @brief Assignment operator  
121 - *  
122 - * @param[in] source Diffuse source cube.  
123 - * @return Diffuse source cube.  
124 - ***************************************************************************/  
125 -GCTACubeSourceDiffuse& GCTACubeSourceDiffuse::operator=(const GCTACubeSourceDiffuse& source)  
126 -{  
127 - // Execute only if object is not identical  
128 - if (this != &source) {  
129 -  
130 - // Copy base class members  
131 - this->GCTACubeSource::operator=(source);  
132 -  
133 - // Free members  
134 - free_members();  
135 -  
136 - // Initialise private members  
137 - init_members();  
138 -  
139 - // Copy members  
140 - copy_members(source);  
141 -  
142 - } // endif: object was not identical  
143 -  
144 - // Return this object  
145 - return *this;  
146 -}  
147 -  
148 -  
149 -/*==========================================================================  
150 - = =  
151 - = Public methods =  
152 - = =  
153 - ==========================================================================*/  
154 -  
155 -/***********************************************************************//**  
156 - * @brief Clear instance  
157 - *  
158 - * This method properly resets the object to an initial state.  
159 - ***************************************************************************/  
160 -void GCTACubeSourceDiffuse::clear(void)  
161 -{  
162 - // Free class members  
163 - free_members();  
164 - this->GCTACubeSource::free_members();  
165 -  
166 - // Initialise members  
167 - this->GCTACubeSource::init_members();  
168 - init_members();  
169 -  
170 - // Return  
171 - return;  
172 -}  
173 -  
174 -  
175 -/***********************************************************************//**  
176 - * @brief Clone instance  
177 - *  
178 - * @return Deep copy of diffuse source cube.  
179 - ***************************************************************************/  
180 -GCTACubeSourceDiffuse* GCTACubeSourceDiffuse::clone(void) const  
181 -{  
182 - // Return deep copy  
183 - return new GCTACubeSourceDiffuse(*this);  
184 -}  
185 -  
186 -  
187 -/***********************************************************************//**  
188 - * @brief Set diffuse source cube for a given observation  
189 - *  
190 - * @param[in] name Model name.  
191 - * @param[in] model Spatial model.  
192 - * @param[in] obs Observation.  
193 - *  
194 - * @exception GException::invalid_value  
195 - * Event or response cube not available.  
196 - * @exception GException::invalid_argument  
197 - * Spatial model is not of type GModelSpatialDiffuse.  
198 - *  
199 - * Sets the diffuse source cube assuming no energy dispersion and assuming  
200 - * a negligible variation of the effective area over the size of the  
201 - * point spread function.  
202 - ***************************************************************************/  
203 -void GCTACubeSourceDiffuse::set(const std::string& name,  
204 - const GModelSpatial& model,  
205 - const GObservation& obs)  
206 -{  
207 - // Debug option: initialise statistics  
208 - #if defined(G_DEBUG_SET)  
209 - int n_pixels_computed = 0;  
210 - std::cout << "GCTACubeSourceDiffuse::set entred." << std::endl;  
211 - #ifdef _OPENMP  
212 - double t_start = omp_get_wtime();  
213 - #else  
214 - clock_t t_start = clock();  
215 - #endif  
216 - #endif  
217 -  
218 - // Get pointer on CTA event cube  
219 - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());  
220 - if (cube == NULL) {  
221 - std::string msg = "Observation does not contain a CTA event cube.";  
222 - throw GException::invalid_value(G_SET, msg);  
223 - }  
224 -  
225 - // Get pointer on CTA response cube  
226 - const GCTAResponseCube* rsp = dynamic_cast<const GCTAResponseCube*>(obs.response());  
227 - if (rsp == NULL) {  
228 - std::string msg = "Observation does not contain a CTA response cube.";  
229 - throw GException::invalid_value(G_SET, msg);  
230 - }  
231 -  
232 - // Get pointer on diffuse model  
233 - const GModelSpatialDiffuse* spatial = dynamic_cast<const GModelSpatialDiffuse*>(&model);  
234 - if (spatial == NULL) {  
235 - std::string msg = "Spatial model is not of type GModelSpatialDiffuse.";  
236 - throw GException::invalid_argument(G_SET, msg);  
237 - }  
238 -  
239 - // Get diffuse source attributes  
240 - m_name = name;  
241 - GTime obsTime = cube->time();  
242 -  
243 - // Setup empty skymap  
244 - m_cube = cube->counts();  
245 - m_cube = 0.0;  
246 -  
247 - // Get energy boundaries and store them in node array  
248 - GEbounds ebounds = cube->ebounds();  
249 - for (int ieng = 0; ieng < ebounds.size(); ++ieng) {  
250 - m_logE.append(ebounds.elogmean(ieng).log10TeV());  
251 - }  
252 -  
253 - // Get livetime (in seconds) and deadtime correction factor  
254 - double livetime = rsp->exposure().livetime();  
255 - double deadc = rsp->exposure().deadc();  
256 -  
257 - // Get Psf radius (in degrees)  
258 - double delta_max = rsp->psf().delta_max() * gammalib::rad2deg;  
259 -  
260 - // Continue only if livetime is >0  
261 - if (livetime > 0.0) {  
262 -  
263 - // Loop over all spatial bins  
264 - for (int pixel = 0; pixel < cube->npix(); ++pixel) {  
265 -  
266 - // Get cube pixel sky direction  
267 - GSkyDir obsDir = cube->counts().inx2dir(pixel);  
268 -  
269 - // Continue only if model contains that sky direction  
270 - if (spatial->contains(obsDir, delta_max)) {  
271 -  
272 - // Loop over all energy layers  
273 - for (int iebin = 0; iebin < cube->ebins(); ++iebin) {  
274 -  
275 - // Get cube layer energy  
276 - const GEnergy& obsEng = cube->energy(iebin);  
277 -  
278 - // Determine exposure. We assume here that the exposure  
279 - // does not vary significantly over the PSF and just  
280 - // compute it at the pixel centre. We furthermore assume  
281 - // no energy dispersion, and thus compute exposure using  
282 - // the observed energy.  
283 - double aeff = rsp->exposure()(obsDir, obsEng);  
284 -  
285 - // Continue only if effective area is positive  
286 - if (aeff > 0.0) {  
287 -  
288 - // Recover effective area from exposure  
289 - aeff /= livetime;  
290 -  
291 - // Compute product of PSF and diffuse map, integrated  
292 - // over the relevant PSF area. We assume no energy  
293 - // dispersion and thus compute the product using the  
294 - // observed energy.  
295 - #if defined(G_PSF_INTEGRATE)  
296 - double psf = this->psf(rsp, &model, obsDir, obsEng, obsTime);  
297 - #else  
298 - double psf = model.eval(GPhoton(obsDir, obsEng, obsTime));  
299 - #endif  
300 -  
301 - // Set cube value  
302 - m_cube(pixel, iebin) = aeff * psf * deadc;  
303 -  
304 - } // endif: effective area was positive  
305 -  
306 - } // endfor: looped over all energy layers  
307 -  
308 - } // endif: pixel was contained in model  
309 -  
310 - // Debug option: update statistics  
311 - #if defined(G_DEBUG_SET)  
312 - n_pixels_computed++;  
313 - #endif  
314 -  
315 - } // endfor: looped over all spatial pixels  
316 -  
317 - } // endif: livetime was positive  
318 -  
319 - // Debug option: show statistics  
320 - #if defined(G_DEBUG_SET)  
321 - #ifdef _OPENMP  
322 - double t_elapse = omp_get_wtime()-t_start;  
323 - #else  
324 - double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;  
325 - #endif  
326 - std::cout << " Maximum theta ................: " << theta_max << " radians" << std::endl;  
327 - std::cout << " Number of spatial pixels used : " << n_pixels_computed << std::endl;  
328 - std::cout << " CPU usage ....................: " << t_elapse << " sec" << std::endl;  
329 - std::cout << "GCTACubeSourceDiffuse::set exit." << std::endl;  
330 - #endif  
331 -  
332 - // Return  
333 - return;  
334 -}  
335 -  
336 -  
337 -/***********************************************************************//**  
338 - * @brief Return instrument response function  
339 - *  
340 - * @param[in] pixel Spatial pixel index [0,...,???]  
341 - * @param[in] srcEng True photon energy  
342 - * @return Instrument response function.  
343 - *  
344 - * Returns the instrument response function for a given true photon energy  
345 - * @p srcEng.  
346 - ***************************************************************************/  
347 -double GCTACubeSourceDiffuse::irf(const int& pixel,  
348 - const GEnergy& srcEng) const  
349 -{  
350 - // Initialise response function  
351 - double irf = 0.0;  
352 -  
353 - // Continue only if there is energy information for the map cube  
354 - if (m_logE.size() > 0) {  
355 -  
356 - // Set energy for interpolation in log-energy  
357 - m_logE.set_value(srcEng.log10TeV());  
358 -  
359 - // Compute map cube intensity for the left and right map  
360 - double left_irf = m_cube(pixel, m_logE.inx_left());  
361 - double right_irf = m_cube(pixel, m_logE.inx_right());  
362 -  
363 - // Perform log-log interpolation  
364 - if (left_irf > 0.0 && right_irf > 0.0) {  
365 - double log_intensity = m_logE.wgt_left() * std::log(left_irf) +  
366 - m_logE.wgt_right() * std::log(right_irf);  
367 - irf = std::exp(log_intensity);  
368 - }  
369 - else if (left_irf > 0.0) {  
370 - irf = left_irf;  
371 - }  
372 - else if (right_irf > 0.0) {  
373 - irf = right_irf;  
374 - }  
375 -  
376 - } // endif: energy information was available  
377 -  
378 - // Return instrument response function  
379 - return irf;  
380 -}  
381 -  
382 -  
383 -/***********************************************************************//**  
384 - * @brief Integrate PSF over diffuse model  
385 - *  
386 - * @param[in] rsp Response function.  
387 - * @param[in] model Diffuse spatial model.  
388 - * @param[in] obsDir Observed event direction.  
389 - * @param[in] srcEng True photon energy.  
390 - * @param[in] srcTime True photon arrival time.  
391 - *  
392 - * Computes the integral  
393 - *  
394 - * \f[  
395 - * \int_0^{\delta_{\rm max}}  
396 - * {\rm PSF}(\delta) \times  
397 - * \int_0^{2\pi} {\rm Map}(\delta, \phi) \sin \delta  
398 - * {\rm d}\phi {\rm d}\delta  
399 - * \f]  
400 - *  
401 - * where \f${\rm Map}(\delta, \phi)\f$ is the diffuse map in the coordinate  
402 - * system of the point spread function, defined by the angle \f$\delta\f$  
403 - * between the true and the measured photon direction and the azimuth angle  
404 - * \f$\phi\f$ around the measured photon direction.  
405 - * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread  
406 - * function.  
407 - ***************************************************************************/  
408 -double GCTACubeSourceDiffuse::psf(const GCTAResponseCube* rsp,  
409 - const GModelSpatial* model,  
410 - const GSkyDir& obsDir,  
411 - const GEnergy& srcEng,  
412 - const GTime& srcTime) const  
413 -{  
414 - // Set minimum and maximum number of iterations for Romberg integration.  
415 - // These values have been determined after careful testing, see  
416 - // https://cta-redmine.irap.omp.eu/issues/2458  
417 - static const int min_iter_delta = 5;  
418 - static const int min_iter_phi = 5;  
419 - static const int max_iter_delta = 8;  
420 - static const int max_iter_phi = 8;  
421 -  
422 - // Initialise PSF  
423 - double psf = 0.0;  
424 -  
425 - // Compute rotation matrix to convert from PSF centred coordinate system  
426 - // spanned by delta and phi into the reference frame of the observed  
427 - // arrival direction given in Right Ascension and Declination.  
428 - GMatrix ry;  
429 - GMatrix rz;  
430 - ry.eulery(obsDir.dec_deg() - 90.0);  
431 - rz.eulerz(-obsDir.ra_deg());  
432 - GMatrix rot = (ry * rz).transpose();  
433 -  
434 - // Get offset angle integration interval in radians  
435 - double delta_min = 0.0;  
436 - double delta_max = 1.1 * rsp->psf().delta_max();  
437 -  
438 - // Get resolution of spatial model  
439 - double resolution = gammalib::resolution(model);  
440 -  
441 -  
442 - // Setup integration kernel. We take here the observed photon arrival  
443 - // direction as the true photon arrival direction because the PSF does  
444 - // not vary significantly over a small region.  
445 - cta_psf_diffuse_kern_delta integrand(rsp, model, &rot,  
446 - obsDir, srcEng, srcTime,  
447 - min_iter_phi, max_iter_phi,  
448 - resolution);  
449 -  
450 - // Set number of radial integration iterations  
451 - int iter = gammalib::iter_rho(delta_max, resolution,  
452 - min_iter_delta, max_iter_delta);  
453 -  
454 - // Setup integration  
455 - GIntegral integral(&integrand);  
456 -  
457 - // Set fixed number of iterations  
458 - integral.fixed_iter(iter);  
459 -  
460 - // Integrate over PSF delta angle  
461 - psf = integral.romberg(delta_min, delta_max);  
462 -  
463 - // Return PSF  
464 - return psf;  
465 -}  
466 -  
467 -  
468 -/***********************************************************************//**  
469 - * @brief Print diffuse source cube information  
470 - *  
471 - * @param[in] chatter Chattiness (defaults to NORMAL).  
472 - * @return String containing diffuse source cube information.  
473 - ***************************************************************************/  
474 -std::string GCTACubeSourceDiffuse::print(const GChatter& chatter) const  
475 -{  
476 - // Initialise result string  
477 - std::string result;  
478 -  
479 - // Continue only if chatter is not silent  
480 - if (chatter != SILENT) {  
481 -  
482 - // Append header  
483 - result.append("=== GCTACubeSourceDiffuse ===");  
484 - result.append("\n"+gammalib::parformat("Source name") + name());  
485 -  
486 - } // endif: chatter was not silent  
487 -  
488 - // Return result  
489 - return result;  
490 -}  
491 -  
492 -  
493 -/*==========================================================================  
494 - = =  
495 - = Private methods =  
496 - = =  
497 - ==========================================================================*/  
498 -  
499 -/***********************************************************************//**  
500 - * @brief Initialise class members  
501 - ***************************************************************************/  
502 -void GCTACubeSourceDiffuse::init_members(void)  
503 -{  
504 - // Initialise members  
505 - m_cube.clear();  
506 - m_logE.clear();  
507 -  
508 - // Return  
509 - return;  
510 -}  
511 -  
512 -  
513 -/***********************************************************************//**  
514 - * @brief Copy class members  
515 - *  
516 - * @param[in] source Point diffuse cube.  
517 - ***************************************************************************/  
518 -void GCTACubeSourceDiffuse::copy_members(const GCTACubeSourceDiffuse& source)  
519 -{  
520 - // Copy members  
521 - m_cube = source.m_cube;  
522 - m_logE = source.m_logE;  
523 -  
524 - // Return  
525 - return;  
526 -}  
527 -  
528 -/***********************************************************************//**  
529 - * @brief Delete class members  
530 - ***************************************************************************/  
531 -void GCTACubeSourceDiffuse::free_members(void)  
532 -{  
533 - // Return  
534 - return;  
535 -}  
inst/cta/src/GCTACubeSourcePoint.cpp deleted
View file @9d9b745
1 -/***************************************************************************  
2 - * GCTACubeSourcePoint.cpp - CTA cube analysis point source class *  
3 - * ----------------------------------------------------------------------- *  
4 - * copyright (C) 2014-2016 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 GCTACubeSourcePoint.cpp  
23 - * @brief CTA cube analysis point source class implementation  
24 - * @author Juergen Knoedlseder  
25 - */  
26 -  
27 -/* __ Includes ___________________________________________________________ */  
28 -#ifdef HAVE_CONFIG_H  
29 -#include <config.h>  
30 -#endif  
31 -#include "GCTACubeSourcePoint.hpp"  
32 -#include "GTools.hpp"  
33 -#include "GModelSpatial.hpp"  
34 -#include "GObservation.hpp"  
35 -#include "GModelSpatialPointSource.hpp"  
36 -#include "GCTAEventCube.hpp"  
37 -#include "GCTAResponseCube.hpp"  
38 -  
39 -/* __ Method name definitions ____________________________________________ */  
40 -#define G_SET "GCTACubeSourcePoint::set(GModelSpatial&, GObservation&)"  
41 -  
42 -/* __ Macros _____________________________________________________________ */  
43 -  
44 -/* __ Coding definitions _________________________________________________ */  
45 -  
46 -/* __ Debug definitions __________________________________________________ */  
47 -  
48 -/* __ Constants __________________________________________________________ */  
49 -  
50 -  
51 -/*==========================================================================  
52 - = =  
53 - = Constructors/destructors =  
54 - = =  
55 - ==========================================================================*/  
56 -  
57 -/***********************************************************************//**  
58 - * @brief Void constructor  
59 - ***************************************************************************/  
60 -GCTACubeSourcePoint::GCTACubeSourcePoint(void) : GCTACubeSource()  
61 -{  
62 - // Initialise class members  
63 - init_members();  
64 -  
65 - // Return  
66 - return;  
67 -}  
68 -  
69 -  
70 -/***********************************************************************//**  
71 - * @brief Copy constructor  
72 - *  
73 - * @param[in] source Point source cube.  
74 - ***************************************************************************/  
75 -GCTACubeSourcePoint::GCTACubeSourcePoint(const GCTACubeSourcePoint& source) :  
76 - GCTACubeSource(source)  
77 -{  
78 - // Initialise class members  
79 - init_members();  
80 -  
81 - // Copy members  
82 - copy_members(source);  
83 -  
84 - // Return  
85 - return;  
86 -}  
87 -  
88 -  
89 -/***********************************************************************//**  
90 - * @brief Destructor  
91 - ***************************************************************************/  
92 -GCTACubeSourcePoint::~GCTACubeSourcePoint(void)  
93 -{  
94 - // Free members  
95 - free_members();  
96 -  
97 - // Return  
98 - return;  
99 -}  
100 -  
101 -  
102 -/*==========================================================================  
103 - = =  
104 - = Operators =  
105 - = =  
106 - ==========================================================================*/  
107 -  
108 -/***********************************************************************//**  
109 - * @brief Assignment operator  
110 - *  
111 - * @param[in] source Point source cube.  
112 - * @return Point source cube.  
113 - ***************************************************************************/  
114 -GCTACubeSourcePoint& GCTACubeSourcePoint::operator=(const GCTACubeSourcePoint& source)  
115 -{  
116 - // Execute only if object is not identical  
117 - if (this != &source) {  
118 -  
119 - // Copy base class members  
120 - this->GCTACubeSource::operator=(source);  
121 -  
122 - // Free members  
123 - free_members();  
124 -  
125 - // Initialise private members  
126 - init_members();  
127 -  
128 - // Copy members  
129 - copy_members(source);  
130 -  
131 - } // endif: object was not identical  
132 -  
133 - // Return this object  
134 - return *this;  
135 -}  
136 -  
137 -  
138 -/*==========================================================================  
139 - = =  
140 - = Public methods =  
141 - = =  
142 - ==========================================================================*/  
143 -  
144 -/***********************************************************************//**  
145 - * @brief Clear instance  
146 - *  
147 - * This method properly resets the object to an initial state.  
148 - ***************************************************************************/  
149 -void GCTACubeSourcePoint::clear(void)  
150 -{  
151 - // Free class members  
152 - free_members();  
153 - this->GCTACubeSource::free_members();  
154 -  
155 - // Initialise members  
156 - this->GCTACubeSource::init_members();  
157 - init_members();  
158 -  
159 - // Return  
160 - return;  
161 -}  
162 -  
163 -  
164 -/***********************************************************************//**  
165 - * @brief Clone instance  
166 - *  
167 - * @return Deep copy of point source cube.  
168 - ***************************************************************************/  
169 -GCTACubeSourcePoint* GCTACubeSourcePoint::clone(void) const  
170 -{  
171 - // Return deep copy  
172 - return new GCTACubeSourcePoint(*this);  
173 -}  
174 -  
175 -  
176 -/***********************************************************************//**  
177 - * @brief Set point source cube for a given observation  
178 - *  
179 - * @param[in] name Model name.  
180 - * @param[in] model Spatial model.  
181 - * @param[in] obs Observation.  
182 - ***************************************************************************/  
183 -void GCTACubeSourcePoint::set(const std::string& name,  
184 - const GModelSpatial& model,  
185 - const GObservation& obs)  
186 -{  
187 - // Get pointer to model source model  
188 - const GModelSpatialPointSource* ptsrc = dynamic_cast<const GModelSpatialPointSource*>(&model);  
189 - if (ptsrc == NULL) {  
190 - std::string msg = "Model is not a spatial point source model.";  
191 - throw GException::invalid_value(G_SET, msg);  
192 - }  
193 -  
194 - // Get pointer on CTA event cube  
195 - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());  
196 - if (cube == NULL) {  
197 - std::string msg = "Observation does not contain a CTA event cube.";  
198 - throw GException::invalid_value(G_SET, msg);  
199 - }  
200 -  
201 - // Get pointer on CTA response cube  
202 - const GCTAResponseCube* rsp = dynamic_cast<const GCTAResponseCube*>(obs.response());  
203 - if (rsp == NULL) {  
204 - std::string msg = "Observation does not contain a CTA response cube.";  
205 - throw GException::invalid_value(G_SET, msg);  
206 - }  
207 -  
208 - // Get point source attributes  
209 - m_name = name;  
210 - m_dir = ptsrc->dir();  
211 - GTime srcTime = cube->time();  
212 -  
213 - // Set PSF deltas in radians  
214 - m_deltas = rsp->psf().deltas();  
215 - for (int i = 0; i < m_deltas.size(); ++i) {  
216 - m_deltas[i] *= gammalib::deg2rad;  
217 - }  
218 -  
219 - // Initialise data members  
220 - m_aeff.assign(cube->ebins(), 0.0);  
221 - m_delta_map.assign(cube->npix(), 0.0);  
222 - m_psf.assign(cube->ebins()*m_deltas.size(), 0.0);  
223 -  
224 - // Get livetime (in seconds) and deadtime correction factor  
225 - double livetime = rsp->exposure().livetime();  
226 - double deadc = rsp->exposure().deadc();  
227 -  
228 - // Compute deadtime corrected effective area  
229 - for (int i = 0; i < cube->ebins(); ++i) {  
230 -  
231 - // Initialise effective area  
232 - double aeff = 0.0;  
233 -  
234 - // Continue only if livetime is >0  
235 - if (livetime > 0.0) {  
236 -  
237 - // Get source energy  
238 - const GEnergy& srcEng = cube->energy(i);  
239 -  
240 - // Get exposure  
241 - aeff = rsp->exposure()(m_dir, srcEng);  
242 -  
243 - // Recover effective area from exposure  
244 - aeff /= livetime;  
245 -  
246 - // Apply deadtime correction  
247 - aeff *= deadc;  
248 -  
249 - } // endif: livetime was positive  
250 -  
251 - // Store deadtime corrected effective area in vector  
252 - m_aeff[i] = aeff;  
253 -  
254 - } // endfor: looped over energies  
255 -  
256 - // Compute distance map  
257 - for (int i = 0; i < cube->npix(); ++i) {  
258 -  
259 - // Get cube pixel sky direction  
260 - GSkyDir obsDir = cube->counts().inx2dir(i);  
261 -  
262 - // Determine angular distance between point source direction and  
263 - // cube pixel sky direction (radians)  
264 - double delta = m_dir.dist(obsDir);  
265 -  
266 - // Store angular distance in vector  
267 - m_delta_map[i] = delta;  
268 -  
269 - } // endfor: looped over spatial cube pixels  
270 -  
271 - // Compute point spread function  
272 - for (int i = 0, inx = 0; i < cube->ebins(); ++i) {  
273 -  
274 - // Get source energy  
275 - const GEnergy& srcEng = cube->energy(i);  
276 -  
277 - // Loop over delta bins  
278 - for (int k = 0; k < m_deltas.size(); ++k, ++inx) {  
279 -  
280 - // Get PSF  
281 - m_psf[inx] = rsp->psf()(m_dir, m_deltas[k], srcEng);  
282 -  
283 - } // endfor: looped over delta bins  
284 -  
285 - } // endfor: looped over PSF  
286 -  
287 - // Return  
288 - return;  
289 -}  
290 -  
291 -  
292 -/***********************************************************************//**  
293 - * @brief Set point source cube for a given observation  
294 - *  
295 - * @param[in] ieng Energy layer.  
296 - * @param[in] delta Distance between true and measured photon direction  
297 - * (radians).  
298 - ***************************************************************************/  
299 -double GCTACubeSourcePoint::psf(const int& ieng, const double& delta) const  
300 -{  
301 - // Set node array interpolation values  
302 - m_deltas.set_value(delta);  
303 -  
304 - // Compute offset  
305 - int offset = ieng * m_deltas.size();  
306 -  
307 - // Compute PSF by bi-linear interpolation  
308 - double psf = m_deltas.wgt_left() * m_psf[offset+m_deltas.inx_left()] +  
309 - m_deltas.wgt_right() * m_psf[offset+m_deltas.inx_right()];  
310 -  
311 - // Reset negative PSF values  
312 - if (psf < 0.0) {  
313 - psf = 0.0;  
314 - }  
315 -  
316 - // Return  
317 - return psf;  
318 -}  
319 -  
320 -  
321 -/***********************************************************************//**  
322 - * @brief Print point source cube information  
323 - *  
324 - * @param[in] chatter Chattiness (defaults to NORMAL).  
325 - * @return String containing point source cube information.  
326 - ***************************************************************************/  
327 -std::string GCTACubeSourcePoint::print(const GChatter& chatter) const  
328 -{  
329 - // Initialise result string  
330 - std::string result;  
331 -  
332 - // Continue only if chatter is not silent  
333 - if (chatter != SILENT) {  
334 -  
335 - // Append header  
336 - result.append("=== GCTACubeSourcePoint ===");  
337 - result.append("\n"+gammalib::parformat("Source name") + name());  
338 -  
339 - } // endif: chatter was not silent  
340 -  
341 - // Return result  
342 - return result;  
343 -}  
344 -  
345 -  
346 -/*==========================================================================  
347 - = =  
348 - = Private methods =  
349 - = =  
350 - ==========================================================================*/  
351 -  
352 -/***********************************************************************//**  
353 - * @brief Initialise class members  
354 - ***************************************************************************/  
355 -void GCTACubeSourcePoint::init_members(void)  
356 -{  
357 - // Initialise members  
358 - m_dir.clear();  
359 - m_aeff.clear();  
360 - m_delta_map.clear();  
361 - m_psf.clear();  
362 - m_deltas.clear();  
363 -  
364 - // Return  
365 - return;  
366 -}  
367 -  
368 -  
369 -/***********************************************************************//**  
370 - * @brief Copy class members  
371 - *  
372 - * @param[in] source Point source cube.  
373 - ***************************************************************************/  
374 -void GCTACubeSourcePoint::copy_members(const GCTACubeSourcePoint& source)  
375 -{  
376 - // Copy members  
377 - m_dir = source.m_dir;  
378 - m_aeff = source.m_aeff;  
379 - m_delta_map = source.m_delta_map;  
380 - m_psf = source.m_psf;  
381 - m_deltas = source.m_deltas;  
382 -  
383 - // Return  
384 - return;  
385 -}  
386 -  
387 -/***********************************************************************//**  
388 - * @brief Delete class members  
389 - ***************************************************************************/  
390 -void GCTACubeSourcePoint::free_members(void)  
391 -{  
392 - // Return  
393 - return;  
394 -}