GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMSupport.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMSupport.cpp - COMPTEL support functions *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2021 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 GCOMSupport.hpp
23  * @brief Implementation of support function used by COMPTEL classes
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 "GWcs.hpp"
33 #include "GWcsCAR.hpp"
34 #include "GSkyMap.hpp"
35 #include "GCOMSupport.hpp"
36 
37 /* __ Coding definitions _________________________________________________ */
38 
39 /* __ Debug definitions __________________________________________________ */
40 
41 
42 /***********************************************************************//**
43  * @brief Changes Mercator's projection to cartesian projection
44  *
45  * @param[in,out] map Skymap.
46  *
47  * Changes the World Coordinate System of a sky map from Mercartor's
48  * projection to a cartesian projection. This transformation is needed to
49  * correct the wrong WCS headers that are found in the COMPTEL files
50  * distributed through the HEASARC web site.
51  *
52  * The method operates only on sky maps that use Mercator's projection.
53  * Nothing is done for all other sky map. This method can thus be
54  * transparently be applied to the sky maps read from COMPTEL data files.
55  ***************************************************************************/
57 {
58  // Get WCS poiunter
59  const GWcs* wcs = dynamic_cast<const GWcs*>(map.projection());
60 
61  // Apply kluge only if the skymap is in Mercator's projection
62  if (wcs != NULL && wcs->code() == "MER") {
63 
64  // Extract original WCS definition
65  double crval1 = wcs->crval(0);
66  double crval2 = wcs->crval(1);
67  double crpix1 = wcs->crpix(0);
68  double crpix2 = wcs->crpix(1);
69  double cdelt1 = wcs->cdelt(0);
70  double cdelt2 = wcs->cdelt(1);
71  int naxis1 = map.nx();
72  int naxis2 = map.ny();
73 
74  // Compute adjustment of WCS definition to map centre
75  double adjust1 = (double(naxis1) + 1.0) / 2.0 - crpix1;
76  double adjust2 = (double(naxis2) + 1.0) / 2.0 - crpix2;
77 
78  // Adjust WCS definition to map centre
79  crval1 += adjust1;
80  crpix1 += adjust1;
81  crval2 += adjust2;
82  crpix2 += adjust2;
83 
84  // Put latitude to zero to assure correct projection
85  crpix2 -= crval2;
86  crval2 = 0.0;
87 
88  // Allocate WCS
89  GWcsCAR car(wcs->coordsys(), crval1, crval2,
90  crpix1, crpix2,
91  cdelt1, cdelt2);
92 
93  // Set new WCS
94  map.projection(car);
95 
96  } // endif: skymap was in Mercator projection
97 
98  // Return
99  return;
100 }
101 
102 
103 /***********************************************************************//**
104  * @brief Return D1 energy deposit
105  *
106  * @param[in] energy Input photon energy (MeV).
107  * @param[in] phigeo Geometrical scatter angle (deg).
108  * @return D1 energy deposit in MeV.
109  ***************************************************************************/
110 double gammalib::com_energy1(const double& energy, const double& phigeo)
111 {
112  // Compute D2 energy deposit
113  double e2 = com_energy2(energy, phigeo);
114 
115  // Return D1 energy deposit
116  return (energy - e2);
117 }
118 
119 
120 /***********************************************************************//**
121  * @brief Return D2 energy deposit
122  *
123  * @param[in] energy Input photon energy (MeV).
124  * @param[in] phigeo Geometrical scatter angle (deg).
125  * @return D2 energy deposit in MeV.
126  ***************************************************************************/
127 double gammalib::com_energy2(const double& energy, const double& phigeo)
128 {
129  // Compute 1-cos(phigeo)
130  double one_minus_cos = 1.0 - std::cos(phigeo * gammalib::deg2rad);
131 
132  // Compute D2 energy deposit
133  double e2 = energy / (one_minus_cos * energy / gammalib::mec2 + 1.0);
134 
135  // Return D2 energy deposit
136  return e2;
137 }
138 
139 
140 /***********************************************************************//**
141  * @brief Return D2 module exclusion region X position
142  *
143  * @param[in] id2 Module identifier [1,...,14].
144  * @return Exclusion region X position (cm).
145  *
146  * Returns the D2 module exclusion region X position for a given module.
147  * The method does not check the validity of the module identifier.
148  *
149  * The values have been implemented from the MPE-FPM-4 file.
150  ***************************************************************************/
151 const double& gammalib::com_exd2x(const int& id2)
152 {
153  // Set D2 module exclusion region X positions
154  static const double exd2x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
155  0.0, 0.0, 0.0, -51.7, 0.0, 9.0, -34.7};
156 
157  // Return value
158  return exd2x[id2-1];
159 }
160 
161 
162 /***********************************************************************//**
163  * @brief Return D2 module exclusion region Y position
164  *
165  * @param[in] id2 Module identifier [1,...,14].
166  * @return Exclusion region Y position (cm).
167  *
168  * Returns the D2 module exclusion region Y position for a given module.
169  * The method does not check the validity of the module identifier.
170  *
171  * The values have been implemented from the MPE-FPM-4 file.
172  ***************************************************************************/
173 const double& gammalib::com_exd2y(const int& id2)
174 {
175  // Set D2 module exclusion region Y positions
176  static const double exd2y[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
177  0.0, 0.0, 0.0, 8.7, 0.0, +41.2, +49.0};
178 
179  // Return value
180  return exd2y[id2-1];
181 }
182 
183 
184 /***********************************************************************//**
185  * @brief Return D2 module exclusion region radius
186  *
187  * @param[in] id2 Module identifier [1,...,14].
188  * @return Exclusion region radius (cm).
189  *
190  * Returns the D2 module exclusion region radius for a given module.
191  * The method does not check the validity of the module identifier.
192  *
193  * The values have been implemented from the MPE-FPM-4 file.
194  ***************************************************************************/
195 const double& gammalib::com_exd2r(const int& id2)
196 {
197  // Set D2 module exclusion region radii
198  static const double exd2r[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
199  0.0, 0.0, 0.0, 9.0, 0.0, 9.0, 9.0};
200 
201  // Return value
202  return exd2r[id2-1];
203 }
const double mec2
Definition: GTools.hpp:52
Sky map class.
Definition: GSkyMap.hpp:89
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:379
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:465
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
double com_energy2(const double &energy, const double &phigeo)
Return D2 energy deposit.
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:887
Gammalib tools definition.
Sky map class definition.
const double & com_exd2r(const int &id2)
Return D2 module exclusion region radius.
virtual std::string code(void) const =0
Implementation of support function used by COMPTEL classes.
const double deg2rad
Definition: GMath.hpp:43
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:911
const double & com_exd2y(const int &id2)
Return D2 module exclusion region Y position.
void com_wcs_mer2car(GSkyMap &map)
Changes Mercator&#39;s projection to cartesian projection.
Definition: GCOMSupport.cpp:56
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:935
const double & com_exd2x(const int &id2)
Return D2 module exclusion region X position.
Plate carree (CAR) projection class definition.
virtual std::string coordsys(void) const
Returns coordinate system.
Plate carree (CAR) projection class definition.
Definition: GWcsCAR.hpp:43
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:363
double com_energy1(const double &energy, const double &phigeo)
Return D1 energy deposit.
Abstract world coordinate system base class definition.