GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMInstChars.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMInstChars.cpp - COMPTEL Instrument Characteristics class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2019 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 GCOMInstChars.cpp
23  * @brief COMPTEL Instrument Characteristics class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <algorithm>
32 #include "GMath.hpp"
33 #include "GEnergy.hpp"
34 #include "GCOMInstChars.hpp"
35 #include "GCOMSupport.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 
39 /* __ Macros _____________________________________________________________ */
40 
41 /* __ Coding definitions _________________________________________________ */
42 #define TEST_KN
43 
44 /* __ Debug definitions __________________________________________________ */
45 //#define G_DEBUG_READ_SELFVETO
46 
47 /* __ Constants __________________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Void constructor
58  *
59  * Creates empty COMPTEL instrument characteristics.
60  ***************************************************************************/
62 {
63  // Initialise members
64  init_members();
65 
66  // Return
67  return;
68 }
69 
70 
71 /***********************************************************************//**
72  * @brief Copy constructor
73  *
74  * @param[in] ict COMPTEL instrument characteristics.
75  **************************************************************************/
77 {
78  // Initialise members
79  init_members();
80 
81  // Copy members
82  copy_members(ict);
83 
84  // Return
85  return;
86 }
87 
88 
89 /***********************************************************************//**
90  * @brief Response constructor
91  *
92  * @param[in] caldb Calibration database.
93  * @param[in] ictname ICT response name.
94  *
95  * Create COMPTEL instrument characteristics by loading an SDA file from a
96  * calibration database.
97  ***************************************************************************/
98 GCOMInstChars::GCOMInstChars(const GCaldb& caldb, const std::string& ictname)
99 {
100  // Initialise members
101  init_members();
102 
103  // Set calibration database
104  this->caldb(caldb);
105 
106  // Load instrument characteristics
107  this->load(ictname);
108 
109  // Return
110  return;
111 }
112 
113 
114 /***********************************************************************//**
115  * @brief Destructor
116  *
117  * Destroys instance of COMPTEL instrument characteristics.
118  ***************************************************************************/
120 {
121  // Free members
122  free_members();
123 
124  // Return
125  return;
126 }
127 
128 
129 /*==========================================================================
130  = =
131  = Operators =
132  = =
133  ==========================================================================*/
134 
135 /***********************************************************************//**
136  * @brief Assignment operator
137  *
138  * @param[in] ict COMPTEL instrument characteristics.
139  * @return COMPTEL instrument characteristics.
140  ***************************************************************************/
142 {
143  // Execute only if object is not identical
144  if (this != &ict) {
145 
146  // Free members
147  free_members();
148 
149  // Initialise members
150  init_members();
151 
152  // Copy members
153  copy_members(ict);
154 
155  } // endif: object was not identical
156 
157  // Return this object
158  return *this;
159 }
160 
161 
162 /*==========================================================================
163  = =
164  = Public methods =
165  = =
166  ==========================================================================*/
167 
168 /***********************************************************************//**
169  * @brief Clear instance
170  *
171  * Clears COMPTEL instrument characteristics object by resetting all members
172  * to an initial state. Any information that was present in the object before
173  * will be lost.
174  ***************************************************************************/
176 {
177  // Free class members
178  free_members();
179 
180  // Initialise members
181  init_members();
182 
183  // Return
184  return;
185 }
186 
187 
188 /***********************************************************************//**
189  * @brief Clone instance
190  *
191  * @return Pointer to deep copy of COMPTEL instrument characteristics.
192  ***************************************************************************/
194 {
195  return new GCOMInstChars(*this);
196 }
197 
198 
199 /***********************************************************************//**
200  * @brief Load COMPTEL instrument characteristics.
201  *
202  * @param[in] ictname COMPTEL instrument characteristics.
203  *
204  * Loads the COMPTEL instrument characteristics with specified name
205  * @p ictname. The method first searchs for an appropriate response in the
206  * calibration database. If no appropriate response is found, the method
207  * takes the database root path and response name to build the full path to
208  * the response file, and tries to load the response from these paths.
209  ***************************************************************************/
210 void GCOMInstChars::load(const std::string& ictname)
211 {
212  // Clear instance but conserve calibration database
213  GCaldb caldb = m_caldb;
214  clear();
215  m_caldb = caldb;
216 
217  // First attempt reading the response using the GCaldb interface
218  GFilename filename = m_caldb.filename("","","ICT","","",ictname);
219 
220  // If filename is empty then build filename from CALDB root path and
221  // response name
222  if (filename.is_empty()) {
223  filename = gammalib::filepath(m_caldb.rootdir(), ictname);
224  if (!filename.exists()) {
225  GFilename testname = filename + ".fits";
226  if (testname.exists()) {
227  filename = testname;
228  }
229  }
230  }
231 
232  // Open FITS file
233  GFits fits(filename);
234 
235  // Get ICT tables
236  const GFitsTable& d1inter = *fits.table("D1INTER");
237  const GFitsTable& d1dens = *fits.table("D1DENS");
238  const GFitsTable& d1pos = *fits.table("D1POS");
239  const GFitsTable& d1rad = *fits.table("D1RAD");
240  const GFitsTable& d1thick = *fits.table("D1THICK");
241  const GFitsTable& d2inter = *fits.table("D2INTER");
242  const GFitsTable& d2dens = *fits.table("D2DENS");
243  const GFitsTable& d2pos = *fits.table("D2POS");
244  const GFitsTable& d2rad = *fits.table("D2RAD");
245  const GFitsTable& d2thick = *fits.table("D2THICK");
246  const GFitsTable& thbar = *fits.table("THBAR");
247  const GFitsTable& delz = *fits.table("DELZ");
248  const GFitsTable& alu = *fits.table("ALU");
249  const GFitsTable& aldens = *fits.table("ALDENS");
250  const GFitsTable& althick = *fits.table("ALTHICK");
251  const GFitsTable& aboved1 = *fits.table("ABOVED1");
252  const GFitsTable& abdens = *fits.table("ABDENS");
253  const GFitsTable& abthick = *fits.table("ABTHICK");
254  const GFitsTable& vetodome = *fits.table("VETODOME");
255  const GFitsTable& vetodens = *fits.table("VETODENS");
256  const GFitsTable& v1thick = *fits.table("V1THICK");
257  const GFitsTable& vthick = *fits.table("VTHICK");
258  const GFitsTable& selfveto = *fits.table("SELFVETO");
259  const GFitsTable& d1multi = *fits.table("D1MULTI");
260  const GFitsTable& d2multi = *fits.table("D2MULTI");
261 
262  // Read information from ICT tables
270  read_pos(d1pos, m_d1pos_x, m_d1pos_y);
271  read_pos(d2pos, m_d2pos_x, m_d2pos_y);
272  read_selfveto(selfveto);
273 
274  // Read constants
275  m_d1dens = d1dens["DENSITY"]->real(0);
276  m_d1rad = d1rad["RADIUS"]->real(0);
277  m_d1thick = d1thick["THICKNESS"]->real(0);
278  m_d2dens = d2dens["DENSITY"]->real(0);
279  m_d2rad = d2rad["RADIUS"]->real(0);
280  m_d2thick = d2thick["THICKNESS"]->real(0);
281  m_thbar = thbar["ANGLE"]->real(0);
282  m_delz = delz["DISTANCE"]->real(0);
283  m_aldens = aldens["DENSITY"]->real(0);
284  m_althick = althick["THICKNESS"]->real(0);
285  m_abdens = abdens["DENSITY"]->real(0);
286  m_abthick = abthick["THICKNESS"]->real(0);
287  m_vetodens = vetodens["DENSITY"]->real(0);
288  m_v1thick = v1thick["THICKNESS"]->real(0);
289  m_vthick = vthick["THICKNESS"]->real(0);
290 
291  // Close ICT FITS file
292  fits.close();
293 
294  // Return
295  return;
296 }
297 
298 
299 /***********************************************************************//**
300  * @brief Return D1 interaction probability
301  *
302  * @param[in] energy Input photon energy (MeV).
303  * @return D1 interaction probability.
304  *
305  * Computes the D1 interaction probability as function of energy using
306  *
307  * \f[
308  * P(E) = 1 - \exp \left(-\mu_m(E) \, l \right)
309  * \f]
310  *
311  * where
312  * \f$-\mu(E)\f$ is the energy-dependent D1 attenuation coefficient in units
313  * of \f$1/cm\f$ that is interpolated using a log-log interpolation of the
314  * ICT table values, and \f$l\f$ is the thickness of the D1 module in
315  * \f$cm\f$.
316  ***************************************************************************/
317 double GCOMInstChars::prob_D1inter(const double& energy) const
318 {
319  // Initialise probability
320  double prob = 0.0;
321 
322  // Continue only if there are energies
323  if (m_d1inter_energies.size() > 0) {
324 
325  // Get log of energy
326  double logE = std::log(energy);
327 
328  // Get interaction coefficient
329  double logc = m_d1inter_energies.interpolate(logE, m_d1inter_coeffs);
330  double coeff = std::exp(logc) * m_d1thick;
331 
332  // Compute interaction probability
333  prob = 1.0 - std::exp(-coeff);
334 
335  }
336 
337  // Return probability
338  return prob;
339 }
340 
341 
342 /***********************************************************************//**
343  * @brief Return probability that no multihit occured
344  *
345  * @param[in] energy Input photon energy (MeV).
346  * @return Probability that no multihit occured.
347  *
348  * Returns the probability that there is no multihit. The probability is
349  * directly interpolated using a log-log interpolation from the D1 values
350  * that are given in the ICT table. The D2 values are not used.
351  ***************************************************************************/
352 double GCOMInstChars::prob_no_multihit(const double& energy) const
353 {
354  // Initialise probability
355  double prob = 1.0;
356 
357  // Continue only if there are energies
358  if (m_d1multi_energies.size() > 0) {
359 
360  // Get log of energy
361  double logE = std::log(energy);
362 
363  // Compute probability that there is no multihit
365 
366  // Make sure that probability does not exceed 1
367  if (prob > 1.0) {
368  prob = 1.0;
369  }
370 
371  }
372 
373  // Return probability
374  return prob;
375 }
376 
377 
378 /***********************************************************************//**
379  * @brief Return transmission above D1
380  *
381  * @param[in] energy Input photon energy (MeV).
382  * @return Transmission above D1.
383  *
384  * Computes the transmission of material above D1 as function of energy
385  * using
386  *
387  * \f[
388  * T(E) = \exp \left(-\mu(E) l \right)
389  * \f]
390  *
391  * where
392  * \f$\mu(E)\f$ is the energy-dependent interaction coefficient for material
393  * above D1 in units of \f$cm^{-1}\f$ that is interpolated using a log-log
394  * interpolation of the ICT table values, and
395  * \f$l\f$ is the thickness of the material above D1 in \f$cm\f$.
396  ***************************************************************************/
397 double GCOMInstChars::trans_D1(const double& energy) const
398 {
399  // Initialise transmission
400  double transmission = 1.0;
401 
402  // Continue only if there are energies
403  if (m_aboved1_energies.size() > 0) {
404 
405  // Get log of energy
406  double logE = std::log(energy);
407 
408  // Get attenuation coefficient
409  double logc = m_aboved1_energies.interpolate(logE, m_aboved1_coeffs);
410  double coeff = std::exp(logc) * m_abthick;
411 
412  // Compute transmission
413  transmission = std::exp(-coeff);
414 
415  }
416 
417  // Return transmission
418  return transmission;
419 }
420 
421 
422 /***********************************************************************//**
423  * @brief Return V1 veto dome transmission
424  *
425  * @param[in] energy Input photon energy (MeV).
426  * @return V1 veto dome transmission.
427  *
428  * Computes the V1 veto dome transmission as function of energy
429  * using
430  *
431  * \f[
432  * T(E) = \exp \left(-\mu(E) l \right)
433  * \f]
434  *
435  * where
436  * \f$\mu(E)\f$ is the energy-dependent interaction coefficient of the V1
437  * veto dome in units of \f$cm^{-1}\f$ that is interpolated using a log-log
438  * interpolation of the ICT table values, and
439  * \f$l\f$ is the thickness of the V1 veto dome in \f$cm\f$.
440  ***************************************************************************/
441 double GCOMInstChars::trans_V1(const double& energy) const
442 {
443  // Initialise transmission
444  double transmission = 1.0;
445 
446  // Continue only if there are energies
447  if (m_veto_energies.size() > 0) {
448 
449  // Get log of energy
450  double logE = std::log(energy);
451 
452  // Get attenuation coefficient
453  double logc = m_veto_energies.interpolate(logE, m_veto_coeffs);
454  double coeff = std::exp(logc) * m_v1thick;
455 
456  // Compute transmission
457  transmission = std::exp(-coeff);
458 
459  }
460 
461  // Return transmission
462  return transmission;
463 }
464 
465 
466 /***********************************************************************//**
467  * @brief Return D2 interaction probability
468  *
469  * @param[in] energy Input photon energy (MeV).
470  * @param[in] phigeo Geometrical scatter angle (deg).
471  * @return D2 interaction probability.
472  *
473  * Computes the D2 interaction probability as function of energy using
474  *
475  * \f[
476  * P(E) = 1 - \exp \left(\frac{-\mu_m(E) \, l}
477  * {\cos \varphi_{\rm geo}} \right)
478  * \f]
479  *
480  * where
481  * \f$-\mu(E)\f$ is the energy-dependent D2 attenuation coefficient in units
482  * of \f$1/cm\f$ that is interpolated using a log-log interpolation of the
483  * ICT table values, \f$l\f$ is the thickness of the D2 module in \f$cm\f$,
484  * and \f$\varphi_{\rm geo}\f$ is the geometrical scatter angle.
485  ***************************************************************************/
486 double GCOMInstChars::prob_D2inter(const double& energy, const double& phigeo) const
487 {
488  // Initialise probability
489  double prob = 0.0;
490 
491  // Continue only if there are energies
492  if (m_d2inter_energies.size() > 0) {
493 
494  // Compute log of D2 energy deposit
495  double logE2 = std::log(gammalib::com_energy2(energy, phigeo));
496 
497  // Get interaction coefficient
498  double logc = m_d2inter_energies.interpolate(logE2, m_d2inter_coeffs);
499  double length = m_d2thick / std::cos(phigeo * gammalib::deg2rad);
500  double coeff = std::exp(logc) * length;
501 
502  // Compute interaction probability
503  prob = 1.0 - std::exp(-coeff);
504 
505  }
506 
507  // Return probability
508  return prob;
509 }
510 
511 
512 /***********************************************************************//**
513  * @brief Return transmission of material between D1 and D2
514  *
515  * @param[in] energy Input photon energy (MeV).
516  * @param[in] phigeo Geometrical scatter angle (deg).
517  * @return Transmission of material between D1 and D2.
518  *
519  * Computes the transmission of material between D1 and D2 as function of
520  * energy using
521  *
522  * \f[
523  * T(E) = \exp \left(\frac{-\mu(E) \, \rho \, l}
524  * {\cos \varphi_{\rm geo}} \right)
525  * \f]
526  *
527  * where
528  * \f$\mu(E)\f$ is the energy-dependent interaction coefficient of the
529  * material between D1 and D2 in units of \f$cm^{-1}\f$ that is interpolated
530  * using a log-log interpolation of the ICT table values,
531  * \f$l\f$ is the thickness of the material between D1 and D2 in \f$cm\f$,
532  * and \f$\varphi_{\rm geo}\f$ is the geometrical scatter angle.
533  ***************************************************************************/
534 double GCOMInstChars::trans_D2(const double& energy, const double& phigeo) const
535 {
536  // Initialise transmission
537  double transmission = 1.0;
538 
539  // Continue only if there are energies
540  if (m_alu_energies.size() > 0) {
541 
542  // Compute log of D2 energy deposit
543  double logE2 = std::log(gammalib::com_energy2(energy, phigeo));
544 
545  // Get attenuation coefficient
546  double logc = m_alu_energies.interpolate(logE2, m_alu_coeffs);
547  double length = m_althick / std::cos(phigeo * gammalib::deg2rad);
548  double coeff = std::exp(logc) * length;
549 
550  // Compute transmission
551  transmission = std::exp(-coeff);
552 
553  }
554 
555  // Return transmission
556  return transmission;
557 }
558 
559 
560 /***********************************************************************//**
561  * @brief Return V2+V3 veto dome transmission
562  *
563  * @param[in] energy Input photon energy (MeV).
564  * @param[in] phigeo Geometrical scatter angle (deg).
565  * @return V2+V3 veto dome transmission.
566  *
567  * Computes the V2+V3 veto dome transmission as function of energy
568  * using
569  *
570  * \f[
571  * T(E) = \exp \left(\frac{-\mu(E) \, \rho \, l}
572  * {\cos \varphi_{\rm geo}} \right)
573  * \f]
574  *
575  * where
576  * \f$\mu(E)\f$ is the energy-dependent interaction coefficient of the V2
577  * and V3 veto domes in units of \f$cm^{-1}\f$ that is interpolated using a
578  * log-log interpolation of the ICT table values,
579  * \f$l\f$ is the thickness of the V2+V3 veto domes in \f$cm\f$, and
580  * \f$\varphi_{\rm geo}\f$ is the geometrical scatter angle.
581  ***************************************************************************/
582 double GCOMInstChars::trans_V23(const double& energy, const double& phigeo) const
583 {
584  // Initialise transmission
585  double transmission = 1.0;
586 
587  // Continue only if there are energies
588  if (m_veto_energies.size() > 0) {
589 
590  // Compute log of D2 energy deposit
591  double logE2 = std::log(gammalib::com_energy2(energy, phigeo));
592 
593  // Get attenuation coefficient
594  double logc = m_veto_energies.interpolate(logE2, m_veto_coeffs);
595  double length = m_vthick / std::cos(phigeo * gammalib::deg2rad);
596  double coeff = std::exp(logc) * length;
597 
598  // Compute transmission
599  transmission = std::exp(-coeff);
600 
601  }
602 
603  // Return transmission
604  return transmission;
605 }
606 
607 
608 /***********************************************************************//**
609  * @brief Return probability that photon was not self vetoed
610  *
611  * @param[in] energy Input photon energy (MeV).
612  * @param[in] zenith Zenith angle (deg).
613  * @return Probability that photon was not self vetoed.
614  *
615  * Returns the probability that the photon was not self vetoed. The
616  * probability is directly bi-linearly interpolated from the values that
617  * are given in the ICT table.
618  ***************************************************************************/
619 double GCOMInstChars::prob_no_selfveto(const double& energy, const double& zenith) const
620 {
621  // Initialise probability
622  double prob = 1.0;
623 
624  // Get number of energies and zenith angles in array
625  int n_energies = m_selfveto_energies.size();
626  int n_zeniths = m_selfveto_zeniths.size();
627 
628  // Continue only if there are energies and zenith angles
629  if ((n_energies > 0) && (n_zeniths > 0)) {
630 
631  // Set energy interpolation
633 
634  // Set zenith angle interpolation
636 
637  // Set array indices for bi-linear interpolation
638  int inx_zenith_left = m_selfveto_zeniths.inx_left() * n_energies;
639  int inx_zenith_right = m_selfveto_zeniths.inx_right() * n_energies;
640  int inx1 = m_selfveto_energies.inx_left() + inx_zenith_left;
641  int inx2 = m_selfveto_energies.inx_left() + inx_zenith_right;
642  int inx3 = m_selfveto_energies.inx_right() + inx_zenith_left;
643  int inx4 = m_selfveto_energies.inx_right() + inx_zenith_right;
644 
645  // Set weighting factors for bi-linear interpolation
650 
651  // Perform bi-linear interpolation
652  prob = wgt1 * m_selfveto_coeffs[inx1] +
653  wgt2 * m_selfveto_coeffs[inx2] +
654  wgt3 * m_selfveto_coeffs[inx3] +
655  wgt4 * m_selfveto_coeffs[inx4];
656 
657  // Make sure that probability does is within [0,1]
658  if (prob < 0.0) {
659  prob = 0.0;
660  }
661  else if (prob > 1.0) {
662  prob = 1.0;
663  }
664 
665  // And now convert into a no self veto probability
666  prob = 1.0 - prob;
667 
668  }
669 
670  // Return probability
671  return prob;
672 }
673 
674 
675 /***********************************************************************//**
676  * @brief Return multi-scatter correction
677  *
678  * @param[in] energy Input photon energy (MeV).
679  * @param[in] phigeo Geometrical scatter angle (deg).
680  * @return Correction factor due to multi-scatter.
681  *
682  * Returns the fraction of photons that have undergone a single scatter
683  * and which leave the D1 module unattenuated (no second interaction when
684  * traversing the remaining path inside the same module).
685  *
686  * RES already calculates the fraction of photons which undergo a single
687  * scatter for VERTICALLY indicent photons, based on the Klein-Nishina
688  * cross-section and the composition of NE213A. Therefore the above mentioned
689  * transmission can be applied as a multiplicative correction per phigeo.
690  * However, in order to calculate that correction, an integral over the
691  * module must be performed, which properly takes into account the radiative
692  * transfer inside the module at all R and z. We again use the assumption
693  * of vertical photon incidence, which simplifies the calculation. The
694  * integral is done over the location of the first interaction.
695  *
696  * Note that the transmission calculated is conservative : in reality it will
697  * be a bit higher because of a fraction of the photons which undergo a
698  * second scatter have a final scatter angle after escaping the module within
699  * the phigeo-range of the PSF. These events will, however, mainly be located
700  * at large phibar (large D1E deposit).
701  *
702  * The code implementation is based on the COMPASS RESPSIT2 function
703  * MULTIS.F (release ?, 27-NOV-92).
704  ***************************************************************************/
705 double GCOMInstChars::multi_scatter(const double& energy, const double& phigeo) const
706 {
707  // Set integration step size
708  const int nz = 10; //!< Number of steps in z-direction
709  const int nphi = 20; //!< Number of azimuthal steps
710 
711  // Compute stuff
712  double alpha = energy / gammalib::mec2;
713  double sth = std::sin(phigeo * gammalib::deg2rad);
714  double cth = std::cos(phigeo * gammalib::deg2rad);
715 
716  // Get the attenuation coefficient for NE213A
717  double mu = 1.0 / ne213a_mfpath(energy);
718 
719  // Compute the energy after scattering and the corresponding attenuation
720  double e_scattered = energy / (1.0 + alpha * (1.0 - cth));
721  double mu_scattered = 1.0 / ne213a_mfpath(e_scattered);
722 
723  // Compute the full vertical optical depth of the module
724  double tau0 = mu * m_d1thick;
725 
726  // Compute the fraction of photons that are interacting
727  double total_interactions = 1.0 - std::exp(-tau0);
728 
729  // Perform integration per rho-value (symmetry!) over
730  // (1) z (the geometrical depth in the module)
731  // (2) phi (azimuth of scatter)
732  // at geometrical scatter angle phigeo.
733  double deltaz = m_d1thick / double(nz);
734  double deltau = deltaz * mu;
735  double dltphi = gammalib::pi / double(nphi);
736 
737  // Initialise integration loop. The binning in radius is based on
738  // constant surface per cylinder
739  double rho = 1.0;
740  double delrho = 2.0 * rho;
741  double rholow = rho - 0.5 * delrho; // Should be 0
742  double rhoupp = rho + 0.5 * delrho; // Should be 2
743  double surface = 2.0 * delrho * rho; // Should be 4
744 
745  // Compute the integration normalisation value
746  double kappa = deltau / double(nphi);
747 
748  // Integration loop over the module radius
749  double contribution_rho = 0.0;
750  double total_surface = surface;
751  int klast = 0; // Will be >0 if the loop should be exited
752  for (int krho = 0; krho < 100; ++krho) {
753 
754  // If we have reached the end then exit the loop now
755  if (klast > 0) {
756  break;
757  }
758 
759  // Determine radius limits. If D1 module radius is reached or
760  // exceeded, limit the upper rho value to the D1 module radius,
761  // re-compute the surface, and signal to exit the loop in the
762  // next round.
763  if (krho > 0) {
764  rholow = rhoupp;
765  rhoupp = std::sqrt(surface + rholow*rholow);
766  rho = 0.5 * (rhoupp + rholow);
767  if (rhoupp > m_d1rad) {
768  rhoupp = m_d1rad;
769  rho = 0.5 * (rhoupp + rholow);
770  surface = (rhoupp*rhoupp - rholow*rholow);
771  klast = krho;
772  }
773  total_surface += surface;
774  }
775 
776  // Initialise azimuthal results
777  std::vector<double> r(nphi, 0.0);
778 
779  // Compute remaining radius as function of azimuth angle if the
780  // first interaction was at radius rho.
781  for (int lphi = 0; lphi < nphi; ++lphi) {
782  double phi = dltphi * (lphi + 0.5);
783  double term = rho * std::sin(phi);
784  r[lphi] = -rho * std::cos(phi) +
785  std::sqrt(m_d1rad*m_d1rad - term*term);
786  }
787 
788  // Perform integration over depth
789  double contribution_z = 0.0;
790  for (int kz = 0; kz < nz; ++kz) {
791 
792  // Calculate depth
793  double z = (kz + 0.5) * deltaz;
794  double tau = z * mu;
795 
796  // Compute the remaining length in the module after the first
797  // interaction. Test for foward/backward scattering although
798  // this should not be relevant for any reasonable phigeo values.
799  double length0 = 1000.0;
800  if (sth < 0.99) { // True for all reasonable phigeo's
801  if (cth < 0.0) { // False for all reasonable phigeo's
802  length0 = z/std::abs(cth);
803  }
804  else {
805  length0 = (m_d1thick-z)/cth; // All reasonable phigeo's
806  }
807  }
808 
809  // Perform integration over azimuth
810  double contribution_phi = 0.0;
811  for (int kphi = 0; kphi < nphi; ++kphi) {
812 
813  // Compute the actual remaining length for a given azimuth
814  // angle phi. Limit the length to the remaining radius.
815  double length = length0;
816  if (length * sth > r[kphi]) {
817  length = r[kphi] / sth;
818  }
819 
820  // Now compute the probability that the photon is not
821  // absorbed during the remaining path through the detector.
822  contribution_phi += std::exp(-length * mu_scattered);
823 
824  } // endfor: looped over azimuth
825 
826  // Average contribution must later be divided by nphi for
827  // phi-pixels and multiplied by deltau; this is done by
828  // multiplying with kappa. Do this at end, to save computation
829  contribution_z += contribution_phi * std::exp(-tau);
830 
831  } // endfor: looped over depth
832 
833  // Add the average contribution of radius rho
834  contribution_rho += contribution_z * kappa * surface;
835 
836  } // endfor: looped over radius
837 
838  // Compute average
839  double transmission = contribution_rho / (total_surface * total_interactions);
840 
841  // Return transmission
842  return transmission;
843 }
844 
845 
846 /***********************************************************************//**
847  * @brief Return PSD correction
848  *
849  * @param[in] energy Input photon energy (MeV).
850  * @param[in] phigeo Geometrical scatter angle (deg).
851  * @return Correction factor due to PSD selection.
852  *
853  * Returns the D1 energy dependent PSD correction as described in
854  * COM-RP-ROL-DRG-016 and COM-RP-ROL-DRG-35. It applies to a standard PSD
855  * selection of 0-110.
856  *
857  * The acceptance probability fit formula
858  *
859  * \f[
860  * P_{\rm acc} = 1 - \frac{1}{a_1 \times E_1^{a_2} + 1}
861  * \f]
862  *
863  * where \f$a_1=1727.9\f$, \f$a_2=2.53\f$ and \f$E_1\f$ is the D1 energy
864  * deposit in MeV. Coefficients are taken from Boron calibration data
865  * (ZA=1.14) to remain consistent with Rob van Dijk's SIMPSF corrections.
866  *
867  * The code implementation is based on the COMPASS RESPSIT2 function
868  * PSDACP.F (release 1.0, 11-DEC-92).
869  ***************************************************************************/
870 double GCOMInstChars::psd_correction(const double& energy, const double& phigeo) const
871 {
872  // Set constants
873  const double a1 = 1727.9;
874  const double a2 = 2.530;
875 
876  // Compute D1 energy deposit
877  double e1 = gammalib::com_energy1(energy, phigeo);
878 
879  // Original COMPASS code
880  double psdacp = 1.0 - (1.0 / (a1 * std::pow(e1, a2) + 1.0));
881 
882  // Return fraction
883  return psdacp;
884 }
885 
886 
887 /***********************************************************************//**
888  * @brief Print COMPTEL instrument characteristics
889  *
890  * @param[in] chatter Chattiness.
891  * @return String containing COMPTEL instrument characteristics.
892  ***************************************************************************/
893 std::string GCOMInstChars::print(const GChatter& chatter) const
894 {
895  // Initialise result string
896  std::string result;
897 
898  // Continue only if chatter is not silent
899  if (chatter != SILENT) {
900 
901  // Append header
902  result.append("=== GCOMInstChars ===");
903 
904  // Append D1INTER information
905  result.append("\n"+gammalib::parformat("D1 interaction coeffs."));
906  if (m_d1inter_energies.size() > 0) {
907  int last = m_d1inter_energies.size()-1;
908  result.append(gammalib::str(std::exp(m_d1inter_energies[0]))+ " - ");
909  result.append(gammalib::str(std::exp(m_d1inter_energies[last]))+ " MeV");
910  result.append(" [");
912  result.append(" - ");
914  result.append("] cm^2/g");
915  }
916  else {
917  result.append("not defined");
918  }
919 
920  // Append D2INTER information
921  result.append("\n"+gammalib::parformat("D2 interaction coeffs."));
922  if (m_d2inter_energies.size() > 0) {
923  int last = m_d2inter_energies.size()-1;
924  result.append(gammalib::str(std::exp(m_d2inter_energies[0]))+ " - ");
925  result.append(gammalib::str(std::exp(m_d2inter_energies[last]))+ " MeV");
926  result.append(" [");
928  result.append(" - ");
930  result.append("] cm^2/g");
931  }
932  else {
933  result.append("not defined");
934  }
935 
936  // Append ALU information
937  result.append("\n"+gammalib::parformat("Al interaction coeffs."));
938  if (m_alu_energies.size() > 0) {
939  result.append(gammalib::str(std::exp(m_alu_energies[0]))+ " - ");
940  result.append(gammalib::str(std::exp(m_alu_energies[m_alu_energies.size()-1]))+ " MeV");
941  result.append(" [");
942  result.append(gammalib::str(std::exp(min_coeff(m_alu_coeffs))));
943  result.append(" - ");
944  result.append(gammalib::str(std::exp(max_coeff(m_alu_coeffs))));
945  result.append("] 1/cm");
946  }
947  else {
948  result.append("not defined");
949  }
950 
951  // Append ABOVED1 information
952  result.append("\n"+gammalib::parformat("Above D1 atten. coeffs."));
953  if (m_aboved1_energies.size() > 0) {
954  result.append(gammalib::str(std::exp(m_aboved1_energies[0]))+ " - ");
955  result.append(gammalib::str(std::exp(m_aboved1_energies[m_aboved1_energies.size()-1]))+ " MeV");
956  result.append(" [");
958  result.append(" - ");
960  result.append("] 1/cm");
961  }
962  else {
963  result.append("not defined");
964  }
965 
966  // Append VETODOME information
967  result.append("\n"+gammalib::parformat("Vetodome atten. coeffs."));
968  if (m_veto_energies.size() > 0) {
969  result.append(gammalib::str(std::exp(m_veto_energies[0]))+ " - ");
970  result.append(gammalib::str(std::exp(m_veto_energies[m_veto_energies.size()-1]))+ " MeV");
971  result.append(" [");
972  result.append(gammalib::str(std::exp(min_coeff(m_veto_coeffs))));
973  result.append(" - ");
974  result.append(gammalib::str(std::exp(max_coeff(m_veto_coeffs))));
975  result.append("] cm^2/g");
976  }
977  else {
978  result.append("not defined");
979  }
980 
981  // Append SELFVETO information
982  result.append("\n"+gammalib::parformat("Selfveto probabilities"));
983  if (m_selfveto_energies.size() > 0) {
984  result.append(gammalib::str(m_selfveto_energies[0])+ " - ");
985  result.append(gammalib::str(m_selfveto_energies[m_selfveto_energies.size()-1])+ " MeV x ");
986  if (m_selfveto_zeniths.size() > 0) {
987  result.append(gammalib::str(m_selfveto_zeniths[0])+ " - ");
988  result.append(gammalib::str(m_selfveto_zeniths[m_selfveto_zeniths.size()-1])+ " deg");
989  }
990  result.append(" [");
991  result.append(gammalib::str(min_coeff(m_selfveto_coeffs)));
992  result.append(" - ");
993  result.append(gammalib::str(max_coeff(m_selfveto_coeffs)));
994  result.append("]");
995  }
996  else {
997  result.append("not defined");
998  }
999 
1000  // Append D1MULTI information
1001  result.append("\n"+gammalib::parformat("D1 multihit probabilities"));
1002  if (m_d1multi_energies.size() > 0) {
1003  result.append(gammalib::str(std::exp(m_d1multi_energies[0]))+ " - ");
1004  result.append(gammalib::str(std::exp(m_d1multi_energies[m_d1multi_energies.size()-1]))+ " MeV");
1005  result.append(" [");
1007  result.append(" - ");
1009  result.append("]");
1010  }
1011  else {
1012  result.append("not defined");
1013  }
1014 
1015  // Append D2MULTI information
1016  result.append("\n"+gammalib::parformat("D2 multihit probabilities"));
1017  if (m_d2multi_energies.size() > 0) {
1018  result.append(gammalib::str(std::exp(m_d2multi_energies[0]))+ " - ");
1019  result.append(gammalib::str(std::exp(m_d2multi_energies[m_d2multi_energies.size()-1]))+ " MeV");
1020  result.append(" [");
1022  result.append(" - ");
1024  result.append("]");
1025  }
1026  else {
1027  result.append("not defined");
1028  }
1029 
1030  // Append D1POS information
1031  result.append("\n"+gammalib::parformat("D1 positions"));
1032  result.append(gammalib::str(m_d1pos_x.size())+"x & ");
1033  result.append(gammalib::str(m_d1pos_y.size())+"y");
1034 
1035  // Append D2POS information
1036  result.append("\n"+gammalib::parformat("D2 positions"));
1037  result.append(gammalib::str(m_d2pos_x.size())+"x & ");
1038  result.append(gammalib::str(m_d2pos_y.size())+"y");
1039 
1040  // Append D1DENS information
1041  result.append("\n"+gammalib::parformat("D1 density"));
1042  result.append(gammalib::str(m_d1dens)+ " g/cm^-3");
1043 
1044  // Append D1RAD information
1045  result.append("\n"+gammalib::parformat("D1 radius"));
1046  result.append(gammalib::str(m_d1rad)+ " cm");
1047 
1048  // Append D1THICK information
1049  result.append("\n"+gammalib::parformat("D1 thickness"));
1050  result.append(gammalib::str(m_d1thick)+ " cm");
1051 
1052  // Append D2DENS information
1053  result.append("\n"+gammalib::parformat("D2 density"));
1054  result.append(gammalib::str(m_d2dens)+ " g/cm^-3");
1055 
1056  // Append D2RAD information
1057  result.append("\n"+gammalib::parformat("D2 radius"));
1058  result.append(gammalib::str(m_d2rad)+ " cm");
1059 
1060  // Append D2THICK information
1061  result.append("\n"+gammalib::parformat("D2 thickness"));
1062  result.append(gammalib::str(m_d2thick)+ " cm");
1063 
1064  // Append THBAR information
1065  result.append("\n"+gammalib::parformat("Average D2 incident angle"));
1066  result.append(gammalib::str(m_thbar)+ " deg");
1067 
1068  // Append DELZ information
1069  result.append("\n"+gammalib::parformat("Distance between D1 and D2"));
1070  result.append(gammalib::str(m_delz)+ " cm");
1071 
1072  // Append ALDENS information
1073  result.append("\n"+gammalib::parformat("Dens. of Al plate above D2"));
1074  result.append(gammalib::str(m_aldens)+ " g/cm^-3");
1075 
1076  // Append ALTHICK information
1077  result.append("\n"+gammalib::parformat("Thk. of Al plate above D2"));
1078  result.append(gammalib::str(m_althick)+ " cm");
1079 
1080  // Append ABDENS information
1081  result.append("\n"+gammalib::parformat("Density above D1"));
1082  result.append(gammalib::str(m_abdens)+ " g/cm^-3");
1083 
1084  // Append ABTHICK information
1085  result.append("\n"+gammalib::parformat("Thickness above D1"));
1086  result.append(gammalib::str(m_abthick)+ " cm");
1087 
1088  // Append VETODENS information
1089  result.append("\n"+gammalib::parformat("Density of veto domes"));
1090  result.append(gammalib::str(m_vetodens)+ " g/cm^-3");
1091 
1092  // Append V1THICK information
1093  result.append("\n"+gammalib::parformat("Thickness of V1"));
1094  result.append(gammalib::str(m_v1thick)+ " cm");
1095 
1096  // Append VTHICK information
1097  result.append("\n"+gammalib::parformat("Thickness of V2+V3"));
1098  result.append(gammalib::str(m_vthick)+ " cm");
1099 
1100  // Append calibration database
1101  result.append("\n"+m_caldb.print(chatter));
1102 
1103  // Append information
1104 
1105  } // endif: chatter was not silent
1106 
1107  // Return result
1108  return result;
1109 }
1110 
1111 
1112 /*==========================================================================
1113  = =
1114  = Private methods =
1115  = =
1116  ==========================================================================*/
1117 
1118 /***********************************************************************//**
1119  * @brief Initialise class members
1120  ***************************************************************************/
1122 {
1123  // Initialise members
1124  m_caldb.clear();
1126  m_d1inter_coeffs.clear();
1128  m_d2inter_coeffs.clear();
1130  m_alu_coeffs.clear();
1132  m_aboved1_coeffs.clear();
1134  m_veto_coeffs.clear();
1137  m_selfveto_coeffs.clear();
1139  m_d1multi_coeffs.clear();
1141  m_d2multi_coeffs.clear();
1142  m_d1pos_x.clear();
1143  m_d1pos_y.clear();
1144  m_d2pos_x.clear();
1145  m_d2pos_y.clear();
1146  m_d1dens = 0.0;
1147  m_d1rad = 0.0;
1148  m_d1thick = 0.0;
1149  m_d2dens = 0.0;
1150  m_d2rad = 0.0;
1151  m_d2thick = 0.0;
1152  m_thbar = 0.0;
1153  m_delz = 0.0;
1154  m_aldens = 0.0;
1155  m_althick = 0.0;
1156  m_abdens = 0.0;
1157  m_abthick = 0.0;
1158  m_vetodens = 0.0;
1159  m_v1thick = 0.0;
1160  m_vthick = 0.0;
1161 
1162  // Return
1163  return;
1164 }
1165 
1166 
1167 /***********************************************************************//**
1168  * @brief Copy class members
1169  *
1170  * @param[in] ict COMPTEL instrument characteristics.
1171  ***************************************************************************/
1173 {
1174  // Copy attributes
1175  m_caldb = ict.m_caldb;
1181  m_alu_coeffs = ict.m_alu_coeffs;
1193  m_d1pos_x = ict.m_d1pos_x;
1194  m_d1pos_y = ict.m_d1pos_y;
1195  m_d2pos_x = ict.m_d2pos_x;
1196  m_d2pos_y = ict.m_d2pos_y;
1197  m_d1dens = ict.m_d1dens;
1198  m_d1rad = ict.m_d1rad;
1199  m_d1thick = ict.m_d1thick;
1200  m_d2dens = ict.m_d2dens;
1201  m_d2rad = ict.m_d2rad;
1202  m_d2thick = ict.m_d2thick;
1203  m_thbar = ict.m_thbar;
1204  m_delz = ict.m_delz;
1205  m_aldens = ict.m_aldens;
1206  m_althick = ict.m_althick;
1207  m_abdens = ict.m_abdens;
1208  m_abthick = ict.m_abthick;
1209  m_vetodens = ict.m_vetodens;
1210  m_v1thick = ict.m_v1thick;
1211  m_vthick = ict.m_vthick;
1212 
1213  // Return
1214  return;
1215 }
1216 
1217 
1218 /***********************************************************************//**
1219  * @brief Delete class members
1220  ***************************************************************************/
1222 {
1223  // Return
1224  return;
1225 }
1226 
1227 
1228 /***********************************************************************//**
1229  * @brief Read energy dependent coefficients.
1230  *
1231  * @param[in] table FITS table.
1232  * @param[out] energies Energy node array.
1233  * @param[out] coeffs Coefficients.
1234  *
1235  * Read energy dependent coefficients from FITS table and store their natural
1236  * logarithm in the @p energies and @p coeffs vectors.
1237  ***************************************************************************/
1239  GNodeArray& energies,
1240  std::vector<double>& coeffs)
1241 {
1242  // Initialise energies and coefficients
1243  energies.clear();
1244  coeffs.clear();
1245 
1246  // Extract number of entries in table
1247  int num = table.nrows();
1248 
1249  // If there are entries then read them
1250  if (num > 0) {
1251 
1252  // By default use "COEFFICIENT" column for coefficients, but if
1253  // table contains a "PROBABILITY" column then use that column
1254  std::string coeff_name = "COEFFICIENT";
1255  if (table.contains("PROBABILITY")) {
1256  coeff_name = "PROBABILITY";
1257  }
1258 
1259  // Get column pointers
1260  const GFitsTableCol* ptr_energy = table["ENERGY"];
1261  const GFitsTableCol* ptr_coeff = table[coeff_name];
1262 
1263  // Copy data from table into vectors. Skip any non-positive values
1264  // since we store the logarithms for a log-log interpolation.
1265  for (int i = 0; i < num; ++i) {
1266  double energy = ptr_energy->real(i);
1267  double coeff = ptr_coeff->real(i);
1268  if ((energy > 0.0) && (coeff > 0.0)) {
1269  energies.append(std::log(energy));
1270  coeffs.push_back(std::log(coeff));
1271  }
1272  }
1273 
1274  } // endif: there were entries
1275 
1276  // Return
1277  return;
1278 }
1279 
1280 
1281 /***********************************************************************//**
1282  * @brief Read module positions.
1283  *
1284  * @param[in] table FITS table.
1285  * @param[out] x X-positions of module.
1286  * @param[out] y Y-positions of module.
1287  *
1288  * Read module positions from FITS table.
1289  ***************************************************************************/
1290 void GCOMInstChars::read_pos(const GFitsTable& table, std::vector<double>& x,
1291  std::vector<double>& y)
1292 {
1293  // Initialise positions
1294  x.clear();
1295  y.clear();
1296 
1297  // Extract number of entries in table
1298  int num = table.nrows();
1299 
1300  // If there are entries then read them
1301  if (num > 0) {
1302 
1303  // Get column pointers
1304  const GFitsTableCol* ptr_detnum = table["DETNUM"];
1305  const GFitsTableCol* ptr_x = table["X"];
1306  const GFitsTableCol* ptr_y = table["Y"];
1307 
1308  // Copy data from table into vectors
1309  for (int i = 0; i < num; ++i) {
1310  int detnum = ptr_detnum->integer(i) - 1;
1311  x.push_back(ptr_x->real(detnum));
1312  y.push_back(ptr_y->real(detnum));
1313  }
1314 
1315  } // endif: there were entries
1316 
1317  // Return
1318  return;
1319 }
1320 
1321 
1322 /***********************************************************************//**
1323  * @brief Read selfveto coefficients.
1324  *
1325  * @param[in] table FITS table.
1326  *
1327  * Read selfveto coefficients from FITS table. The selfveto coefficients
1328  * depend on energy and zenith angle, and the input data are not necessarily
1329  * regularly sampled. Therefore the method will form a regular grid from
1330  * the provided values and compute the interpolated values from that grid.
1331  ***************************************************************************/
1333 {
1334  // Initialise selfveto vectors
1337  m_selfveto_coeffs.clear();
1338 
1339  // Extract number of entries in table
1340  int num = table.nrows();
1341 
1342  // If there are entries then read them
1343  if (num > 0) {
1344 
1345  // Get column pointers
1346  const GFitsTableCol* ptr_energies = table["ENERGY"];
1347  const GFitsTableCol* ptr_zeniths = table["ZENITH"];
1348  const GFitsTableCol* ptr_coeffs = table["PROBABILITY"];
1349 
1350  // Initialise energies and zenith vectors for sorting before creating
1351  // node arrays
1352  std::vector<double> energies;
1353  std::vector<double> zeniths;
1354 
1355  // Determine all energies and zenith angles in table
1356  for (int i = 0; i < num; ++i) {
1357 
1358  // Add energy if it is not yet in array
1359  double energy = ptr_energies->real(i);
1360  bool add = true;
1361  for (int k = 0; k < energies.size(); ++k) {
1362  if (energies[k] == energy) {
1363  add = false;
1364  break;
1365  }
1366  }
1367  if (add) {
1368  energies.push_back(energy);
1369  }
1370 
1371  // Add zenith angle if it is not yet in array
1372  double zenith = ptr_zeniths->real(i);
1373  add = true;
1374  for (int k = 0; k < m_selfveto_zeniths.size(); ++k) {
1375  if (m_selfveto_zeniths[k] == zenith) {
1376  add = false;
1377  break;
1378  }
1379  }
1380  if (add) {
1381  zeniths.push_back(zenith);
1382  }
1383 
1384  } // endfor: collected all energies and zenith angles
1385 
1386  // Sort energies and zenith angles
1387  std::sort(energies.begin(), energies.end());
1388  std::sort(zeniths.begin(), zeniths.end());
1389 
1390  // Set node arrays
1391  m_selfveto_energies.nodes(energies);
1392  m_selfveto_zeniths.nodes(zeniths);
1393 
1394  // Get array size
1395  int neng = m_selfveto_energies.size();
1396  int nzenith = m_selfveto_zeniths.size();
1397 
1398  // Fill all coefficients with -1, meaning that there is no value
1399  m_selfveto_coeffs.assign(neng*nzenith, -1.0);
1400 
1401  // Fill all coefficients
1402  for (int i = 0; i < num; ++i) {
1403 
1404  // Find nearest energy
1405  int ieng = 0;
1406  double energy = ptr_energies->real(i);
1407  double delta = std::abs(energy - m_selfveto_energies[0]);
1408  for (int k = 1; k < neng; ++k) {
1409  double d = std::abs(energy - m_selfveto_energies[k]);
1410  if (d < delta) {
1411  delta = d;
1412  ieng = k;
1413  }
1414  }
1415 
1416  // Find nearest zenith angle
1417  int izenith = 0;
1418  double zenith = ptr_zeniths->real(i);
1419  delta = std::abs(zenith - m_selfveto_zeniths[0]);
1420  for (int k = 1; k < nzenith; ++k) {
1421  double d = std::abs(zenith - m_selfveto_zeniths[k]);
1422  if (d < delta) {
1423  delta = d;
1424  izenith = k;
1425  }
1426  }
1427 
1428  // Compute index
1429  int index = ieng + neng*izenith;
1430 
1431  // Set coefficient
1432  m_selfveto_coeffs[index] = ptr_coeffs->real(i);
1433 
1434  } // endfor: filled all coefficients
1435 
1436  // Now loop over all energies and fill the missing zenith angles
1437  for (int ieng = 0; ieng < neng; ++ieng) {
1438 
1439  // Debug: print header
1440  #if defined(G_DEBUG_READ_SELFVETO)
1441  std::cout << "Energy: " << m_selfveto_energies[ieng] << std::endl;
1442  std::cout << " Zenith angles:";
1443  #endif
1444 
1445  // Set up arrays for interpolation
1446  GNodeArray nodes;
1447  std::vector<double> values;
1448  for (int izenith = 0; izenith < nzenith; ++izenith) {
1449  int index = ieng + neng*izenith;
1450  if (m_selfveto_coeffs[index] > -0.5) {
1451  nodes.append(m_selfveto_zeniths[izenith]);
1452  values.push_back(m_selfveto_coeffs[index]);
1453  #if defined(G_DEBUG_READ_SELFVETO)
1454  std::cout << " " << m_selfveto_zeniths[izenith];
1455  std::cout << "(" << m_selfveto_coeffs[index] << ")";
1456  #endif
1457  }
1458  }
1459 
1460  // Debug: print header
1461  #if defined(G_DEBUG_READ_SELFVETO)
1462  std::cout << std::endl;
1463  std::cout << " Interpolated zenith angles:";
1464  #endif
1465 
1466  // Interpolate coefficients for zenith angles
1467  for (int izenith = 0; izenith < nzenith; ++izenith) {
1468  int index = ieng + neng*izenith;
1469  if (m_selfveto_coeffs[index] < -0.5) {
1470  double coeffs = nodes.interpolate(m_selfveto_zeniths[izenith], values);
1471  if (coeffs < 0.0) {
1472  coeffs = 0.0;
1473  }
1474  else if (coeffs > 1.0) {
1475  coeffs = 1.0;
1476  }
1477  m_selfveto_coeffs[index] = coeffs;
1478  #if defined(G_DEBUG_READ_SELFVETO)
1479  std::cout << " " << m_selfveto_zeniths[izenith];
1480  std::cout << "(" << m_selfveto_coeffs[index] << ")";
1481  #endif
1482  }
1483  }
1484 
1485  // Debug: print line feed
1486  #if defined(G_DEBUG_READ_SELFVETO)
1487  std::cout << std::endl;
1488  #endif
1489 
1490  } // endfor: looped over all energies
1491 
1492  } // endif: there were entries
1493 
1494  // Return
1495  return;
1496 }
1497 
1498 
1499 /***********************************************************************//**
1500  * @brief Return NE213A mean free path
1501  *
1502  * @param[in] energy Energy (MeV)
1503  * @return Mean free path (cm)
1504  *
1505  * The mean free path as a function of energy, linearly interpolated between
1506  * data-points. Table of mean free paths for NE213A taken from L. Kuiper
1507  * calculated in JOB ROL-SPDSD1-43. This includes (1) Compton scattering and
1508  * (2) Pair production.
1509  ***************************************************************************/
1510 double GCOMInstChars::ne213a_mfpath(const double& energy) const
1511 {
1512  // NE213A mean free path values for energies from 1,2,...,30 MeV
1513  const double mfpath[] = {16.19, 23.21, 29.01, 34.04, 38.46,
1514  42.35, 45.64, 48.81, 51.53, 54.17,
1515  56.15, 58.09, 59.96, 61.74, 63.41,
1516  64.62, 65.80, 66.94, 68.04, 69.09,
1517  69.86, 70.61, 71.33, 72.03, 72.70,
1518  73.33, 73.93, 74.50, 75.03, 75.51};
1519 
1520  // Normalisation constant for low-energy extrapolation
1521  const double norm = mfpath[1] * kn_cross_section(2.0 / gammalib::mec2);
1522 
1523  // Initialise mfpath
1524  double result = 0.0;
1525 
1526  // Do a Klein-Nishina extrapolation for low energies
1527  if (energy < 2.0) {
1528  result = norm / kn_cross_section(energy / gammalib::mec2);
1529  }
1530 
1531  // ... otherwise do a linear interpolation
1532  else {
1533  int index = int(energy) - 1;
1534  if (index >= 29) {
1535  index = 28;
1536  }
1537  double slope = mfpath[index+1] - mfpath[index];
1538  result = mfpath[index] + slope * (energy - double(index + 1));
1539  }
1540 
1541  // Return result
1542  return result;
1543 }
1544 
1545 
1546 /***********************************************************************//**
1547  * @brief Return integrated Klein-Nishina cross section
1548  *
1549  * @param[in] k \f$E/ m_e c^2\f$
1550  * @return Integrated Klein-Nishina cross section
1551  *
1552  * Computes
1553  *
1554  * \f[
1555  * \sigma_{\rm KN}(k) = \frac{1+k}{k^2}
1556  * \left[
1557  * \frac{2(1+k)}{1+2k} - \frac{\ln(1+2k)}{k}
1558  * \right] + \frac{\ln(1+2k)}{2k} -
1559  * \frac{1+3k}{(1+2k)^2}
1560  * \f]
1561  *
1562  * where \f$k = E/ m_e c^2\f$.
1563  ***************************************************************************/
1564 double GCOMInstChars::kn_cross_section(const double& k) const
1565 {
1566  // Compute integrated Klein-Nishina cross section
1567  double k1 = k + 1.0;
1568  double k2p1 = 1.0 + 2.0*k;
1569  double logk2p1 = std::log(k2p1);
1570  double kn = k1/(k*k) * (2.0*k1/k2p1 - logk2p1/k) +
1571  logk2p1/(2.0*k) - (1.0+3.0*k)/(k2p1*k2p1);
1572 
1573  // Return
1574  return kn;
1575 }
1576 
1577 
1578 /***********************************************************************//**
1579  * @brief Return minimum coefficient.
1580  *
1581  * @param[in] coeffs Coefficients.
1582  * @return Minimum coefficient.
1583  *
1584  * Returns minimum coefficient.
1585  ***************************************************************************/
1586 double GCOMInstChars::min_coeff(const std::vector<double>& coeffs) const
1587 {
1588  // Initialise minimum
1589  double min = 0.0;
1590 
1591  // Continue only if there are coefficients
1592  if (coeffs.size() > 0) {
1593 
1594  // Set first value as minimum
1595  min = coeffs[0];
1596 
1597  // Search for minimum
1598  for (int i = 1; i < coeffs.size(); ++i) {
1599  if (coeffs[i] < min) {
1600  min = coeffs[i];
1601  }
1602  }
1603 
1604  }
1605 
1606  // Return minimum
1607  return min;
1608 }
1609 
1610 
1611 /***********************************************************************//**
1612  * @brief Return maximum coefficient.
1613  *
1614  * @param[in] coeffs Coefficients.
1615  * @return Maximum coefficient.
1616  *
1617  * Returns maximum coefficient.
1618  ***************************************************************************/
1619 double GCOMInstChars::max_coeff(const std::vector<double>& coeffs) const
1620 {
1621  // Initialise maximum
1622  double max = 0.0;
1623 
1624  // Continue only if there are coefficients
1625  if (coeffs.size() > 0) {
1626 
1627  // Set first value as maximum
1628  max = coeffs[0];
1629 
1630  // Search for maximum
1631  for (int i = 1; i < coeffs.size(); ++i) {
1632  if (coeffs[i] > max) {
1633  max = coeffs[i];
1634  }
1635  }
1636 
1637  }
1638 
1639  // Return maximum
1640  return max;
1641 }
const double mec2
Definition: GTools.hpp:52
COMPTEL Instrument Characteristics interface definition.
double m_delz
Distance between D1 and D2 levels (cm)
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
std::vector< double > m_d2pos_x
D2 x-position (cm)
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
Node array class.
Definition: GNodeArray.hpp:60
void read_selfveto(const GFitsTable &table)
Read selfveto coefficients.
double prob_no_multihit(const double &energy) const
Return probability that no multihit occured.
Energy value class definition.
std::vector< double > m_d2inter_coeffs
D2 interaction coefficients.
GNodeArray m_d2inter_energies
D2 interaction coefficient energies (MeV)
double trans_D2(const double &energy, const double &phigeo) const
Return transmission of material between D1 and D2.
double prob_D2inter(const double &energy, const double &phigeo) const
Return D2 interaction probability.
const double pi
Definition: GMath.hpp:35
double trans_D1(const double &energy) const
Return transmission above D1.
double max_coeff(const std::vector< double > &coeffs) const
Return maximum coefficient.
double m_abdens
Density above D1 (g/cm^-3)
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
double prob_D1inter(const double &energy) const
Return D1 interaction probability.
std::vector< double > m_d2multi_coeffs
D2 multihit attenuation coefficients (probability)
GNodeArray m_alu_energies
Al interaction coefficient energies (MeV)
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double m_vetodens
Density of veto domes (g/cm^-3)
GCaldb m_caldb
Calibration database.
double com_energy2(const double &energy, const double &phigeo)
Return D2 energy deposit.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
bool contains(const std::string &colname) const
Checks the presence of a column in table.
Definition: GFitsTable.cpp:716
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
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
void free_members(void)
Delete class members.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void read_pos(const GFitsTable &table, std::vector< double > &x, std::vector< double > &y)
Read module positions.
Interface for the COMPTEL Instrument Characteristics class.
GNodeArray m_aboved1_energies
Above D1 attenuation coefficient energies (MeV)
GNodeArray m_d1multi_energies
D1 multihit attenuation coefficient energies (MeV)
double m_aldens
Density of aluminium plate above D2 (g/cm^-3)
FITS file class.
Definition: GFits.hpp:63
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
double min(const GVector &vector)
Computes minimum vector element.
Definition: GVector.cpp:843
double m_d2dens
D2 density (g/cm^-3)
void copy_members(const GCOMInstChars &ict)
Copy class members.
Implementation of support function used by COMPTEL classes.
double m_vthick
Thickness of V2 and V3 veto domes together (cm)
GCOMInstChars * clone(void) const
Clone instance.
double psd_correction(const double &energy, const double &phigeo) const
Return PSD correction.
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
std::vector< double > m_d1pos_y
D1 y-position (cm)
double ne213a_mfpath(const double &energy) const
Return NE213A mean free path.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
Calibration database class.
Definition: GCaldb.hpp:66
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
std::vector< double > m_d2pos_y
D2 y-position (cm)
double m_abthick
Thickness above D1 (cm)
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
double m_v1thick
Thickness of V1 veto dome (cm)
void load(const std::string &ictname)
Load COMPTEL instrument characteristics.
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
void read_coeffs(const GFitsTable &table, GNodeArray &energies, std::vector< double > &coeffs)
Read energy dependent coefficients.
std::vector< double > m_d1pos_x
D1 x-position (cm)
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
double m_d1dens
D1 density (g/cm^-3)
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL instrument characteristics.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void init_members(void)
Initialise class members.
GNodeArray m_selfveto_energies
Selfveto energies (MeV)
GNodeArray m_d2multi_energies
D2 multihit attenuation coefficient energies (MeV)
GChatter
Definition: GTypemaps.hpp:33
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
const GCaldb & caldb(void) const
Return calibration database.
double m_d2thick
D2 thickness (cm)
std::vector< double > m_d1multi_coeffs
D1 multihit attenuation coefficients (probability)
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:872
void nodes(const int &num, const double *array)
Set node array.
Definition: GNodeArray.cpp:325
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
std::vector< double > m_aboved1_coeffs
Above D1 attenuation coefficients.
GNodeArray m_selfveto_zeniths
Selfveto zenith angle (deg)
double prob_no_selfveto(const double &energy, const double &zenith) const
Return probability that photon was not self vetoed.
GCOMInstChars(void)
Void constructor.
virtual int integer(const int &row, const int &inx=0) const =0
virtual double real(const int &row, const int &inx=0) const =0
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
double m_d1thick
D1 thickness (cm)
void clear(void)
Clear instance.
double m_d1rad
D1 radius (cm)
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
double kn_cross_section(const double &k) const
Return integrated Klein-Nishina cross section.
std::vector< double > m_selfveto_coeffs
Selfveto coefficients (probability)
double m_althick
Thickness of aluminium plate above D2 (cm)
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:386
double trans_V1(const double &energy) const
Return V1 veto dome transmission.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
double m_thbar
Average D2 incident angle (deg)
double m_d2rad
D2 radius (cm)
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
double com_energy1(const double &energy, const double &phigeo)
Return D1 energy deposit.
GNodeArray m_d1inter_energies
D1 interaction coefficient energies (MeV)
double trans_V23(const double &energy, const double &phigeo) const
Return V2+V3 veto dome transmission.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
std::vector< double > m_alu_coeffs
Al interaction coefficients.
GCOMInstChars & operator=(const GCOMInstChars &ict)
Assignment operator.
std::vector< double > m_d1inter_coeffs
D1 interaction coefficients.
void clear(void)
Clear calibration database.
Definition: GCaldb.cpp:194
std::vector< double > m_veto_coeffs
Veto dome attenuation coefficients.
double min_coeff(const std::vector< double > &coeffs) const
Return minimum coefficient.
std::string rootdir(void) const
Return path to CALDB root directory.
Definition: GCaldb.cpp:259
GNodeArray m_veto_energies
Veto dome attenuation coefficient energies (MeV)
Mathematical function definitions.
double multi_scatter(const double &energy, const double &phigeo) const
Return multi-scatter correction.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
~GCOMInstChars(void)
Destructor.