GammaLib  2.0.0
 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-2022 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 //!< Supress integration warnings
45 //#define G_NORMALISE_RESPONSE_VECTOR //!< Normalise response vector
46 
47 /* __ Debug definitions __________________________________________________ */
48 //#define G_DEBUG_UPDATE_RESPONSE_VECTOR
49 
50 /* __ Constants __________________________________________________________ */
51 
52 
53 /*==========================================================================
54  = =
55  = Constructors/destructors =
56  = =
57  ==========================================================================*/
58 
59 /***********************************************************************//**
60  * @brief Void constructor
61  *
62  * Creates an empty COMPTEL D2 module response.
63  ***************************************************************************/
65 {
66  // Initialise members
67  init_members();
68 
69  // Return
70  return;
71 }
72 
73 
74 /***********************************************************************//**
75  * @brief Copy constructor
76  *
77  * @param[in] rsp COMPTEL D2 module response.
78  **************************************************************************/
80 {
81  // Initialise members
82  init_members();
83 
84  // Copy members
85  copy_members(rsp);
86 
87  // Return
88  return;
89 }
90 
91 
92 /***********************************************************************//**
93  * @brief Response constructor
94  *
95  * @param[in] caldb Calibration database.
96  * @param[in] sdbname SDB response name.
97  *
98  * Create COMPTEL D2 module response by loading an SDB file from a
99  * calibration database.
100  ***************************************************************************/
101 GCOMD2Response::GCOMD2Response(const GCaldb& caldb, const std::string& sdbname)
102 {
103  // Initialise members
104  init_members();
105 
106  // Set calibration database
107  this->caldb(caldb);
108 
109  // Load D2 module response
110  this->load(sdbname);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Destructor
119  *
120  * Destroys instance of COMPTEL response object.
121  ***************************************************************************/
123 {
124  // Free members
125  free_members();
126 
127  // Return
128  return;
129 }
130 
131 
132 /*==========================================================================
133  = =
134  = Operators =
135  = =
136  ==========================================================================*/
137 
138 /***********************************************************************//**
139  * @brief Assignment operator
140  *
141  * @param[in] rsp COMPTEL D2 module response.
142  * @return COMPTEL D2 module response.
143  ***************************************************************************/
145 {
146  // Execute only if object is not identical
147  if (this != &rsp) {
148 
149  // Free members
150  free_members();
151 
152  // Initialise members
153  init_members();
154 
155  // Copy members
156  copy_members(rsp);
157 
158  } // endif: object was not identical
159 
160  // Return this object
161  return *this;
162 }
163 
164 
165 /***********************************************************************//**
166  * @brief D2 module response evaluation operator
167  *
168  * @param[in] etrue True energy (MeV).
169  * @param[in] ereco Reconstructed energy (MeV).
170  * @return COMPTEL D2 module response.
171  *
172  * Computes the D2 response for reconstructed energy @p ereco in case that
173  * the true energy was @p etrue. A zero response is returned if the
174  * reconstucted energy is outside the validity limit, defined by
175  *
176  * \f[
177  * [{\tt m\_emin} - 0.75 {\tt m\_ewidth}, {\tt m\_emax}]
178  * \f]
179  *
180  * The code implementation is based on the COMPASS RESD1 function
181  * RESRS209.RESD2.F (release 1.0, 14-OCT-91).
182  ***************************************************************************/
183 double GCOMD2Response::operator()(const double& etrue, const double& ereco) const
184 {
185  // Initialise response and area with zero
186  double response = 0.0;
187 
188  // Continue only if a response was loaded
189  if (!m_energies.is_empty()) {
190 
191  // Update response parameters
192  update_cache(etrue);
193 
194  // Compute response only if reconstructed energy is within validity
195  // limits
196  if ((ereco >= (m_emin - 0.75 * m_ewidth)) && (ereco <= m_emax)) {
197 
198  // Update response vector
199  update_response_vector(etrue);
200 
201  // Get response value from response vector
202  response = m_rsp_energies.interpolate(ereco, m_rsp_values);
203  if (response < 0.0) {
204  response = 0.0;
205  }
206 
207  // Apply Gaussian shaped threshold according to COMPASS function
208  // RESD22. Note that thrinf and thrsig are computed in INITIL.F
209  double thrinf = m_emin + 0.5 * m_ewidth;
210  double thrsig = m_ewidth/(2.0 * std::sqrt(2.0 * gammalib::ln2));
211  if (ereco < thrinf) {
212  double fac = (ereco - thrinf)/thrsig;
213  fac *= fac;
214  if (fac > 50.0) {
215  response = 0.0;
216  }
217  else {
218  response *= std::exp(-0.5 * fac);
219  }
220  }
221 
222  } // endif: reconstructed energy within validity limits
223 
224  } // endif: response was loaded
225 
226  // Return response
227  return response;
228 }
229 
230 
231 /*==========================================================================
232  = =
233  = Public methods =
234  = =
235  ==========================================================================*/
236 
237 /***********************************************************************//**
238  * @brief Clear instance
239  *
240  * Clears COMPTEL D2 module response object by resetting all members to an
241  * initial state. Any information that was present in the object before will
242  * be lost.
243  ***************************************************************************/
245 {
246  // Free class members
247  free_members();
248 
249  // Initialise members
250  init_members();
251 
252  // Return
253  return;
254 }
255 
256 
257 /***********************************************************************//**
258  * @brief Clone instance
259  *
260  * @return Pointer to deep copy of COMPTEL D2 module response.
261  ***************************************************************************/
263 {
264  return new GCOMD2Response(*this);
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Load COMPTEL D2 module response.
270  *
271  * @param[in] sdbname COMPTEL D2 module response name.
272  *
273  * Loads the COMPTEL D2 module response with specified name @p sdbname. The
274  * method first searchs for an appropriate response in the calibration
275  * database. If no appropriate response is found, the method takes the
276  * database root path and response name to build the full path to the
277  * response file, and tries to load the response from these paths.
278  ***************************************************************************/
279 void GCOMD2Response::load(const std::string& sdbname)
280 {
281  // Clear instance but conserve calibration database
282  GCaldb caldb = m_caldb;
283  clear();
284  m_caldb = caldb;
285 
286  // First attempt reading the response using the GCaldb interface
287  GFilename filename = m_caldb.filename("","","SDB","","",sdbname);
288 
289  // If filename is empty then build filename from CALDB root path and
290  // response name
291  if (filename.is_empty()) {
292  filename = gammalib::filepath(m_caldb.rootdir(), sdbname);
293  if (!filename.exists()) {
294  GFilename testname = filename + ".fits";
295  if (testname.exists()) {
296  filename = testname;
297  }
298  }
299  }
300 
301  // Open FITS file
302  GFits fits(filename);
303 
304  // Get SDB table
305  const GFitsTable& sdb = *fits.table(1);
306 
307  // Read SDB
308  read(sdb);
309 
310  // Close SDB FITS file
311  fits.close();
312 
313  // Return
314  return;
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Read COMPTEL D2 module response.
320  *
321  * @param[in] table FITS table.
322  *
323  * Read the COMPTEL D2 module response from a SDB FITS table.
324  ***************************************************************************/
326 {
327  // Initialise COMPTEL D2 module response vectors
328  m_energies.clear();
329  m_positions.clear();
330  m_sigmas.clear();
331  m_amplitudes.clear();
332  m_escapes1.clear();
333  m_escapes2.clear();
334  m_tails.clear();
335  m_backgrounds.clear();
336  m_emins.clear();
337  m_ewidths.clear();
338  m_emaxs.clear();
339 
340  // Extract number of entries in table
341  int num = table.nrows();
342 
343  // If there are entries then read them
344  if (num > 0) {
345 
346  // Get column pointers
347  const GFitsTableCol* ptr_energy = table["ENERGY"];
348  const GFitsTableCol* ptr_position = table["POSITION"];
349  const GFitsTableCol* ptr_sigma = table["WIDTH"];
350  const GFitsTableCol* ptr_amplitude = table["AMPLITUDE"];
351  const GFitsTableCol* ptr_escape1 = table["ESCAPE1"];
352  const GFitsTableCol* ptr_escape2 = table["ESCAPE2"];
353  const GFitsTableCol* ptr_tail = table["TAIL"];
354  const GFitsTableCol* ptr_background = table["BACKGROUND"];
355  const GFitsTableCol* ptr_emin = table["EMIN"];
356  const GFitsTableCol* ptr_ewidth = table["EWIDTH"];
357  const GFitsTableCol* ptr_emax = table["EMAX"];
358 
359  // Copy data from table into vectors
360  for (int i = 0; i < num; ++i) {
361  m_energies.append(ptr_energy->real(i));
362  m_positions.push_back(ptr_position->real(i));
363  m_sigmas.push_back(ptr_sigma->real(i));
364  m_amplitudes.push_back(ptr_amplitude->real(i));
365  m_escapes1.push_back(ptr_escape1->real(i));
366  m_escapes2.push_back(ptr_escape2->real(i));
367  m_tails.push_back(ptr_tail->real(i));
368  m_backgrounds.push_back(ptr_background->real(i));
369  m_emins.push_back(ptr_emin->real(i));
370  m_ewidths.push_back(ptr_ewidth->real(i));
371  m_emaxs.push_back(ptr_emax->real(i));
372  }
373 
374  } // endif: there were entries
375 
376  // Return
377  return;
378 }
379 
380 
381 /***********************************************************************//**
382  * @brief Write COMPTEL D2 module response.
383  *
384  * @param[in] table FITS table.
385  *
386  * Write the COMPTEL D2 module response into a SDB FITS table.
387  ***************************************************************************/
389 {
390  // Set extension name
391  table.extname("SDB");
392 
393  // Set number of entries
394  int num = m_energies.size();
395 
396  // Allocate columns
397  GFitsTableFloatCol col_energy("ENERGY", num);
398  GFitsTableDoubleCol col_position("POSITION", num);
399  GFitsTableDoubleCol col_width("WIDTH", num);
400  GFitsTableDoubleCol col_amplitude("AMPLITUDE", num);
401  GFitsTableDoubleCol col_escape1("ESCAPE1", num);
402  GFitsTableDoubleCol col_escape2("ESCAPE2", num);
403  GFitsTableDoubleCol col_tail("TAIL", num);
404  GFitsTableDoubleCol col_background("BACKGROUND", num);
405  GFitsTableDoubleCol col_emin("EMIN", num);
406  GFitsTableDoubleCol col_ewidth("EWIDTH", num);
407  GFitsTableDoubleCol col_emax("EMAX", num);
408 
409  // Set column units
410  col_energy.unit("MeV");
411  col_position.unit("MeV");
412  col_width.unit("MeV");
413  col_emin.unit("MeV");
414  col_ewidth.unit("MeV");
415  col_emax.unit("MeV");
416 
417  // Copy data from vectors into table
418  for (int i = 0; i < num; ++i) {
419  col_energy(i) = m_energies[i];
420  col_position(i) = m_positions[i];
421  col_width(i) = m_sigmas[i];
422  col_amplitude(i) = m_amplitudes[i];
423  col_escape1(i) = m_escapes1[i];
424  col_escape2(i) = m_escapes2[i];
425  col_tail(i) = m_tails[i];
426  col_background(i) = m_backgrounds[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_escape1);
438  table.append(col_escape2);
439  table.append(col_tail);
440  table.append(col_background);
441  table.append(col_emin);
442  table.append(col_ewidth);
443  table.append(col_emax);
444 
445  // Return
446  return;
447 }
448 
449 
450 /***********************************************************************//**
451  * @brief Print COMPTEL D2 module response information
452  *
453  * @param[in] chatter Chattiness.
454  * @return String containing COMPTEL D2 module response information.
455  ***************************************************************************/
456 std::string GCOMD2Response::print(const GChatter& chatter) const
457 {
458  // Initialise result string
459  std::string result;
460 
461  // Continue only if chatter is not silent
462  if (chatter != SILENT) {
463 
464  // Append header
465  result.append("=== GCOMD2Response ===");
466 
467  // Append D2 module response information
468  result.append("\n"+gammalib::parformat("Energy range"));
469  if (m_energies.size() > 0) {
470  result.append(gammalib::str(m_energies[0])+ " - ");
471  result.append(gammalib::str(m_energies[m_energies.size()-1])+ " MeV");
472  }
473  else {
474  result.append("not defined");
475  }
476  result.append("\n"+gammalib::parformat("Entries"));
477  result.append(gammalib::str(m_energies.size()));
478 
479  // Append calibration database
480  result.append("\n"+m_caldb.print(chatter));
481 
482  // Append information
483 
484  } // endif: chatter was not silent
485 
486  // Return result
487  return result;
488 }
489 
490 
491 /*==========================================================================
492  = =
493  = Private methods =
494  = =
495  ==========================================================================*/
496 
497 /***********************************************************************//**
498  * @brief Initialise class members
499  ***************************************************************************/
501 {
502  // Initialise members
503  m_caldb.clear();
504  m_energies.clear();
505  m_positions.clear();
506  m_sigmas.clear();
507  m_amplitudes.clear();
508  m_escapes1.clear();
509  m_escapes2.clear();
510  m_tails.clear();
511  m_backgrounds.clear();
512  m_emins.clear();
513  m_ewidths.clear();
514  m_emaxs.clear();
515 
516  // Initialise pre-computation cache
517  m_energy = 0.0;
518  m_position = 0.0;
519  m_sigma = 0.0;
520  m_amplitude = 0.0;
521  m_escape1 = 0.0;
522  m_escape2 = 0.0;
523  m_tail = 0.0;
524  m_background = 0.0;
525  m_emin = 0.0;
526  m_ewidth = 0.0;
527  m_emax = 0.0;
528  m_pos_escape1 = 0.0;
529  m_pos_escape2 = 0.0;
530  m_wgt_photo = 0.0;
531  m_wgt_escape1 = 0.0;
532  m_wgt_escape2 = 0.0;
533  m_compton_edge = 0.0;
534 
535  // Pre-computed response vector
536  m_rsp_etrue = 0.0;
538  m_rsp_values.clear();
539 
540  // Return
541  return;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Copy class members
547  *
548  * @param[in] rsp COMPTEL response.
549  ***************************************************************************/
551 {
552  // Copy attributes
553  m_caldb = rsp.m_caldb;
554  m_energies = rsp.m_energies;
555  m_positions = rsp.m_positions;
556  m_sigmas = rsp.m_sigmas;
558  m_escapes1 = rsp.m_escapes1;
559  m_escapes2 = rsp.m_escapes2;
560  m_tails = rsp.m_tails;
562  m_emins = rsp.m_emins;
563  m_ewidths = rsp.m_ewidths;
564  m_emaxs = rsp.m_emaxs;
565 
566  // Copy pre-computation cache
567  m_energy = rsp.m_energy;
568  m_position = rsp.m_position;
569  m_sigma = rsp.m_sigma;
570  m_amplitude = rsp.m_amplitude;
571  m_escape1 = rsp.m_escape1;
572  m_escape2 = rsp.m_escape2;
573  m_tail = rsp.m_tail;
575  m_emin = rsp.m_emin;
576  m_ewidth = rsp.m_ewidth;
577  m_emax = rsp.m_emax;
580  m_wgt_photo = rsp.m_wgt_photo;
584 
585  // Copy pre-computed response vector
586  m_rsp_etrue = rsp.m_rsp_etrue;
589 
590  // Return
591  return;
592 }
593 
594 
595 /***********************************************************************//**
596  * @brief Delete class members
597  ***************************************************************************/
599 {
600  // Return
601  return;
602 }
603 
604 
605 /***********************************************************************//**
606  * @brief Update computation cache
607  *
608  * @param[in] etrue True energy (MeV).
609  *
610  * The method assumes that there is a valid D2 module response.
611  ***************************************************************************/
612 void GCOMD2Response::update_cache(const double& etrue) const
613 {
614  // Update only if the true energy has changed
615  if (etrue != m_energy) {
616 
617  // Set true energy
618  m_energy = etrue;
619 
620  // If true energy is below lowest energy or above largest energy
621  // then set response to zero
622  if ((etrue < m_energies[0]) ||
623  (etrue > m_energies[m_energies.size()-1])) {
624  m_position = 0.0;
625  m_sigma = 0.0;
626  m_amplitude = 0.0;
627  m_escape1 = 0.0;
628  m_escape2 = 0.0;
629  m_tail = 0.0;
630  m_background = 0.0;
631  m_emin = 0.0;
632  m_ewidth = 0.0;
633  m_emax = 0.0;
634  m_pos_escape1 = 0.0;
635  m_pos_escape2 = 0.0;
636  m_wgt_escape1 = 0.0;
637  m_wgt_escape2 = 0.0;
638  m_compton_edge = 0.0;
639  }
640 
641  // ... otherwise interpolate response parameters
642  else {
643 
644  // Interpolate response parameters
655 
656  // Derive escape peak positions
659 
660  // Derive inverse standard deviations of photo and escape peaks
662  if (m_wgt_photo > 0.0) {
663  m_wgt_photo = 1.0 / m_wgt_photo;
664  }
666  if (m_wgt_escape1 > 0.0) {
668  }
670  if (m_wgt_escape2 > 0.0) {
672  }
673 
674  // Derive Compton edge
675  m_compton_edge = m_position / (1.0 + 0.5 * gammalib::mec2 / m_position);
676 
677  // Weight Gaussian amplitudes to assure proper normalization
681 
682  } // endif: true energy is in valid range
683 
684  } // endif: true energy has changed
685 
686  // Return
687  return;
688 }
689 
690 
691 /***********************************************************************//**
692  * @brief Update response vector
693  *
694  * @param[in] etrue True energy (MeV).
695  *
696  * Updates the vector that stores the normalised D2 module response as
697  * function of energy. The vector is computed using
698  *
699  * \f[
700  * R_{\rm D2}(E|E_0) =
701  * B_1 \exp \left( -\frac{1}{2} \frac{(E_0-E)^2}{\sigma^2(E_0)} \right) +
702  * 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) +
703  * 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) +
704  * B_4 \int_{E'} KN_{\rm mod}(E'|E,E_0) dE' +
705  * B_5 \int_{E'} B_{\rm c}(E'|E,E_0) dE'
706  * \f]
707  *
708  * where
709  * \f$B1\f$ is the amplitude of the photo peak,
710  * \f$B2\f$ is the amplitude of the first escape peak,
711  * \f$B3\f$ is the amplitude of the second escape peak,
712  * \f$B4\f$ is the amplitude of the Compton tail,
713  * \f$B5\f$ is the amplitude of the Compton background,
714  * \f$E_0\f$ is the position of the photo peak, and
715  * \f$\sigma(E)\f$ is the energy dependent width of the photo peak. The
716  * constant \f$n\f$ is chosen so that
717  *
718  * \f[
719  * \int_{E} R_{\rm D2}(E|E_0) dE = 1
720  * \f]
721  *
722  * The method assumes that there is a valid D2 module response.
723  *
724  * The code implementation is based on the COMPASS RESRS209 function INITIL.F
725  * (release 1.0, 14-Oct-91).
726  ***************************************************************************/
727 void GCOMD2Response::update_response_vector(const double& etrue) const
728 {
729  // Set constants
730  const double prcs = 1.0e-15;
731 
732  // Update only if the true energy has changed
733  if (etrue != m_rsp_etrue) {
734 
735  // Set true energy
736  m_rsp_etrue = etrue;
737 
738  // Update cache to determine the spectral parameters at the true
739  // energy
740  update_cache(etrue);
741 
742  // Continue only if position is positive
743  if (m_position > 0.0) {
744 
745  // Clear response vectors (this is very very important, otherwise
746  // some old elements reside and we just append to them, and things
747  // get out of order!!!)
749  m_rsp_values.clear();
750 
751  // Initialise response vector
752  double ebin = 0.001 * m_position;
753  int nbin = int(1.35 * m_position / ebin);
754  if (nbin < 1) {
755  nbin = 1;
756  }
757  m_rsp_energies.reserve(nbin);
758  m_rsp_values.reserve(nbin);
759 
760  // Initialise response vector normalisation
761  #if defined(G_NORMALISE_RESPONSE_VECTOR)
762  double norm = 0.0;
763  #endif
764 
765  // Get start bin for convolution of Compton background
766  int istart = int((etrue - 5.0 * m_sigma) / ebin);
767 
768  // Fill response vector
769  for (int i = 0; i < nbin; ++i) {
770 
771  // Compute bin energy
772  double ereco = (i + 0.5) * ebin;
773 
774  // Store bin energy
775  m_rsp_energies.append(ereco);
776 
777  // Initialise response value
778  double value = 0.0;
779 
780  // Add photo peak if amplitude is positive and reconstructed
781  // energy is within 5 sigma of the Gaussian position
782  if (m_amplitude > 0.0) {
783  double arg = (m_position-ereco) * m_wgt_photo;
784  if (std::abs(arg) < 5.0) {
785  value += m_amplitude * std::exp(-0.5 * arg * arg);
786  }
787  }
788 
789  // Add first escape peak if amplitude is positive and
790  // reconstructed energy is within 5 sigma of the Gaussian
791  // position
792  if (m_pos_escape1 > 0.0 && m_escape1 > 0.0) {
793  double arg = (m_pos_escape1-ereco) * m_wgt_escape1;
794  if (std::abs(arg) < 5.0) {
795  value += m_escape1 * std::exp(-0.5 * arg * arg);
796  }
797  }
798 
799  // Add second escape peak if amplitude is positive and
800  // reconstructed energy is within 5 sigma of the Gaussian
801  // position
802  if (m_pos_escape2 > 0.0 && m_escape2 > 0.0) {
803  double arg = (m_pos_escape2-ereco) * m_wgt_escape2;
804  if (std::abs(arg) < 5.0) {
805  value += m_escape2 * std::exp(-0.5 * arg * arg);
806  }
807  }
808 
809  // Compute sigma for continua
810  double sigma = 0.0;
811  if ((m_tail > 0.0) || (m_background > 0.0)) {
812  double e = (ereco > m_emin) ? ereco : m_emin;
813  sigma = m_energies.interpolate(e, m_sigmas);
814  }
815 
816  // Add Compton tail if amplitude is positive
817  if (m_tail > 0.0) {
818 
819  // Compute convolution limits for energy bin
820  double emin = ereco - 3.0 * sigma;
821  double emax = ereco + 3.0 * sigma;
822 
823  // Assure that maximum convolution energy is not above the
824  // Compton edge energy and that minimum convolution energy
825  // is below precision
826  if (emax > m_compton_edge) {
827  emax = m_compton_edge;
828  }
829  if (emin < prcs) {
830  emin = prcs;
831  }
832 
833  // Continue only if the minimum convolution energy is
834  // below the Compton edge energy
835  if (emin <= m_compton_edge - prcs) {
836 
837  // Setup integration kernel for Compton tail
838  kn_gauss_kernel integrand(ereco, m_position, m_compton_edge, sigma);
839 
840  // Setup integral
841  GIntegral integral(&integrand);
842 
843  // Set precision
844  integral.eps(1.0e-5);
845 
846  // No warnings
847  #if defined(G_UPDATE_RESPONSE_VECTOR_NO_WARNINGS)
848  integral.silent(true);
849  #endif
850 
851  // Perform integration over Gaussian width
852  value += m_tail * integral.romberg(emin, emax);
853 
854  } // endif: minimum energy below Compton edge
855 
856  } // endif: added Compton tail
857 
858  // Add Compton background if amplitude is positive
859  if (m_background > 0.0) {
860 
861  // If bin is before part of convolution with photo
862  // peak then set contribution to the background level
863  if (i < istart) {
864  value += m_background;
865  }
866 
867  // ... otherwise convolve Compton background with
868  // photopeak
869  else {
870 
871  // Compute convolution limits
872  double emin = ereco - 3.0 * sigma;
873  double emax = ereco + 3.0 * sigma;
874 
875  // Assure that maximum convolution energy is not above
876  // the true energy and that minimum convolution energy
877  // is below precision
878  if (emax > etrue) {
879  emax = etrue;
880  }
881  if (emin < prcs) {
882  emin = prcs;
883  }
884 
885  // Continue only if the minimum convolution energy is
886  // below the photo peak
887  if (emin <= etrue-prcs) {
888 
889  // Setup integration kernel for Compton background
890  bkg_gauss_kernel integrand(ereco, m_position, sigma);
891 
892  // Setup integral
893  GIntegral integral(&integrand);
894 
895  // Set precision
896  integral.eps(1.0e-5);
897 
898  // No warnings
899  #if defined(G_UPDATE_RESPONSE_VECTOR_NO_WARNINGS)
900  integral.silent(true);
901  #endif
902 
903  // Perform integration over Gaussian width
904  value += m_background * integral.romberg(emin, emax);
905 
906  } // endif: maximum energy below photo peak
907 
908  } // endelse: energy bin significantly before photo peak
909 
910  } // endif: added Compton background
911 
912  // Store response value
913  m_rsp_values.push_back(value);
914 
915  // Increment normalisation
916  #if defined(G_NORMALISE_RESPONSE_VECTOR)
917  norm += value;
918  #endif
919 
920  } // endfor: looped over response vector bins
921 
922  // Normalise response vector
923  #if defined(G_NORMALISE_RESPONSE_VECTOR)
924  if (norm > 0.0) {
925  norm = 1.0 / (norm * ebin);
926  for (int i = 0; i < nbin; ++i) {
927  m_rsp_values[i] *= norm;
928  }
929  }
930  #endif
931 
932  // Debug
933  #if defined(G_DEBUG_UPDATE_RESPONSE_VECTOR)
934  std::cout << "etrue=" << etrue;
935  std::cout << " ebin=" << ebin;
936  std::cout << " nbin=" << nbin;
937  #if defined(G_NORMALISE_RESPONSE_VECTOR)
938  std::cout << " norm=" << norm;
939  #endif
940  std::cout << std::endl;
941  #endif
942 
943  } // endif: photo-peak position was positive
944 
945  } // endif: true energy has changed
946 
947  // Return
948  return;
949 }
950 
951 
952 /***********************************************************************//**
953  * @brief Computes modified Klein-Nishina cross section multiplied with
954  * Gaussian kernel
955  *
956  * @param[in] e Energy (MeV).
957  * @return Modified Klein-Nishina cross section multiplied with Gaussian
958  * kernel
959  *
960  * Computes the modified Klein-Nishina cross section multiplied with a
961  * Gaussian kernel using
962  *
963  * \f[
964  * KN_{\rm mod}(E|E_2,\hat{E_2}) =
965  * \sigma_{\rm KN}(E|\hat{E_2}) \times f(E|\hat{E_2}) \times
966  * \frac{1}{\sigma^2(E_2) \sqrt{2\pi}}
967  * \exp \left( -\frac{1}{2} \frac{(E-E_2)^2}{\sigma^2(E_2)} \right)
968  * \f]
969  *
970  * where
971  *
972  * \f[
973  * \sigma_{\rm KN}(E|\hat{E_2}) =
974  * \left( \frac{E/\hat{E_2}}{1-E/\hat{E_2}}
975  * \frac{m_e c^2}{\hat{E_2}} \right)^2 -
976  * \frac{E}{\hat{E_2}} + \frac{1}{1-E/\hat{E_2}}
977  * \f]
978  *
979  * is the Klein-Nishina cross section, and
980  *
981  * \f[
982  * f(E|\hat{E_2}) = \left \{
983  * \begin{array}{l l}
984  * \displaystyle
985  * \exp \left( -\mu(E|\hat{E_2}) \, l(\hat{E_2}) \right), &
986  * \mbox{if $\hat{E_2} > 12.14$ MeV} \\
987  * \displaystyle
988  * 0, & \mbox{otherwise} \\
989  * \end{array}
990  * \right .
991  * \f]
992  *
993  * is the probability that a photon, which has been Compton scattered, has
994  * no second interaction before escaping (aborption, compton scattering, pair
995  * creation), where
996  *
997  * \f[
998  * \mu(E|\hat{E_2}) = 0.72 \, e^{-1.28 (\hat{E_2} - E)^{0.35}} +
999  * 0.01 \, (\hat{E_2} - E) +
1000  * 0.014 \, (\hat{E_2} - E)^{-2.5}
1001  * \f]
1002  *
1003  * is the total linear attenuation coefficient in NaI for all processes,
1004  * which is an empirical function which describes the values given by
1005  * Harshaw, and
1006  *
1007  * \f[
1008  * l(\hat{E_2}) = 2.9 \log( \hat{E_2} - 11.14)
1009  * \f]
1010  *
1011  * is an empirical path length in the D2 module (energies are in MeV).
1012  * \f$\hat{E_2}\f$ is the position of the photo peak,
1013  * \f$E_2\f$ is the energy deposit measured in D2,
1014  * \f$\sigma(E_2)\f$ is the standard deviation at the measured energy,
1015  * and \f$E\f$ is the energy over which the convolution is performed.
1016  *
1017  * The code implementation is based on the COMPASS RESRS209 function KLNSUB.F
1018  * (release 1.0, 14-Oct-91) and CHANCT.F (release 2.0, 26-JAN-93).
1019  ***************************************************************************/
1021 {
1022  // Compute terms
1023  double a = e/m_e0;
1024  double term = a/(1.0-a) * gammalib::mec2/m_e0 - 1.0;
1025 
1026  // Compute Klein-Nishina cross section
1027  double value = term * term - a + 1.0/(1.0-a);
1028 
1029  // If the incident energy is above 12.14 MeV then multiply-in the
1030  // high-energy correction that was introduced by Rob van Dijk
1031  if (m_e0 > 12.14) {
1032 
1033  // Compute energy difference between photopeak and requested
1034  // energy
1035  double d = m_e0 - e;
1036 
1037  // Set chance to 1 and hence value to zero if energy difference
1038  // is not positive
1039  if (d <= 0.0) {
1040  value = 0.0;
1041  }
1042 
1043  // ... otherwise compute chance that a photon which has been
1044  // Compton scattered has a second interaction before escaping
1045  // (absorption, Compton scattering, pair creation). Note that
1046  // this chance is not a physical chance, it is like a fudge
1047  // factor
1048  else {
1049 
1050  // Compute linear attenuation coefficient
1051  double mu = 0.72 * std::exp(-1.28 * std::pow(d,0.35)) +
1052  0.01 * d + 0.014 * std::pow(d,-2.5);
1053 
1054  // Compute path length through module
1055  double l = 2.9 * std::log(m_e0 - 11.14);
1056 
1057  // Multiply-in correction factor
1058  value *= std::exp(-mu * l);
1059 
1060  }
1061 
1062  } // endif: Multiplied-in high-energy correction
1063 
1064  // Compute Gaussian kernel
1065  double arg = (e - m_ereco) * m_wgt;
1066  double gauss = std::exp(-0.5*arg*arg) * m_wgt * gammalib::inv_sqrt2pi;
1067 
1068  // Multiply-in Gaussian
1069  value *= gauss;
1070 
1071  // Return value
1072  return value;
1073 }
1074 
1075 
1076 /***********************************************************************//**
1077  * @brief Computes Compton background multiplied with Gaussian
1078  *
1079  * @param[in] e Energy (MeV).
1080  * @return Gaussian kernel value.
1081  *
1082  * Computes a Gaussian kernel value using
1083  *
1084  * \f[
1085  * B_{\rm c}(E|E',E_0) = \frac{1}{\sigma(E') \sqrt{2\pi}}
1086  * \exp \left( -\frac{1}{2}
1087  * \frac{(E-E')^2}{\sigma^2(E')}
1088  * \right)
1089  * \f]
1090  *
1091  * where
1092  * \f$E'\f$ is the reconstructed energy,
1093  * \f$\sigma(E')\f$ is the standard deviation at the reconstructed energy,
1094  * and \f$E\f$ is the energy over which the convolution is performed.
1095  *
1096  * The code implementation is based on the COMPASS RESRS209 function FUNC2.F
1097  * (release 1.0, 14-Oct-91).
1098  ***************************************************************************/
1100 {
1101  // Compute Gaussian
1102  double arg = (e - m_ereco) * m_wgt;
1103  double gauss = std::exp(-0.5*arg*arg) * m_wgt * gammalib::inv_sqrt2pi;
1104 
1105  // Return value
1106  return gauss;
1107 }
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:192
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:482
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:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
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:1253
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:524
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:220
const double ln2
Definition: GMath.hpp:45
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:46
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:206
double m_position
Position of photo peak (MeV)
void free_members(void)
Delete class members.
double m_amplitude
Amplitude of photo peak.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
Calibration database class.
Definition: GCaldb.hpp:66
double eval(const double &e)
Computes modified Klein-Nishina cross section multiplied with Gaussian kernel.
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.
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:1274
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:245
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:208
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.
const double inv_sqrt2pi
Definition: GMath.hpp:41
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:657
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:1232
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:393
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:1143
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:1342
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: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.
double m_rsp_etrue
True energy of response vector.
double operator()(const double &etrue, const double &ereco) const
D2 module response evaluation operator.