GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMD1Response.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMD1Response.cpp - COMPTEL D1 module response class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-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 GCOMD1Response.cpp
23  * @brief COMPTEL D1 module response class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GMath.hpp"
32 #include "GCOMD1Response.hpp"
33 #include "GFitsTable.hpp"
34 #include "GFitsBinTable.hpp"
35 #include "GFitsTableFloatCol.hpp"
36 #include "GFitsTableDoubleCol.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 
40 /* __ Macros _____________________________________________________________ */
41 
42 /* __ Coding definitions _________________________________________________ */
43 
44 /* __ Debug definitions __________________________________________________ */
45 
46 /* __ Constants __________________________________________________________ */
47 
48 
49 /*==========================================================================
50  = =
51  = Constructors/destructors =
52  = =
53  ==========================================================================*/
54 
55 /***********************************************************************//**
56  * @brief Void constructor
57  *
58  * Creates an empty COMPTEL D1 module response.
59  ***************************************************************************/
61 {
62  // Initialise members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief Copy constructor
72  *
73  * @param[in] rsp COMPTEL D1 module response.
74  **************************************************************************/
76 {
77  // Initialise members
78  init_members();
79 
80  // Copy members
81  copy_members(rsp);
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Response constructor
90  *
91  * @param[in] caldb Calibration database.
92  * @param[in] sdaname SDA response name.
93  *
94  * Create COMPTEL D1 module response by loading an SDA file from a
95  * calibration database.
96  ***************************************************************************/
97 GCOMD1Response::GCOMD1Response(const GCaldb& caldb, const std::string& sdaname)
98 {
99  // Initialise members
100  init_members();
101 
102  // Set calibration database
103  this->caldb(caldb);
104 
105  // Load D1 module response
106  this->load(sdaname);
107 
108  // Return
109  return;
110 }
111 
112 
113 /***********************************************************************//**
114  * @brief Destructor
115  *
116  * Destroys instance of COMPTEL response object.
117  ***************************************************************************/
119 {
120  // Free members
121  free_members();
122 
123  // Return
124  return;
125 }
126 
127 
128 /*==========================================================================
129  = =
130  = Operators =
131  = =
132  ==========================================================================*/
133 
134 /***********************************************************************//**
135  * @brief Assignment operator
136  *
137  * @param[in] rsp COMPTEL D1 module response.
138  * @return COMPTEL D1 module response.
139  ***************************************************************************/
141 {
142  // Execute only if object is not identical
143  if (this != &rsp) {
144 
145  // Free members
146  free_members();
147 
148  // Initialise members
149  init_members();
150 
151  // Copy members
152  copy_members(rsp);
153 
154  } // endif: object was not identical
155 
156  // Return this object
157  return *this;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief D1 module response evaluation operator
163  *
164  * @param[in] etrue True energy (MeV).
165  * @param[in] ereco Reconstructed energy (MeV).
166  * @return COMPTEL D1 module response.
167  *
168  * Computes
169  *
170  * \f[
171  * R_{\rm D1}(E|E_0) = n \times \left[
172  * A_1 \exp \left( -\frac{1}{2} \frac{(E_0-E)^2}{\sigma^2(E_0)} \right)
173  * \right]
174  * \f]
175  *
176  * where
177  * \f$B1\f$ is the amplitude of the photo peak,
178  * \f$\sigma(E_0)\f$ is width of the photo peak, and
179  * \f$E_0\f$ is the position of the photo peak.
180  * The constant \f$n\f$ is chosen so that
181  *
182  * \f[
183  * \int_{E} R_{\rm D1}(E|E_0) dE = 1
184  * \f]
185  ***************************************************************************/
186 double GCOMD1Response::operator()(const double& etrue, const double& ereco) const
187 {
188  // Initialise response with zero
189  double response = 0.0;
190 
191  // Continue only if a response was loaded
192  if (!m_energies.is_empty()) {
193 
194  // Update response evaluation cache
195  update_cache(etrue);
196 
197  // Continue only if amplitude is positive
198  if (m_amplitude > 0.0) {
199 
200  // Compute D1 module response. Here is where the real magic
201  // happens. Only consider the response within 5 sigma.
202  double arg = (m_position-ereco) / m_sigma;
203  if (std::abs(arg) < 5.0) {
204  response = m_amplitude * std::exp(-0.5 * arg * arg);
205  }
206 
207  }
208 
209  // Apply threshold. We use a width of 0.1 MeV for the high-energy
210  // threshold width to assure a smooth transition to zero.
211  if (ereco < m_emin) {
212  double arg = (m_emin-ereco) / m_ewidth;
213  response *= std::exp(-0.5 * arg * arg);
214  }
215  else if (ereco > m_emax) {
216  double arg = (m_emax-ereco) / 0.1;
217  response *= std::exp(-0.5 * arg * arg);
218  }
219 
220  } // endif: response was loaded
221 
222  // Return response
223  return response;
224 }
225 
226 
227 /*==========================================================================
228  = =
229  = Public methods =
230  = =
231  ==========================================================================*/
232 
233 /***********************************************************************//**
234  * @brief Clear instance
235  *
236  * Clears COMPTEL D1 module response object by resetting all members to an
237  * initial state. Any information that was present in the object before will
238  * be lost.
239  ***************************************************************************/
241 {
242  // Free class members
243  free_members();
244 
245  // Initialise members
246  init_members();
247 
248  // Return
249  return;
250 }
251 
252 
253 /***********************************************************************//**
254  * @brief Clone instance
255  *
256  * @return Pointer to deep copy of COMPTEL D1 module response.
257  ***************************************************************************/
259 {
260  return new GCOMD1Response(*this);
261 }
262 
263 
264 /***********************************************************************//**
265  * @brief Load COMPTEL D1 module response.
266  *
267  * @param[in] sdaname COMPTEL D1 module response name.
268  *
269  * Loads the COMPTEL D1 module response with specified name @p sdaname. The
270  * method first searchs for an appropriate response in the calibration
271  * database. If no appropriate response is found, the method takes the
272  * database root path and response name to build the full path to the
273  * response file, and tries to load the response from these paths.
274  ***************************************************************************/
275 void GCOMD1Response::load(const std::string& sdaname)
276 {
277  // Clear instance but conserve calibration database
278  GCaldb caldb = m_caldb;
279  clear();
280  m_caldb = caldb;
281 
282  // First attempt reading the response using the GCaldb interface
283  GFilename filename = m_caldb.filename("","","SDA","","",sdaname);
284 
285  // If filename is empty then build filename from CALDB root path and
286  // response name
287  if (filename.is_empty()) {
288  filename = gammalib::filepath(m_caldb.rootdir(), sdaname);
289  if (!filename.exists()) {
290  GFilename testname = filename + ".fits";
291  if (testname.exists()) {
292  filename = testname;
293  }
294  }
295  }
296 
297  // Open FITS file
298  GFits fits(filename);
299 
300  // Get SDA table
301  const GFitsTable& sda = *fits.table(1);
302 
303  // Read SDA
304  read(sda);
305 
306  // Close SDA FITS file
307  fits.close();
308 
309  // Return
310  return;
311 }
312 
313 
314 /***********************************************************************//**
315  * @brief Read COMPTEL D1 module response.
316  *
317  * @param[in] table FITS table.
318  *
319  * Read the COMPTEL D1 module response from a SDA FITS table.
320  ***************************************************************************/
322 {
323  // Initialise COMPTEL D1 module response vectors
324  m_energies.clear();
325  m_positions.clear();
326  m_sigmas.clear();
327  m_amplitudes.clear();
328  m_emins.clear();
329  m_ewidths.clear();
330  m_emaxs.clear();
331 
332  // Extract number of entries in table
333  int num = table.nrows();
334 
335  // If there are entries then read them
336  if (num > 0) {
337 
338  // Get column pointers
339  const GFitsTableCol* ptr_energy = table["ENERGY"];
340  const GFitsTableCol* ptr_position = table["POSITION"];
341  const GFitsTableCol* ptr_sigma = table["WIDTH"];
342  const GFitsTableCol* ptr_amplitude = table["AMPLITUDE"];
343  const GFitsTableCol* ptr_emin = table["EMIN"];
344  const GFitsTableCol* ptr_ewidth = table["EWIDTH"];
345  const GFitsTableCol* ptr_emax = table["EMAX"];
346 
347  // Copy data from table into vectors
348  for (int i = 0; i < num; ++i) {
349  m_energies.append(ptr_energy->real(i));
350  m_positions.push_back(ptr_position->real(i));
351  m_sigmas.push_back(ptr_sigma->real(i));
352  m_amplitudes.push_back(ptr_amplitude->real(i));
353  m_emins.push_back(ptr_emin->real(i));
354  m_ewidths.push_back(ptr_ewidth->real(i));
355  m_emaxs.push_back(ptr_emax->real(i));
356  }
357 
358  } // endif: there were entries
359 
360  // Return
361  return;
362 }
363 
364 
365 /***********************************************************************//**
366  * @brief Write COMPTEL D1 module response.
367  *
368  * @param[in] table FITS table.
369  *
370  * Write the COMPTEL D1 module response into a SDA FITS table.
371  ***************************************************************************/
373 {
374  // Set extension name
375  table.extname("SDA");
376 
377  // Set number of entries
378  int num = m_energies.size();
379 
380  // Allocate columns
381  GFitsTableFloatCol col_energy("ENERGY", num);
382  GFitsTableDoubleCol col_position("POSITION", num);
383  GFitsTableDoubleCol col_width("WIDTH", num);
384  GFitsTableDoubleCol col_amplitude("AMPLITUDE", num);
385  GFitsTableDoubleCol col_emin("EMIN", num);
386  GFitsTableDoubleCol col_ewidth("EWIDTH", num);
387  GFitsTableDoubleCol col_emax("EMAX", num);
388 
389  // Set column units
390  col_energy.unit("MeV");
391  col_position.unit("MeV");
392  col_width.unit("MeV");
393  col_emin.unit("MeV");
394  col_ewidth.unit("MeV");
395  col_emax.unit("MeV");
396 
397  // Copy data from vectors into table
398  for (int i = 0; i < num; ++i) {
399  col_energy(i) = m_energies[i];
400  col_position(i) = m_positions[i];
401  col_width(i) = m_sigmas[i];
402  col_amplitude(i) = m_amplitudes[i];
403  col_emin(i) = m_emins[i];
404  col_ewidth(i) = m_ewidths[i];
405  col_emax(i) = m_emaxs[i];
406  }
407 
408  // Append columns to table
409  table.append(col_energy);
410  table.append(col_position);
411  table.append(col_width);
412  table.append(col_amplitude);
413  table.append(col_emin);
414  table.append(col_ewidth);
415  table.append(col_emax);
416 
417  // Return
418  return;
419 }
420 
421 
422 /***********************************************************************//**
423  * @brief Print COMPTEL D1 module response information
424  *
425  * @param[in] chatter Chattiness.
426  * @return String containing COMPTEL D1 module response information.
427  ***************************************************************************/
428 std::string GCOMD1Response::print(const GChatter& chatter) const
429 {
430  // Initialise result string
431  std::string result;
432 
433  // Continue only if chatter is not silent
434  if (chatter != SILENT) {
435 
436  // Append header
437  result.append("=== GCOMD1Response ===");
438 
439  // Append D1 module response information
440  result.append("\n"+gammalib::parformat("Energy range"));
441  if (m_energies.size() > 0) {
442  result.append(gammalib::str(m_energies[0])+ " - ");
443  result.append(gammalib::str(m_energies[m_energies.size()-1])+ " MeV");
444  }
445  else {
446  result.append("not defined");
447  }
448  result.append("\n"+gammalib::parformat("Entries"));
449  result.append(gammalib::str(m_energies.size()));
450 
451  // Append calibration database
452  result.append("\n"+m_caldb.print(chatter));
453 
454  // Append information
455 
456  } // endif: chatter was not silent
457 
458  // Return result
459  return result;
460 }
461 
462 
463 /*==========================================================================
464  = =
465  = Private methods =
466  = =
467  ==========================================================================*/
468 
469 /***********************************************************************//**
470  * @brief Initialise class members
471  ***************************************************************************/
473 {
474  // Initialise members
475  m_caldb.clear();
476  m_energies.clear();
477  m_positions.clear();
478  m_sigmas.clear();
479  m_amplitudes.clear();
480  m_emins.clear();
481  m_ewidths.clear();
482  m_emaxs.clear();
483 
484  // Initialise pre-computation cache
485  m_energy = 0.0;
486  m_position = 0.0;
487  m_sigma = 0.0;
488  m_amplitude = 0.0;
489  m_emin = 0.0;
490  m_ewidth = 0.0;
491  m_emax = 0.0;
492 
493  // Return
494  return;
495 }
496 
497 
498 /***********************************************************************//**
499  * @brief Copy class members
500  *
501  * @param[in] rsp COMPTEL response.
502  ***************************************************************************/
504 {
505  // Copy attributes
506  m_caldb = rsp.m_caldb;
507  m_energies = rsp.m_energies;
508  m_positions = rsp.m_positions;
509  m_sigmas = rsp.m_sigmas;
511  m_emins = rsp.m_emins;
512  m_ewidths = rsp.m_ewidths;
513  m_emaxs = rsp.m_emaxs;
514 
515  // Copy pre-computation cache
516  m_energy = rsp.m_energy;
517  m_position = rsp.m_position;
518  m_sigma = rsp.m_sigma;
519  m_amplitude = rsp.m_amplitude;
520  m_emin = rsp.m_emin;
521  m_ewidth = rsp.m_ewidth;
522  m_emax = rsp.m_emax;
523 
524  // Return
525  return;
526 }
527 
528 
529 /***********************************************************************//**
530  * @brief Delete class members
531  ***************************************************************************/
533 {
534  // Return
535  return;
536 }
537 
538 
539 /***********************************************************************//**
540  * @brief Update computation cache
541  *
542  * @param[in] etrue True energy (MeV).
543  *
544  * The method assumes that there is a valid D1 module response.
545  ***************************************************************************/
546 void GCOMD1Response::update_cache(const double& etrue) const
547 {
548  // Update only if the true energy has changed
549  if (etrue != m_energy) {
550 
551  // Set true energy
552  m_energy = etrue;
553 
554  // If true energy is below lowest energy or above largest energy
555  // then set response to zero
556  if ((etrue < m_energies[0]) ||
557  (etrue > m_energies[m_energies.size()-1])) {
558  m_position = 0.0;
559  m_sigma = 0.0;
560  m_amplitude = 0.0;
561  m_emin = 0.0;
562  m_ewidth = 0.0;
563  m_emax = 0.0;
564  }
565 
566  // ... otherwise interpolate response parameters
567  else {
568 
569  // Interpolate response parameters
576 
577  // Re-compute Gaussian amplitude to assure normalization
579 
580  }
581 
582  } // endif: true energy has changed
583 
584  // Return
585  return;
586 }
void read(const GFitsTable &table)
Read COMPTEL D1 module response.
const double sqrt_twopi
Definition: GMath.hpp:54
void unit(const std::string &unit)
Set column unit.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
void load(const std::string &sdaname)
Load COMPTEL D1 module response.
~GCOMD1Response(void)
Destructor.
FITS table double column class interface definition.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
std::vector< double > m_emaxs
Upper energy limit of D1.
const GCaldb & caldb(void) const
Return calibration database.
GCOMD1Response * clone(void) const
Clone instance.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
GFilename filename(const std::string &detector, const std::string &filter, const std::string &codename, const std::string &date, const std::string &time, const std::string &expr)
Return calibration file name based on selection parameters.
Definition: GCaldb.cpp:507
void write(GFitsBinTable &table)
Write COMPTEL D1 module response.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
FITS table float column class interface definition.
FITS file class.
Definition: GFits.hpp:63
GCaldb m_caldb
Calibration database.
bool is_empty(void) const
Signals if there are no nodes in node array.
Definition: GNodeArray.hpp:202
GNodeArray m_energies
Input energies.
std::vector< double > m_ewidths
Lower energy threshold width of D1.
void free_members(void)
Delete class members.
void init_members(void)
Initialise class members.
Calibration database class.
Definition: GCaldb.hpp:66
void update_cache(const double &etrue) const
Update computation cache.
Interface for the COMPTEL D1 module response class.
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
std::vector< double > m_sigmas
Photo peak width in MeV.
GChatter
Definition: GTypemaps.hpp:33
GCOMD1Response(void)
Void constructor.
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL D1 module response information.
std::vector< double > m_amplitudes
Photo peak amplitude.
COMPTEL D1 module response class interface definition.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void copy_members(const GCOMD1Response &rsp)
Copy class members.
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
virtual double real(const int &row, const int &inx=0) const =0
std::vector< double > m_positions
Photo peak position in MeV.
FITS binary table class.
GCOMD1Response & operator=(const GCOMD1Response &rsp)
Assignment operator.
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
FITS binary table class definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
FITS table float column.
void clear(void)
Clear instance.
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:386
double operator()(const double &etrue, const double &ereco) const
D1 module response evaluation operator.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
std::vector< double > m_emins
Lower energy threshold of D1.
void clear(void)
Clear calibration database.
Definition: GCaldb.cpp:194
FITS table double column.
std::string rootdir(void) const
Return path to CALDB root directory.
Definition: GCaldb.cpp:259
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.