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