GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMD2Response.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMD2Response.cpp - COMPTEL D2 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 GCOMD2Response.cpp
23  * @brief COMPTEL D2 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 "GIntegral.hpp"
33 #include "GCOMD2Response.hpp"
34 #include "GFitsTable.hpp"
35 #include "GFitsBinTable.hpp"
36 #include "GFitsTableFloatCol.hpp"
37 #include "GFitsTableDoubleCol.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 #define G_UPDATE_RESPONSE_VECTOR_NO_WARNINGS
45 
46 /* __ Debug definitions __________________________________________________ */
47 //#define G_DEBUG_UPDATE_RESPONSE_VECTOR
48 
49 /* __ Constants __________________________________________________________ */
50 
51 
52 /*==========================================================================
53  = =
54  = Constructors/destructors =
55  = =
56  ==========================================================================*/
57 
58 /***********************************************************************//**
59  * @brief Void constructor
60  *
61  * Creates an empty COMPTEL D2 module response.
62  ***************************************************************************/
64 {
65  // Initialise members
66  init_members();
67 
68  // Return
69  return;
70 }
71 
72 
73 /***********************************************************************//**
74  * @brief Copy constructor
75  *
76  * @param[in] rsp COMPTEL D2 module response.
77  **************************************************************************/
79 {
80  // Initialise members
81  init_members();
82 
83  // Copy members
84  copy_members(rsp);
85 
86  // Return
87  return;
88 }
89 
90 
91 /***********************************************************************//**
92  * @brief Response constructor
93  *
94  * @param[in] caldb Calibration database.
95  * @param[in] sdbname SDA response name.
96  *
97  * Create COMPTEL D2 module response by loading an SDB file from a
98  * calibration database.
99  ***************************************************************************/
100 GCOMD2Response::GCOMD2Response(const GCaldb& caldb, const std::string& sdbname)
101 {
102  // Initialise members
103  init_members();
104 
105  // Set calibration database
106  this->caldb(caldb);
107 
108  // Load D2 module response
109  this->load(sdbname);
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief Destructor
118  *
119  * Destroys instance of COMPTEL response object.
120  ***************************************************************************/
122 {
123  // Free members
124  free_members();
125 
126  // Return
127  return;
128 }
129 
130 
131 /*==========================================================================
132  = =
133  = Operators =
134  = =
135  ==========================================================================*/
136 
137 /***********************************************************************//**
138  * @brief Assignment operator
139  *
140  * @param[in] rsp COMPTEL D2 module response.
141  * @return COMPTEL D2 module response.
142  ***************************************************************************/
144 {
145  // Execute only if object is not identical
146  if (this != &rsp) {
147 
148  // Free members
149  free_members();
150 
151  // Initialise members
152  init_members();
153 
154  // Copy members
155  copy_members(rsp);
156 
157  } // endif: object was not identical
158 
159  // Return this object
160  return *this;
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief D2 module response evaluation operator
166  *
167  * @param[in] etrue True energy (MeV).
168  * @param[in] ereco Reconstructed energy (MeV).
169  * @return COMPTEL D2 module response.
170  *
171  * @todo Add Compton components to response
172  ***************************************************************************/
173 double GCOMD2Response::operator()(const double& etrue, const double& ereco) const
174 {
175  // Initialise response and area with zero
176  double response = 0.0;
177 
178  // Continue only if a response was loaded
179  if (!m_energies.is_empty()) {
180 
181  // Update response vector
182  update_response_vector(etrue);
183 
184  // Get response value from response vector
185  response = m_rsp_energies.interpolate(ereco, m_rsp_values);
186  if (response < 0.0) {
187  response = 0.0;
188  }
189 
190  // Apply threshold. We use a width of 0.1 MeV for the high-energy
191  // threshold width to assure a smooth transition to zero.
192  if (ereco < m_emin) {
193  double arg = (m_emin-ereco) / m_ewidth;
194  response *= std::exp(-0.5 * arg * arg);
195  }
196  else if (ereco > m_emax) {
197  double arg = (m_emax-ereco) / 0.1;
198  response *= std::exp(-0.5 * arg * arg);
199  }
200 
201  } // endif: response was loaded
202 
203  // Return response
204  return response;
205 }
206 
207 
208 /*==========================================================================
209  = =
210  = Public methods =
211  = =
212  ==========================================================================*/
213 
214 /***********************************************************************//**
215  * @brief Clear instance
216  *
217  * Clears COMPTEL D2 module response object by resetting all members to an
218  * initial state. Any information that was present in the object before will
219  * be lost.
220  ***************************************************************************/
222 {
223  // Free class members
224  free_members();
225 
226  // Initialise members
227  init_members();
228 
229  // Return
230  return;
231 }
232 
233 
234 /***********************************************************************//**
235  * @brief Clone instance
236  *
237  * @return Pointer to deep copy of COMPTEL D2 module response.
238  ***************************************************************************/
240 {
241  return new GCOMD2Response(*this);
242 }
243 
244 
245 /***********************************************************************//**
246  * @brief Load COMPTEL D2 module response.
247  *
248  * @param[in] sdbname COMPTEL D2 module response name.
249  *
250  * Loads the COMPTEL D2 module response with specified name @p sdbname. The
251  * method first searchs for an appropriate response in the calibration
252  * database. If no appropriate response is found, the method takes the
253  * database root path and response name to build the full path to the
254  * response file, and tries to load the response from these paths.
255  ***************************************************************************/
256 void GCOMD2Response::load(const std::string& sdbname)
257 {
258  // Clear instance but conserve calibration database
259  GCaldb caldb = m_caldb;
260  clear();
261  m_caldb = caldb;
262 
263  // First attempt reading the response using the GCaldb interface
264  GFilename filename = m_caldb.filename("","","SDB","","",sdbname);
265 
266  // If filename is empty then build filename from CALDB root path and
267  // response name
268  if (filename.is_empty()) {
269  filename = gammalib::filepath(m_caldb.rootdir(), sdbname);
270  if (!filename.exists()) {
271  GFilename testname = filename + ".fits";
272  if (testname.exists()) {
273  filename = testname;
274  }
275  }
276  }
277 
278  // Open FITS file
279  GFits fits(filename);
280 
281  // Get SDB table
282  const GFitsTable& sdb = *fits.table(1);
283 
284  // Read SDB
285  read(sdb);
286 
287  // Close SDB FITS file
288  fits.close();
289 
290  // Return
291  return;
292 }
293 
294 
295 /***********************************************************************//**
296  * @brief Read COMPTEL D2 module response.
297  *
298  * @param[in] table FITS table.
299  *
300  * Read the COMPTEL D2 module response from a SDB FITS table.
301  ***************************************************************************/
303 {
304  // Initialise COMPTEL D2 module response vectors
305  m_energies.clear();
306  m_positions.clear();
307  m_sigmas.clear();
308  m_amplitudes.clear();
309  m_escapes1.clear();
310  m_escapes2.clear();
311  m_tails.clear();
312  m_backgrounds.clear();
313  m_emins.clear();
314  m_ewidths.clear();
315  m_emaxs.clear();
316 
317  // Extract number of entries in table
318  int num = table.nrows();
319 
320  // If there are entries then read them
321  if (num > 0) {
322 
323  // Get column pointers
324  const GFitsTableCol* ptr_energy = table["ENERGY"];
325  const GFitsTableCol* ptr_position = table["POSITION"];
326  const GFitsTableCol* ptr_sigma = table["WIDTH"];
327  const GFitsTableCol* ptr_amplitude = table["AMPLITUDE"];
328  const GFitsTableCol* ptr_escape1 = table["ESCAPE1"];
329  const GFitsTableCol* ptr_escape2 = table["ESCAPE2"];
330  const GFitsTableCol* ptr_tail = table["TAIL"];
331  const GFitsTableCol* ptr_background = table["BACKGROUND"];
332  const GFitsTableCol* ptr_emin = table["EMIN"];
333  const GFitsTableCol* ptr_ewidth = table["EWIDTH"];
334  const GFitsTableCol* ptr_emax = table["EMAX"];
335 
336  // Copy data from table into vectors
337  for (int i = 0; i < num; ++i) {
338  m_energies.append(ptr_energy->real(i));
339  m_positions.push_back(ptr_position->real(i));
340  m_sigmas.push_back(ptr_sigma->real(i));
341  m_amplitudes.push_back(ptr_amplitude->real(i));
342  m_escapes1.push_back(ptr_escape1->real(i));
343  m_escapes2.push_back(ptr_escape2->real(i));
344  m_tails.push_back(ptr_tail->real(i));
345  m_backgrounds.push_back(ptr_background->real(i));
346  m_emins.push_back(ptr_emin->real(i));
347  m_ewidths.push_back(ptr_ewidth->real(i));
348  m_emaxs.push_back(ptr_emax->real(i));
349  }
350 
351  } // endif: there were entries
352 
353  // Return
354  return;
355 }
356 
357 
358 /***********************************************************************//**
359  * @brief Write COMPTEL D2 module response.
360  *
361  * @param[in] table FITS table.
362  *
363  * Write the COMPTEL D2 module response into a SDB FITS table.
364  ***************************************************************************/
366 {
367  // Set extension name
368  table.extname("SDB");
369 
370  // Set number of entries
371  int num = m_energies.size();
372 
373  // Allocate columns
374  GFitsTableFloatCol col_energy("ENERGY", num);
375  GFitsTableDoubleCol col_position("POSITION", num);
376  GFitsTableDoubleCol col_width("WIDTH", num);
377  GFitsTableDoubleCol col_amplitude("AMPLITUDE", num);
378  GFitsTableDoubleCol col_escape1("ESCAPE1", num);
379  GFitsTableDoubleCol col_escape2("ESCAPE2", num);
380  GFitsTableDoubleCol col_tail("TAIL", num);
381  GFitsTableDoubleCol col_background("BACKGROUND", num);
382  GFitsTableDoubleCol col_emin("EMIN", num);
383  GFitsTableDoubleCol col_ewidth("EWIDTH", num);
384  GFitsTableDoubleCol col_emax("EMAX", num);
385 
386  // Set column units
387  col_energy.unit("MeV");
388  col_position.unit("MeV");
389  col_width.unit("MeV");
390  col_emin.unit("MeV");
391  col_ewidth.unit("MeV");
392  col_emax.unit("MeV");
393 
394  // Copy data from vectors into table
395  for (int i = 0; i < num; ++i) {
396  col_energy(i) = m_energies[i];
397  col_position(i) = m_positions[i];
398  col_width(i) = m_sigmas[i];
399  col_amplitude(i) = m_amplitudes[i];
400  col_escape1(i) = m_escapes1[i];
401  col_escape2(i) = m_escapes2[i];
402  col_tail(i) = m_tails[i];
403  col_background(i) = m_backgrounds[i];
404  col_emin(i) = m_emins[i];
405  col_ewidth(i) = m_ewidths[i];
406  col_emax(i) = m_emaxs[i];
407  }
408 
409  // Append columns to table
410  table.append(col_energy);
411  table.append(col_position);
412  table.append(col_width);
413  table.append(col_amplitude);
414  table.append(col_escape1);
415  table.append(col_escape2);
416  table.append(col_tail);
417  table.append(col_background);
418  table.append(col_emin);
419  table.append(col_ewidth);
420  table.append(col_emax);
421 
422  // Return
423  return;
424 }
425 
426 
427 /***********************************************************************//**
428  * @brief Print COMPTEL D2 module response information
429  *
430  * @param[in] chatter Chattiness.
431  * @return String containing COMPTEL D2 module response information.
432  ***************************************************************************/
433 std::string GCOMD2Response::print(const GChatter& chatter) const
434 {
435  // Initialise result string
436  std::string result;
437 
438  // Continue only if chatter is not silent
439  if (chatter != SILENT) {
440 
441  // Append header
442  result.append("=== GCOMD2Response ===");
443 
444  // Append D2 module response information
445  result.append("\n"+gammalib::parformat("Energy range"));
446  if (m_energies.size() > 0) {
447  result.append(gammalib::str(m_energies[0])+ " - ");
448  result.append(gammalib::str(m_energies[m_energies.size()-1])+ " MeV");
449  }
450  else {
451  result.append("not defined");
452  }
453  result.append("\n"+gammalib::parformat("Entries"));
454  result.append(gammalib::str(m_energies.size()));
455 
456  // Append calibration database
457  result.append("\n"+m_caldb.print(chatter));
458 
459  // Append information
460 
461  } // endif: chatter was not silent
462 
463  // Return result
464  return result;
465 }
466 
467 
468 /*==========================================================================
469  = =
470  = Private methods =
471  = =
472  ==========================================================================*/
473 
474 /***********************************************************************//**
475  * @brief Initialise class members
476  ***************************************************************************/
478 {
479  // Initialise members
480  m_caldb.clear();
481  m_energies.clear();
482  m_positions.clear();
483  m_sigmas.clear();
484  m_amplitudes.clear();
485  m_escapes1.clear();
486  m_escapes2.clear();
487  m_tails.clear();
488  m_backgrounds.clear();
489  m_emins.clear();
490  m_ewidths.clear();
491  m_emaxs.clear();
492 
493  // Initialise pre-computation cache
494  m_energy = 0.0;
495  m_position = 0.0;
496  m_sigma = 0.0;
497  m_amplitude = 0.0;
498  m_escape1 = 0.0;
499  m_escape2 = 0.0;
500  m_tail = 0.0;
501  m_background = 0.0;
502  m_emin = 0.0;
503  m_ewidth = 0.0;
504  m_emax = 0.0;
505  m_pos_escape1 = 0.0;
506  m_pos_escape2 = 0.0;
507  m_wgt_photo = 0.0;
508  m_wgt_escape1 = 0.0;
509  m_wgt_escape2 = 0.0;
510  m_compton_edge = 0.0;
511 
512  // Pre-computed response vector
513  m_rsp_etrue = 0.0;
515  m_rsp_values.clear();
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_escapes1 = rsp.m_escapes1;
536  m_escapes2 = rsp.m_escapes2;
537  m_tails = rsp.m_tails;
539  m_emins = rsp.m_emins;
540  m_ewidths = rsp.m_ewidths;
541  m_emaxs = rsp.m_emaxs;
542 
543  // Copy pre-computation cache
544  m_energy = rsp.m_energy;
545  m_position = rsp.m_position;
546  m_sigma = rsp.m_sigma;
547  m_amplitude = rsp.m_amplitude;
548  m_escape1 = rsp.m_escape1;
549  m_escape2 = rsp.m_escape2;
550  m_tail = rsp.m_tail;
552  m_emin = rsp.m_emin;
553  m_ewidth = rsp.m_ewidth;
554  m_emax = rsp.m_emax;
557  m_wgt_photo = rsp.m_wgt_photo;
561 
562  // Copy pre-computed response vector
563  m_rsp_etrue = rsp.m_rsp_etrue;
566 
567  // Return
568  return;
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Delete class members
574  ***************************************************************************/
576 {
577  // Return
578  return;
579 }
580 
581 
582 /***********************************************************************//**
583  * @brief Update computation cache
584  *
585  * @param[in] etrue True energy (MeV).
586  *
587  * The method assumes that there is a valid D2 module response.
588  ***************************************************************************/
589 void GCOMD2Response::update_cache(const double& etrue) const
590 {
591  // Update only if the true energy has changed
592  if (etrue != m_energy) {
593 
594  // Set true energy
595  m_energy = etrue;
596 
597  // If true energy is below lowest energy or above largest energy
598  // then set response to zero
599  if ((etrue < m_energies[0]) ||
600  (etrue > m_energies[m_energies.size()-1])) {
601  m_position = 0.0;
602  m_sigma = 0.0;
603  m_amplitude = 0.0;
604  m_escape1 = 0.0;
605  m_escape2 = 0.0;
606  m_tail = 0.0;
607  m_background = 0.0;
608  m_emin = 0.0;
609  m_ewidth = 0.0;
610  m_emax = 0.0;
611  m_pos_escape1 = 0.0;
612  m_pos_escape2 = 0.0;
613  m_wgt_escape1 = 0.0;
614  m_wgt_escape2 = 0.0;
615  m_compton_edge = 0.0;
616  }
617 
618  // ... otherwise interpolate response parameters
619  else {
620 
621  // Interpolate response parameters
632 
633  // Derive escape peak positions
636 
637  // Derive inverse standard deviations of photo and escape peaks
639  if (m_wgt_photo > 0.0) {
640  m_wgt_photo = 1.0 / m_wgt_photo;
641  }
643  if (m_wgt_escape1 > 0.0) {
645  }
647  if (m_wgt_escape2 > 0.0) {
649  }
650 
651  // Derive Compton edge
652  m_compton_edge = m_position / (1.0 + 0.5 * gammalib::mec2 / m_position);
653 
654  } // endif: true energy is in valid range
655 
656  } // endif: true energy has changed
657 
658  // Return
659  return;
660 }
661 
662 
663 /***********************************************************************//**
664  * @brief Update response vector
665  *
666  * @param[in] etrue True energy (MeV).
667  *
668  * Updates the vector that stores the normalised D2 module response as
669  * function of energy. The vector is computed using
670  *
671  * \f[
672  * R_{\rm D2}(E|E_0) = n \times \left[
673  * B_1 \exp \left( -\frac{1}{2} \frac{(E_0-E)^2}{\sigma^2(E_0)} \right) +
674  * B_2 \exp \left( -\frac{1}{2} \frac{(E_0-m_e c^2-E)^2}{\sigma^2(E_0-m_e c^2)} \right) +
675  * B_3 \exp \left( -\frac{1}{2} \frac{(E_0-2m_e c^2-E)^2}{\sigma^2(E_0-2m_e c^2)} \right) +
676  * B_4 \int_{E'} KN_{\rm mod}(E'|E,E_0) dE' +
677  * B_5 \int_{E'} B_{\rm c}(E'|E,E_0) dE'
678  * \right]
679  * \f]
680  *
681  * where
682  * \f$B1\f$ is the amplitude of the photo peak,
683  * \f$B2\f$ is the amplitude of the first escape peak,
684  * \f$B3\f$ is the amplitude of the second escape peak,
685  * \f$B4\f$ is the amplitude of the Compton tail,
686  * \f$B5\f$ is the amplitude of the Compton background,
687  * \f$E_0\f$ is the position of the photo peak, and
688  * \f$\sigma(E)\f$ is the energy dependent width of the photo peak. The
689  * constant \f$n\f$ is chosen so that
690  *
691  * \f[
692  * \int_{E} R_{\rm D2}(E|E_0) dE = 1
693  * \f]
694  *
695  * The method assumes that there is a valid D2 module response.
696  ***************************************************************************/
697 void GCOMD2Response::update_response_vector(const double& etrue) const
698 {
699  // Update only if the true energy has changed
700  if (etrue != m_rsp_etrue) {
701 
702  // Set true energy
703  m_rsp_etrue = etrue;
704 
705  // Update cache to determine the spectral parameters at the true
706  // energy
707  update_cache(etrue);
708 
709  // Continue only if position is positive
710  if (m_position > 0.0) {
711 
712  // Clear response vectors (this is very very important, otherwise
713  // some old elements reside and we just append to them, and things
714  // get out of order!!!)
716  m_rsp_values.clear();
717 
718  // Initialise response vector
719  double ebin = 0.001 * m_position;
720  int nbin = int(1.35 * m_position / ebin);
721  if (nbin < 1) {
722  nbin = 1;
723  }
724  m_rsp_energies.reserve(nbin);
725  m_rsp_values.reserve(nbin);
726 
727  // Initialise response vector normalisation
728  double norm = 0.0;
729 
730  // Fill response vector
731  for (int i = 0; i < nbin; ++i) {
732 
733  // Compute bin energy
734  double ereco = (i + 0.5) * ebin;
735 
736  // Store bin energy
737  m_rsp_energies.append(ereco);
738 
739  // Initialise response value
740  double value = 0.0;
741 
742  // Add photo peak if amplitude is positive and reconstructed
743  // energy is within 5 sigma of the Gaussian position
744  if (m_amplitude > 0.0) {
745  double arg = (m_position-ereco) * m_wgt_photo;
746  if (std::abs(arg) < 5.0) {
747  value += m_amplitude * std::exp(-0.5 * arg * arg);
748  }
749  }
750 
751  // Add first escape peak if amplitude is positive and
752  // reconstructed energy is within 5 sigma of the Gaussian
753  // position
754  if (m_escape1 > 0.0) {
755  double arg = (m_pos_escape1-ereco) * m_wgt_escape1;
756  if (std::abs(arg) < 5.0) {
757  value += m_escape1 * std::exp(-0.5 * arg * arg);
758  }
759  }
760 
761  // Add second escape peak if amplitude is positive and
762  // reconstructed energy is within 5 sigma of the Gaussian
763  // position
764  if (m_escape2 > 0.0) {
765  double arg = (m_pos_escape2-ereco) * m_wgt_escape2;
766  if (std::abs(arg) < 5.0) {
767  value += m_escape2 * std::exp(-0.5 * arg * arg);
768  }
769  }
770 
771  // Compute sigma for continua. Since this may not work
772  double sigma = 0.0;
773  if ((m_tail > 0.0) || (m_background > 0.0)) {
774  double e = (ereco > m_emin) ? ereco : m_emin;
775  sigma = m_energies.interpolate(e, m_sigmas);
776  }
777 
778  // Add Compton tail if amplitude is positive
779  if (m_tail > 0.0) {
780 
781  // Setup integration kernel for Compton tail
782  kn_gauss_kernel integrand(ereco, m_position, m_compton_edge, sigma);
783 
784  // Setup integral
785  GIntegral integral(&integrand);
786 
787  // Set precision
788  integral.eps(1.0e-5);
789 
790  // No warnings
791  #if defined(G_UPDATE_RESPONSE_VECTOR_NO_WARNINGS)
792  integral.silent(true);
793  #endif
794 
795  // Perform integration over Gaussian width
796  double emin = ereco - 3.0 * sigma;
797  double emax = ereco + 3.0 * sigma;
798  value += m_tail * integral.romberg(emin, emax);
799 
800  } // endif: added Compton tail
801 
802  // Add Compton background if amplitude is positive
803  if (m_background > 0.0) {
804 
805  // Setup integration kernel for Compton background
806  bkg_gauss_kernel integrand(ereco, m_position, sigma);
807 
808  // Setup integral
809  GIntegral integral(&integrand);
810 
811  // Set precision
812  integral.eps(1.0e-5);
813 
814  // No warnings
815  #if defined(G_UPDATE_RESPONSE_VECTOR_NO_WARNINGS)
816  integral.silent(true);
817  #endif
818 
819  // Perform integration over Gaussian width
820  double emin = ereco - 3.0 * sigma;
821  double emax = ereco + 3.0 * sigma;
822  value += m_background * integral.romberg(emin, emax);
823 
824  } // endif: added Compton background
825 
826  // Store response value
827  m_rsp_values.push_back(value);
828 
829  // Increment normalisation
830  norm += value;
831 
832  } // endfor: looped over response vector bins
833 
834  // Normalise response vector
835  if (norm > 0.0) {
836  norm = 1.0 / (norm * ebin);
837  for (int i = 0; i < nbin; ++i) {
838  m_rsp_values[i] *= norm;
839  }
840  }
841 
842  // Debug
843  #if defined(G_DEBUG_UPDATE_RESPONSE_VECTOR)
844  std::cout << "etrue=" << etrue;
845  std::cout << " ebin=" << ebin;
846  std::cout << " nbin=" << nbin;
847  std::cout << " norm=" << norm << std::endl;
848  #endif
849 
850  } // endif: photo-peak position was positive
851 
852  } // endif: true energy has changed
853 
854  // Return
855  return;
856 }
857 
858 
859 /***********************************************************************//**
860  * @brief Computes modified Klein-Nishina cross section multiplied with
861  * Gaussian
862  *
863  * @param[in] e Energy (MeV).
864  * @return Modified Klein-Nishina cross section multiplied with Gaussian
865  *
866  * Computes the modified Klein-Nishina cross section multiplied with a
867  * Gaussian using
868  *
869  * \f[
870  * KN_{\rm mod}(E|E',E_0) = \sigma_{\rm KN}(E|E_0) \times (1-C(E|E_0))
871  * \times \exp \left( -\frac{1}{2}
872  * \frac{(E-E')^2}{\sigma^2(E')}
873  * \right)
874  * \f]
875  *
876  * where
877  *
878  * \f[
879  * \sigma_{\rm KN}(E|E_0) = \theta(E_c-E)
880  * \left[
881  * \left( \frac{E/E_0}{1-E/E_0}
882  * \frac{m_e c^2}{E_0} \right)^2 -
883  * \frac{E}{E_0} + \frac{1}{1-E/E_0} \right]
884  * \f]
885  *
886  * is the Klein-Nishina cross section, and
887  *
888  * \f[
889  * C(E|E_0) = \theta(E-12.14) \times (1 - e^{-\mu(E|E_0) \, l(E_0)})
890  * \f]
891  *
892  * is a modification, where \f$\theta(x)\f$ is a step function that is 1 for
893  * positive \f$x\f$ and 0 for non-positive \f$x\f$,
894  *
895  * \f[
896  * \mu(E|E_0) = 0.72 \, e^{-1.28 (E_0 - E)^{0.35}} +
897  * 0.01 \, (E_0 - E) +
898  * 0.014 \, (E_0 - E)^{-2.5}
899  * \f]
900  *
901  * is the total linear attenuation coefficient in NaI for all processes
902  * and
903  *
904  * \f[
905  * l(E_0) = 2.9 \log( E_0 - 11.14)
906  * \f]
907  *
908  * is the pathlength in the D2 module (energies are in MeV). \f$E_0\f$ is
909  * the position of the photo peak,
910  *
911  * \f[
912  * E_c = \frac{E_0}{1 + \frac{m_e c^2}{2 E_0}}
913  * \f]
914  *
915  * is the Compton edge, \f$E'\f$ is the reconstructed energy,
916  * \f$\sigma(E')\f$ is the standard deviation at the reconstructed energy,
917  * and \f$E\f$ is the energy over which the convolution is performed.
918  ***************************************************************************/
920 {
921  // Initialise result
922  double value = 0.0;
923 
924  // Continue only if the energy is below the Compton edge
925  if (e < m_ec) {
926 
927  // Compute terms
928  double a = e/m_e0;
929  double term = a/(1.0-a) * gammalib::mec2/m_e0 - 1.0;
930 
931  // Compute Klein-Nishina cross section
932  value = term * term - a + 1.0/(1.0-a);
933 
934  // If the incident energy is above 12.14 MeV then multiply-in the
935  // high-energy correction that was introduced by Rob van Dijk
936  if (m_e0 > 12.14) {
937 
938  // Compute linear attenuation coefficient
939  double d = m_e0 - e;
940  double mu = 0.72 * std::exp(-1.28 * std::pow(d,0.35)) +
941  0.01 * d + 0.014 * std::pow(d,-2.5);
942 
943  // Compute path length through module
944  double l = 2.9 * std::log(m_e0 - 11.14);
945 
946  // Multiply-in correction factor
947  value *= std::exp(-mu * l);
948 
949  } // endif: Multiplied-in high-energy correction
950 
951  // Compute Gaussian
952  double arg = (e - m_ereco) * m_wgt;
953  double gauss = std::exp(-0.5*arg*arg);
954 
955  // Multiply-in Gaussian
956  value *= gauss;
957 
958  } // endif: energy was below the Compton edge
959 
960  // Return value
961  return value;
962 }
963 
964 
965 /***********************************************************************//**
966  * @brief Computes Compton background multiplied with Gaussian
967  *
968  * @param[in] e Energy (MeV).
969  * @return Compton background multiplied with Gaussian
970  *
971  * Computes the Compton background multiplied with a Gaussian using
972  *
973  * \f[
974  * B_{\rm c}(E|E',E_0) = \theta(E_0-E)
975  * \exp \left( -\frac{1}{2}
976  * \frac{(E-E')^2}{\sigma^2(E')}
977  * \right)
978  * \f]
979  *
980  * where \f$\theta(x)\f$ is a step function that is 1 for positive \f$x\f$
981  * and 0 for non-positive \f$x\f$, \f$E_0\f$ is the position of the
982  * photo peak, \f$E'\f$ is the reconstructed energy, \f$\sigma(E')\f$ is
983  * the standard deviation at the reconstructed energy, and \f$E\f$
984  * is the energy over which the convolution is performed.
985  ***************************************************************************/
987 {
988  // Initialise background
989  double value = (e < m_e0) ? 1.0 : 0.0;
990 
991  // Compute Gaussian
992  double arg = (e - m_ereco) * m_wgt;
993  double gauss = std::exp(-0.5*arg*arg);
994 
995  // Multiply-in Gaussian
996  value *= gauss;
997 
998  // Return value
999  return value;
1000 }
const double mec2
Definition: GTools.hpp:52
GCOMD2Response & operator=(const GCOMD2Response &rsp)
Assignment operator.
GCOMD2Response(void)
Void constructor.
void unit(const std::string &unit)
Set column unit.
void clear(void)
Clear instance.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
void copy_members(const GCOMD2Response &rsp)
Copy class members.
FITS table double column class interface definition.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
void read(const GFitsTable &table)
Read COMPTEL D2 module response.
double m_wgt_escape2
Inverse of width of first escape peak (1/MeV)
std::vector< double > m_emins
Lower energy threshold of D2.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
void update_cache(const double &etrue) const
Update computation cache.
double m_sigma
Width of photo peak (MeV)
std::vector< double > m_rsp_values
Response vector values.
double emax(void) const
Return maximum D2 input energy (MeV)
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
double m_pos_escape2
Position of second escape peak (MeV)
std::vector< double > m_escapes2
Amplitude of second escape peak.
COMPTEL D2 module response class interface definition.
void write(GFitsBinTable &table)
Write COMPTEL D2 module response.
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 clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void reserve(const int &num)
Reserves space for nodes in node array.
Definition: GNodeArray.hpp:216
std::vector< double > m_amplitudes
Photo peak amplitude.
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
double m_wgt_photo
Inverse of width of photo peak (1/MeV)
double m_escape2
Amplitude of second escape peak.
double sigma(const double &etrue) const
Return photo peak standard deviation.
FITS table float column class interface definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
FITS file class.
Definition: GFits.hpp:63
double m_escape1
Amplitude of first escape peak.
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL D2 module response information.
void init_members(void)
Initialise class members.
bool is_empty(void) const
Signals if there are no nodes in node array.
Definition: GNodeArray.hpp:202
double m_position
Position of photo peak (MeV)
void free_members(void)
Delete class members.
double m_amplitude
Amplitude of photo peak.
Calibration database class.
Definition: GCaldb.hpp:66
double eval(const double &e)
Computes modified Klein-Nishina cross section multiplied with Gaussian.
Interface for the COMPTEL D2 module response class.
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
double m_tail
Amplitude of Compton tail.
void load(const std::string &sdbname)
Load COMPTEL D2 module response.
double m_ec
Compton edge energy (MeV)
std::vector< double > m_sigmas
Photo peak width in MeV.
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
double m_emax
Upper energy limit of D2 (MeV)
double m_energy
Incident total energy (MeV)
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
double m_emin
Amplitude of Compton background
~GCOMD2Response(void)
Destructor.
std::vector< double > m_escapes1
Amplitude of first escape peak.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
std::vector< double > m_emaxs
Upper energy limit of D2.
double m_ewidth
Lower energy threshold width of D2 (MeV)
double emin(void) const
Return minimum D2 input energy (MeV)
GChatter
Definition: GTypemaps.hpp:33
const GCaldb & caldb(void) const
Return calibration database.
double m_wgt_escape1
Inverse of width of first escape peak (1/MeV)
GCOMD2Response * clone(void) const
Clone instance.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
double m_pos_escape1
Position of first escape peak (MeV)
void silent(const bool &silent)
Set silence flag.
Definition: GIntegral.hpp:243
double m_wgt
Inverse of Gaussian standard deviation (1/MeV)
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
GNodeArray m_energies
Input energies.
std::vector< double > m_ewidths
Lower energy threshold width of D2.
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:206
std::vector< double > m_backgrounds
Amplitude of Compton background.
virtual double real(const int &row, const int &inx=0) const =0
double m_compton_edge
Position of Compton edge (MeV)
void update_response_vector(const double &etrue) const
Update response vector.
FITS binary table class.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
GNodeArray m_rsp_energies
Response vector energies.
std::vector< double > m_positions
Photo peak position in MeV.
FITS binary table class definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
FITS table float column.
GCaldb m_caldb
Calibration database.
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:386
std::vector< double > m_tails
Amplitude of Compton tail.
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
double eval(const double &e)
Computes Compton background multiplied with Gaussian.
Integration class interface definition.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
double m_e0
Incident energy (MeV)
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.
double m_rsp_etrue
True energy of response vector.
double operator()(const double &etrue, const double &ereco) const
D2 module response evaluation operator.