GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMIaq.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMIaq.cpp - COMPTEL instrument response representation *
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 GCOMIaq.cpp
23  * @brief COMPTEL instrument response representation class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GFits.hpp"
33 #include "GFilename.hpp"
34 #include "GIntegral.hpp"
35 #include "GEbounds.hpp"
36 #include "GModelSpectral.hpp"
37 #include "GModelSpectralPlaw.hpp"
40 #include "GCOMIaq.hpp"
41 #include "GCOMSupport.hpp"
42 
43 /* __ Method name definitions ____________________________________________ */
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 //#define G_RESPONSE_KERNEL_LIMITS
49 //#define G_APPLY_PSD_CORRECTION_IN_RESPONSE_KERNEL
50 #define G_COMPUTE_IAQ_BIN_NO_WARNINGS
51 #define G_LOCATION_SMEARING_NO_WARNINGS
52 
53 /* __ Debug definitions __________________________________________________ */
54 //#define G_DEBUG_COMPTON_KINEMATICS
55 //#define G_DEBUG_COMPUTE_IAQ_BIN
56 //#define G_DEBUG_KLEIN_NISHINA
57 //#define G_DEBUG_KLEIN_NISHINA_INTEGRAL
58 //#define G_DEBUG_RESPONSE_KERNEL
59 //#define G_DEBUG_WEIGHT_IAQ
60 //#define G_DEBUG_SET_CONTINUUM
61 //#define G_DEBUG_LOCATION_SMEARING
62 //#define G_DEBUG_LOCATION_SPREAD
63 
64 /* __ Constants __________________________________________________________ */
65 
66 
67 /*==========================================================================
68  = =
69  = Constructors/destructors =
70  = =
71  ==========================================================================*/
72 
73 /***********************************************************************//**
74  * @brief Void constructor
75  *
76  * Creates an empty COMPTEL instrument response representation.
77  ***************************************************************************/
79 {
80  // Initialise members
81  init_members();
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Copy constructor
90  *
91  * @param[in] iaq COMPTEL instrument response representation.
92  **************************************************************************/
94 {
95  // Initialise members
96  init_members();
97 
98  // Copy members
99  copy_members(iaq);
100 
101  // Return
102  return;
103 }
104 
105 
106 /***********************************************************************//**
107  * @brief COMPTEL instrument response representation constructor
108  *
109  * @param[in] phigeo_max Maximum geometrical scatter angle (deg).
110  * @param[in] phigeo_bin_size Bin size in geometrical scatter angle (deg).
111  * @param[in] phibar_max Maximum Compton scatter angle (deg).
112  * @param[in] phibar_bin_size Bin size in Compton scatter angle (deg).
113  *
114  * @todo Test input argument validity.
115  **************************************************************************/
116 GCOMIaq::GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
117  const double& phibar_max, const double& phibar_bin_size)
118 {
119  // Initialise members
120  init_members();
121 
122  // Compute the number of bins in phigeo and phibar
123  int nphigeo = int(phigeo_max / phigeo_bin_size + 0.5);
124  int nphibar = int(phibar_max / phibar_bin_size + 0.5);
125 
126  // Store the bin definition. Use the number of bins to re-derive the
127  // maximum angles
128  m_phigeo_bin_size = phigeo_bin_size;
129  m_phibar_bin_size = phibar_bin_size;
130  m_phigeo_max = phigeo_bin_size * double(nphigeo);
131  m_phibar_max = phibar_bin_size * double(nphibar);
132 
133  // Allocate FITS image
134  m_iaq = GFitsImageFloat(nphigeo, nphibar);
135 
136  // Set FITS image keywords
137  m_iaq.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
138  m_iaq.card("CRVAL1", 0.5*m_phigeo_bin_size, "[deg] Geometrical scatter angle for reference bin");
139  m_iaq.card("CDELT1", m_phigeo_bin_size, "[deg] Geometrical scatter angle bin size");
140  m_iaq.card("CRPIX1", 1.0, "Reference bin of geometrical scatter angle");
141  m_iaq.card("CTYPE2", "Phibar", "Compton scatter angle");
142  m_iaq.card("CRVAL2", 0.5*m_phibar_bin_size, "[deg] Compton scatter angle for reference bin");
143  m_iaq.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
144  m_iaq.card("CRPIX2", 1.0, "Reference bin of Compton scatter angle");
145  m_iaq.card("BUNIT", "Probability per bin", "Unit of bins");
146  m_iaq.card("TELESCOP", "GRO", "Name of telescope");
147  m_iaq.card("INSTRUME", "COMPTEL", "Name of instrument");
148  m_iaq.card("ORIGIN", "GammaLib", "Origin of FITS file");
149  m_iaq.card("OBSERVER", "Unknown", "Observer that created FITS file");
150 
151  // Return
152  return;
153 }
154 
155 
156 /***********************************************************************//**
157  * @brief Destructor
158  *
159  * Destroys instance of COMPTEL response object.
160  ***************************************************************************/
162 {
163  // Free members
164  free_members();
165 
166  // Return
167  return;
168 }
169 
170 
171 /*==========================================================================
172  = =
173  = Operators =
174  = =
175  ==========================================================================*/
176 
177 /***********************************************************************//**
178  * @brief Assignment operator
179  *
180  * @param[in] iaq COMPTEL instrument response representation.
181  * @return COMPTEL instrument response representation.
182  ***************************************************************************/
184 {
185  // Execute only if object is not identical
186  if (this != &iaq) {
187 
188  // Free members
189  free_members();
190 
191  // Initialise members
192  init_members();
193 
194  // Copy members
195  copy_members(iaq);
196 
197  } // endif: object was not identical
198 
199  // Return this object
200  return *this;
201 }
202 
203 
204 /*==========================================================================
205  = =
206  = Public methods =
207  = =
208  ==========================================================================*/
209 
210 /***********************************************************************//**
211  * @brief Clear instance
212  *
213  * Clears COMPTEL instrument response representation by resetting all members
214  * to an initial state. Any information that was present in the object before
215  * will be lost.
216  ***************************************************************************/
217 void GCOMIaq::clear(void)
218 {
219  // Free class members
220  free_members();
221 
222  // Initialise members
223  init_members();
224 
225  // Return
226  return;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Clone instance
232  *
233  * @return Pointer to deep copy of COMPTEL instrument response
234  * representation.
235  ***************************************************************************/
237 {
238  return new GCOMIaq(*this);
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Set mono-energetic IAQ
244  *
245  * @param[in] energy Input photon energy.
246  * @param[in] ebounds Boundaries of observed energy.
247  *
248  * The code implemented is based on the COMPASS RESPSIT2 function SPCIAQ.F
249  * (release 1.0, 18-DEC-92).
250  *
251  * @todo Implement geometrical smearing.
252  ***************************************************************************/
253 void GCOMIaq::set(const GEnergy& energy, const GEbounds& ebounds)
254 {
255  // Store the energy boundaries
256  m_ebounds = ebounds;
257 
258  // Initialise COMPTEL response information and instrument characteristics
259  init_response();
260 
261  // Remove any extra header cards
262  remove_cards();
263 
264  // Generate IAQ using Compton kinematics
265  compton_kinematics(energy.MeV());
266 
267  // Multiply IAQ with Klein-Nishina factors
268  klein_nishina(energy.MeV());
269 
270  // Weight IAQ
271  weight_iaq(energy.MeV());
272 
273  // Apply location smearing
275 
276  // Set source photon energy
277  m_iaq.card("ENERGY", energy.MeV(), "[MeV] Source photon energy");
278 
279  // Return
280  return;
281 }
282 
283 
284 /***********************************************************************//**
285  * @brief Set continuum IAQ
286  *
287  * @param[in] spectrum Input spectrum.
288  * @param[in] ebounds Boundaries of observed energy.
289  *
290  * Computes the continuum IAQ based on an input @p spectrum. The method
291  * computes the line IAQ for a m_num_energies logarithmically spaced energies
292  * within an energy range that conceivably covers the observed energy band.
293  *
294  * The code implementation is loosly based on the COMPASS RESPSIT2 functions
295  * SPCIAQ.F, GENWEI.F and DERFAN.F (release 1.0, 18-DEC-92). The code was
296  * much simplified.
297  *
298  * @todo Implement geometrical smearing.
299  ***************************************************************************/
300 void GCOMIaq::set(const GModelSpectral& spectrum, const GEbounds& ebounds)
301 {
302  // Store the energy boundaries
303  m_ebounds = ebounds;
304 
305  // Initialise COMPTEL response information and instrument characteristics
306  init_response();
307 
308  // Remove any extra header cards
309  remove_cards();
310 
311  // Get minimum and maximum input energy for continuum IAQ
312  double energy_min = 0.7 * m_ebounds.emin().MeV();
313  double energy_max = (m_response_d1.emax() < m_response_d2.emax()) ?
315 
316  // Debug
317  #if defined(G_DEBUG_SET_CONTINUUM)
318  std::cout << "Minimum energy=" << energy_min << std::endl;
319  std::cout << "Maximum energy=" << energy_max << std::endl;
320  #endif
321 
322  // Setup array of energy boundaries
323  GEbounds ebds(m_num_energies, GEnergy(energy_min, "MeV"),
324  GEnergy(energy_max, "MeV"));
325 
326  // Compute flux over measured total energy interval
327  double flux_total = spectrum.flux(m_ebounds.emin(), m_ebounds.emax());
328 
329  // Store empty copy of IAQ
330  GFitsImageFloat iaq = m_iaq;
331  for (int k = 0; k < iaq.npix(); ++k) {
332  iaq(k) = 0.0;
333  }
334 
335  // Loop over energy bins
336  for (int i = 0; i < ebds.size(); ++i) {
337 
338  // Get log mean energy for bin
339  GEnergy energy = ebds.elogmean(i);
340 
341  // Get weight for this bin
342  double weight = spectrum.flux(ebds.emin(i), ebds.emax(i)) /
343  flux_total;
344 
345  // Initialise sum
346  double sum = 0.0;
347 
348  // Continue only if we have weight
349  if (weight > 0.0) {
350 
351  // Generate IAQ using Compton kinematics
352  compton_kinematics(energy.MeV());
353 
354  // Multiply IAQ with Klein-Nishina factors
355  klein_nishina(energy.MeV());
356 
357  // Weight IAQ
358  weight_iaq(energy.MeV());
359 
360  // Copy weigthed IAQ values
361  for (int k = 0; k < iaq.npix(); ++k) {
362  double value = m_iaq(k) * weight;
363  iaq(k) += value;
364  sum += value;
365  }
366 
367  } // endif: we has weight
368 
369  // Debug
370  #if defined(G_DEBUG_SET_CONTINUUM)
371  std::cout << "Energy=" << energy;
372  std::cout << " weight=" << weight;
373  std::cout << " sum=" << sum << std::endl;
374  #endif
375 
376  } // endfor: looped over energy bins
377 
378  // Store IAQ matrix
379  m_iaq = iaq;
380 
381  // Apply location smearing
383 
384  // Set parameters
385  m_iaq.card("SPECTRUM", spectrum.classname(), "Source spectrum");
386  const GModelSpectralPlaw* plaw =
387  dynamic_cast<const GModelSpectralPlaw*>(&spectrum);
388  const GModelSpectralPlawPhotonFlux* pplaw =
389  dynamic_cast<const GModelSpectralPlawPhotonFlux*>(&spectrum);
390  const GModelSpectralPlawEnergyFlux* eplaw =
391  dynamic_cast<const GModelSpectralPlawEnergyFlux*>(&spectrum);
392  if (plaw != NULL) {
393  m_iaq.card("PLAWINX", plaw->index(), "Power law spectral index");
394  }
395  else if (pplaw != NULL) {
396  m_iaq.card("PLAWINX", pplaw->index(), "Power law spectral index");
397  }
398  else if (eplaw != NULL) {
399  m_iaq.card("PLAWINX", eplaw->index(), "Power law spectral index");
400  }
401  m_iaq.card("NENG", m_num_energies, "Number of incident energies");
402  m_iaq.card("EIMIN", energy_min, "[MeV] Minimum incident energy");
403  m_iaq.card("EIMAX", energy_max, "[MeV] Maximum incident energy");
404 
405  // Return
406  return;
407 }
408 
409 
410 /***********************************************************************//**
411  * @brief Save COMPTEL instrument response representation into FITS file
412  *
413  * @param[in] filename FITS file name.
414  * @param[in] clobber Overwrite existing file?
415  *
416  * Saves the COMPTEL instrument response representation into a FITS file. If
417  * the file exists already and the @p clobber parameter is true, the method
418  * will overwrite the content of the existing file. Otherwise, an exception
419  * is thrown.
420  ***************************************************************************/
421 void GCOMIaq::save(const GFilename& filename, const bool& clobber) const
422 {
423  // Work on a copy of the image to set the extension name
424  GFitsImageFloat iaq = m_iaq;
425 
426  // Set extension name
427  if (filename.has_extname()) {
428  iaq.extname(filename.extname());
429  }
430 
431  // Set FITS image keywords
432  iaq.card("EMIN", m_ebounds.emin().MeV(), "[MeV] Minimum measured photon energy");
433  iaq.card("EMAX", m_ebounds.emax().MeV(), "[MeV] Maximum measured photon energy");
434  iaq.card("E1MIN", m_e1min, "[MeV] Minimum D1 energy deposit");
435  iaq.card("E1MAX", m_e1max, "[MeV] Maximum D1 energy deposit");
436  iaq.card("E2MIN", m_e1min, "[MeV] Minimum D2 energy deposit");
437  iaq.card("E2MAX", m_e1max, "[MeV] Maximum D2 energy deposit");
438  iaq.card("ZENITH", m_zenith, "[deg] Zenith angle for location smearing");
439  iaq.card("PSDCORR", m_psd_correct, "PSD correction for 0-110");
440  iaq.card("PHIRES", m_phibar_resolution, "[deg] Phibar resolution in computation");
441 
442  // Initialise empty FITS file
443  GFits fits;
444 
445  // Append IAQ to FITS file
446  fits.append(iaq);
447 
448  // Save FITS file (without extension name which was extracted earlier
449  // and set in the HDU)
450  fits.saveto(filename.url(), clobber);
451 
452  // Return
453  return;
454 }
455 
456 
457 /***********************************************************************//**
458  * @brief Print COMPTEL instrument response representation information
459  *
460  * @param[in] chatter Chattiness.
461  * @return String containing COMPTEL instrument response representation
462  * information.
463  ***************************************************************************/
464 std::string GCOMIaq::print(const GChatter& chatter) const
465 {
466  // Initialise result string
467  std::string result;
468 
469  // Continue only if chatter is not silent
470  if (chatter != SILENT) {
471 
472  // Append header
473  result.append("=== GCOMIaq ===");
474 
475  // Append information
476  result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
477  if (m_iaq.naxis() > 0) {
478  result.append(gammalib::str(m_iaq.naxes(0)));
479  }
480  else {
481  result.append("0");
482  }
483  result.append("\n"+gammalib::parformat("Number of Phibar bins"));
484  if (m_iaq.naxis() > 1) {
485  result.append(gammalib::str(m_iaq.naxes(1)));
486  }
487  else {
488  result.append("0");
489  }
490  result.append("\n"+gammalib::parformat("Maximum Phigeo value"));
491  result.append(gammalib::str(m_phigeo_max)+" deg");
492  result.append("\n"+gammalib::parformat("Maximum Phibar value"));
493  result.append(gammalib::str(m_phibar_max)+" deg");
494  result.append("\n"+gammalib::parformat("Phigeo bin size"));
495  result.append(gammalib::str(m_phigeo_bin_size)+" deg");
496  result.append("\n"+gammalib::parformat("Phibar bin size"));
497  result.append(gammalib::str(m_phibar_bin_size)+" deg");
498  result.append("\n"+gammalib::parformat("Energy range"));
499  result.append(m_ebounds.emin().print()+" - ");
500  result.append(m_ebounds.emax().print());
501 
502  // Append parameters
503  result.append("\n"+gammalib::parformat("D1 energy range"));
504  result.append(gammalib::str(m_e1min)+" MeV - ");
505  result.append(gammalib::str(m_e1max)+" MeV");
506  result.append("\n"+gammalib::parformat("D2 energy range"));
507  result.append(gammalib::str(m_e2min)+" MeV - ");
508  result.append(gammalib::str(m_e2max)+" MeV");
509  result.append("\n"+gammalib::parformat("Zenith angle"));
510  result.append(gammalib::str(m_zenith)+" deg");
511  result.append("\n"+gammalib::parformat("Phibar resolution"));
512  result.append(gammalib::str(m_phibar_resolution)+" deg");
513  result.append("\n"+gammalib::parformat("PSD correction"));
514  if (m_psd_correct) {
515  result.append("yes");
516  }
517  else {
518  result.append("no");
519  }
520  result.append("\n"+gammalib::parformat("Number of input energies"));
521  result.append(gammalib::str(m_num_energies));
522 
523  // Append IAQ integral
524  double sum = 0.0;
525  for (int i = 0; i < m_iaq.npix(); ++i) {
526  sum += m_iaq(i);
527  }
528  result.append("\n"+gammalib::parformat("Response integral"));
529  result.append(gammalib::str(sum));
530 
531  } // endif: chatter was not silent
532 
533  // Return result
534  return result;
535 }
536 
537 
538 /*==========================================================================
539  = =
540  = Private methods =
541  = =
542  ==========================================================================*/
543 
544 /***********************************************************************//**
545  * @brief Initialise class members
546  ***************************************************************************/
548 {
549  // Initialise members
550  m_iaq.clear();
551  m_ebounds.clear();
554  m_ict.clear();
555  m_phigeo_max = 0.0;
556  m_phibar_max = 0.0;
557  m_phigeo_bin_size = 0.0;
558  m_phibar_bin_size = 0.0;
559 
560  // Initialise parameters
561  m_phibar_resolution = 0.25; //!< Default: 0.25 deg
562  m_e1min = 0.070; //!< Default: 70 keV
563  m_e1max = 20.0; //!< Default: 20 MeV
564  m_e2min = 0.650; //!< Default: 650 keV
565  m_e2max = 30.0; //!< Default: 30 MeV
566  m_num_energies = 50; //!< Default: 50 input energies
567  m_psd_correct = true; //!< Default: use PSD correction
568  m_zenith = 25.0; //!< Default: 25 deg
569 
570  // Return
571  return;
572 }
573 
574 
575 /***********************************************************************//**
576  * @brief Copy class members
577  *
578  * @param[in] iaq COMPTEL instrument response representation.
579  ***************************************************************************/
581 {
582  // Copy attributes
583  m_iaq = iaq.m_iaq;
584  m_ebounds = iaq.m_ebounds;
587  m_ict = iaq.m_ict;
592 
593  // Copy parameters
595  m_e1min = iaq.m_e1min;
596  m_e1max = iaq.m_e1max;
597  m_e2min = iaq.m_e2min;
598  m_e2max = iaq.m_e2max;
601  m_zenith = iaq.m_zenith;
602 
603  // Return
604  return;
605 }
606 
607 
608 /***********************************************************************//**
609  * @brief Delete class members
610  ***************************************************************************/
612 {
613  // Return
614  return;
615 }
616 
617 
618 /***********************************************************************//**
619  * @brief Initialise COMPTEL response function and instrument characteristics
620  ***************************************************************************/
622 {
623  // Initialise COMPTEL response function and instrument characteristics
624  GCaldb caldb("cgro","comptel");
625  m_response_d1 = GCOMD1Response(caldb, "DEFAULT");
626  m_response_d2 = GCOMD2Response(caldb, "DEFAULT");
627  m_ict = GCOMInstChars(caldb, "DEFAULT");
628 
629  // Return
630  return;
631 }
632 
633 
634 /***********************************************************************//**
635  * @brief Remove any extra header cards
636  *
637  * Remove all header cards that may have been attached by the set() methods.
638  ***************************************************************************/
640 {
641  // Remove extra header cards
642  if (m_iaq.has_card("ENERGY")) {
643  const_cast<GFitsHeader&>(m_iaq.header()).remove("ENERGY");
644  }
645  if (m_iaq.has_card("SPECTRUM")) {
646  const_cast<GFitsHeader&>(m_iaq.header()).remove("SPECTRUM");
647  }
648  if (m_iaq.has_card("PLAWINX")) {
649  const_cast<GFitsHeader&>(m_iaq.header()).remove("PLAWINX");
650  }
651  if (m_iaq.has_card("NENG")) {
652  const_cast<GFitsHeader&>(m_iaq.header()).remove("NENG");
653  }
654  if (m_iaq.has_card("EIMIN")) {
655  const_cast<GFitsHeader&>(m_iaq.header()).remove("EIMIN");
656  }
657  if (m_iaq.has_card("EIMAX")) {
658  const_cast<GFitsHeader&>(m_iaq.header()).remove("EIMAX");
659  }
660 
661  // Return
662  return;
663 }
664 
665 
666 /***********************************************************************//**
667  * @brief Generate IAQ matrix based on Compton kinematics
668  *
669  * @param[in] energy Input photon energy (MeV).
670  *
671  * Generates an IAQ matrix based on Compton kinematics only. See the
672  * compute_iaq_bin() method for the formula used to compute each bin of the
673  * IAQ matrix.
674  *
675  * The Phibar dimension of the IAQ is sampled at a resolution that is set
676  * by the m_phibar_resolution member.
677  *
678  * The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F
679  * (release 1.0, 24-FEB-93).
680  ***************************************************************************/
681 void GCOMIaq::compton_kinematics(const double& energy)
682 {
683  // Get IAQ dimensions
684  int n_phigeo = m_iaq.naxes(0);
685  int n_phibar = m_iaq.naxes(1);
686 
687  // Set the phibar step size for the internal computation. The phibar
688  // step size is set by the m_phibar_resolution member. This ensures that
689  // the computation is done with a sufficient resolution in phibar.
690  int n_fine = int(m_phibar_bin_size / m_phibar_resolution + 0.5);
691  if (n_fine < 1) {
692  n_fine = 1;
693  }
694  double dphibar = m_phibar_bin_size / double(n_fine);
695  double dphibar_rad = dphibar * gammalib::deg2rad;
696 
697  // Loop over phigeo
698  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
699 
700  // Get geometrical scatter angle (deg)
701  double phigeo = (double(i_phigeo) + 0.5) * m_phigeo_bin_size;
702 
703  // Compute the true D1 and D2 energy deposits based on the
704  // geometrical scatter angle
705  double etrue2 = gammalib::com_energy2(energy, phigeo);
706  double etrue1 = energy - etrue2;
707 
708  // Debug
709  #if defined(G_DEBUG_COMPTON_KINEMATICS)
710  std::cout << "phigeo=" << phigeo;
711  std::cout << " etrue1=" << etrue1;
712  std::cout << " etrue2=" << etrue2;
713  #endif
714 
715  // Initialise sum
716  double sum = 0.0;
717 
718  // Loop over phibar
719  for (int i_phibar = 0; i_phibar < n_phibar; ++i_phibar) {
720 
721  // Initialise start of phibar for this layer
722  double phibar = double(i_phibar) * m_phibar_bin_size + 0.5 * dphibar;
723 
724  // Initialise response for this phibar layer
725  double response = 0.0;
726 
727  // Loop fine phibar layers
728  for (int i = 0; i < n_fine; ++i, phibar += dphibar) {
729 
730  // Compute IAQ bin
731  response += compute_iaq_bin(etrue1, etrue2, phibar);
732 
733  } // endfor: looper over fine phibar layers
734 
735  // Multiply response by size of fine phibar layers in radians
736  response *= dphibar_rad;
737 
738  // Store result
739  m_iaq(i_phigeo, i_phibar) = response;
740 
741  // Add response to sum
742  sum += response;
743 
744  } // endfor: looped over phibar
745 
746  // Debug
747  #if defined(G_DEBUG_COMPTON_KINEMATICS)
748  std::cout << " Prob=" << sum << std::endl;
749  #endif
750 
751  } // endfor: looped over phigeo
752 
753  // Return
754  return;
755 }
756 
757 
758 /***********************************************************************//**
759  * @brief Computes the IAQ for one bin
760  *
761  * @param[in] etrue1 True D1 energy deposit (MeV).
762  * @param[in] etrue2 True D2 energy deposit (MeV).
763  * @param[in] phibar Compton scatter angle (deg).
764  * @return Response for one IAQ bin
765  *
766  * Compute the integral
767  *
768  * \f[
769  * R(\bar{\varphi}|\hat{E_1},\hat{E_2}) =
770  * \int_{E_1^{\rm min}}^{E_1^{\rm max}}
771  * R(E_1,E_2|\hat{E_1},\hat{E_2}) dE_1
772  * \f]
773  *
774  * where
775  * \f$R(E_1,E_2|\hat{E_1},\hat{E_2})\f$ is computed using
776  * GCOMIaq::response_kernel::eval,
777  *
778  * \f$E_1\f$ and \f$E_2\f$ are the measured D1 and D2 energy deposits,
779  * respectively,
780  * \f$\hat{E_1}\f$ and \f$\hat{E_2}\f$ are the true D1 and D2 energy
781  * deposits, respectively, and \f$\bar{\varphi}\f$ the Compton scatter
782  * angle.
783  *
784  * The code implementation is based on the COMPASS RESPSIT2 functions
785  * RESPSC.F (release 2.0, 14-Mar-91) and RESINT.F (release 3.0, 21-Oct-91).
786  * Since RESPSC.F only defines the threshold while RESINT.F does the real
787  * job, both functions were merged.
788  *
789  * @todo Study impact of integration precision.
790  ***************************************************************************/
791 double GCOMIaq::compute_iaq_bin(const double& etrue1,
792  const double& etrue2,
793  const double& phibar)
794 {
795  // Initialise response
796  double response = 0.0;
797 
798  // Set integration limits
799  double e1min = m_response_d1.emin(etrue1) - 0.5 * m_response_d1.ewidth(etrue1);
800  double e1max = m_response_d1.emax(etrue1);
801 
802  // Set integration interval. Do not integrate more than 3 standard
803  // deviations away from photo peak.
804  double position = m_response_d1.position(etrue1);
805  double sigma = m_response_d1.sigma(etrue1);
806  double e_low = position - 3.0 * sigma;
807  double e_high = position + 3.0 * sigma;
808  double e_min = (e1min < e_low) ? e_low : e1min;
809  double e_max = (e1max > e_high) ? e_high : e1max;
810 
811  // Continue only if the integration interval is positive
812  if (e_min < e_max) {
813 
814  // Setup integration kernel
815  response_kernel integrand(m_response_d1,
817  etrue1,
818  etrue2,
819  phibar,
820  m_ebounds.emin().MeV(),
821  m_ebounds.emax().MeV(),
822  m_e1min,
823  m_e1max,
824  m_e2min,
825  m_e2max);
826 
827  // Setup integral
828  GIntegral integral(&integrand);
829 
830  // Set precision
831  integral.eps(1.0e-4);
832 
833  // No warnings
834  #if defined(G_COMPUTE_IAQ_BIN_NO_WARNINGS)
835  integral.silent(true);
836  #endif
837 
838  // Perform integration
839  response = integral.romberg(e_min, e_max);
840 
841  } // endif: integration interval was positive
842 
843  // Debug
844  #if defined(G_DEBUG_COMPUTE_IAQ_BIN)
845  std::cout << " phibar=" << phibar;
846  std::cout << " e1min=" << e1min;
847  std::cout << " e1max=" << e1max;
848  std::cout << " e_min=" << e_min;
849  std::cout << " e_max=" << e_max;
850  std::cout << " response=" << response << std::endl;
851  #endif
852 
853  // Return response
854  return response;
855 }
856 
857 
858 /***********************************************************************//**
859  * @brief Multiply IAQ matrix by the Klein-Nishina formula
860  *
861  * @param[in] energy Input photon energy (MeV).
862  *
863  * The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F
864  * (release 1.0, 24-FEB-93).
865  ***************************************************************************/
866 void GCOMIaq::klein_nishina(const double& energy)
867 {
868  // Get IAQ dimensions
869  int n_phigeo = m_iaq.naxes(0);
870  int n_phibar = m_iaq.naxes(1);
871 
872  // Debug
873  #if defined(G_DEBUG_KLEIN_NISHINA)
874  double sum = 0.0;
875  #endif
876 
877  // Loop over phigeo
878  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
879 
880  // Get geometrical scatter angle (deg)
881  double phigeo = (double(i_phigeo) + 0.5) * m_phigeo_bin_size;
882 
883  // Get Klein-Nishina value for Phigeo
884  double prob_kn = klein_nishina_bin(energy, phigeo);
885 
886  // Debug
887  #if defined(G_DEBUG_KLEIN_NISHINA)
888  std::cout << "phigeo=" << phigeo << " prob_kn=" << prob_kn << std::endl;
889  sum += prob_kn;
890  #endif
891 
892  // Loop over phibar
893  for (int i_phibar = 0; i_phibar < n_phibar; ++i_phibar) {
894 
895  // Store result
896  m_iaq(i_phigeo, i_phibar) *= prob_kn;
897 
898  } // endfor: looped over phibar
899 
900  } // endfor: looped over phigeo
901 
902  // Debug
903  #if defined(G_DEBUG_KLEIN_NISHINA)
904  std::cout << "Sum of probabilities = " << sum << std::endl;
905  #endif
906 
907  // Return
908  return;
909 }
910 
911 
912 /***********************************************************************//**
913  * @brief Computes Klein-Nishina probability for one bin
914  *
915  * @param[in] energy Input photon energy (MeV).
916  * @param[in] phigeo Geometric scatter angle (deg).
917  * @return Klein-Nishina probability.
918  *
919  * Compute the probability that a photon of @p energy is Compton scattered
920  * in a given \f$\varphi_{\rm geo}\f$ bin using
921  *
922  * \f[
923  * P = \frac{\int_{\varphi_{\rm geo,min}}^{\varphi_{\rm geo,max}}
924  * \sigma_{\rm KN}(E, \varphi_{\rm geo}) d\varphi_{\rm geo}}
925  * {\int_{0}^{\pi}
926  * \sigma_{\rm KN}(E, \varphi_{\rm geo}) d\varphi_{\rm geo}}
927  * \f]
928  *
929  * where \f$\sigma_{\rm KN}(E, \varphi_{\rm geo})\f$ is the Klein-Nishina
930  * cross section.
931  *
932  * The code implementation is based on the COMPASS RESPSIT2 function
933  * RESPSG.F (release 2.0, 26-Apr-90).
934  ***************************************************************************/
935 double GCOMIaq::klein_nishina_bin(const double& energy, const double& phigeo)
936 {
937  // Compute phigeo boundaries in radians. Make sure that they remain in
938  // the range [0,pi].
939  double phigeo_lo = phigeo - 0.5 * m_phigeo_bin_size;
940  double phigeo_hi = phigeo + 0.5 * m_phigeo_bin_size;
941  if (phigeo_lo < 0.0) {
942  phigeo_lo = 0.0;
943  }
944  if (phigeo_hi < 0.0) {
945  phigeo_hi = 0.0;
946  }
947  if (phigeo_lo > 180.0) {
948  phigeo_lo = 180.0;
949  }
950  if (phigeo_hi > 180.0) {
951  phigeo_hi = 180.0;
952  }
953  phigeo_lo *= gammalib::deg2rad;
954  phigeo_hi *= gammalib::deg2rad;
955 
956  // Express input photon energy in units of electron rest mass
957  double alpha = energy / gammalib::mec2;
958 
959  // Calculate the probability of phigeo bin
960  double v_low = 1.0 - 1.0 / (1.0 + (1.0 - std::cos(phigeo_lo)) * alpha);
961  double v_high = 1.0 - 1.0 / (1.0 + (1.0 - std::cos(phigeo_hi)) * alpha);
962  double r_low = klein_nishina_integral(v_low, alpha);
963  double r_high = klein_nishina_integral(v_high, alpha);
964  double prob_bin = (r_high - r_low);
965 
966  // Calculate the probability of [0,pi] range
967  const double v_zero = 0.0;
968  const double v_pi = 1.0 - 1.0 / (1.0 + 2.0 * alpha);
969  const double r_zero = klein_nishina_integral(v_zero, alpha);
970  const double r_pi = klein_nishina_integral(v_pi, alpha);
971  const double prob_tot = (r_pi - r_zero);
972 
973  // Calculate the probability of phigeo bin
974  double prob_phigeo = prob_bin / prob_tot;
975 
976  // Return probability of phigeo bin
977  return prob_phigeo;
978 }
979 
980 
981 /***********************************************************************//**
982  * @brief Computes Klein-Nishina integral
983  *
984  * @param[in] v Integration limit.
985  * @param[in] a Normalized energy.
986  * @return Integrated Klein-Nishina cross section.
987  *
988  * Computes the indefinite integral of the Klein-Nishina cross section
989  * evaluated at the integration limit @p v.
990  *
991  * The code implementation is based on the COMPASS RESPSIT2 function
992  * RESSLV.F (release 1.0, 23-Nov-88). The formula used in that code is
993  * the following:
994  *
995  * W = V*( 1.D0/A + 1.D0 )**2
996  * W = W + ( 2.D0/(A**2) + 2.D0/A - 1.D0 )* DLOG(1.D0 - V)
997  * W = W - (V**2)/2.D0
998  * W = W + 1.D0/(A**2) * 1.D0/(1.D0 - V)
999  *
1000  * It has been check that the original formula gives the same results as the
1001  * version that was implemented.
1002  ***************************************************************************/
1003 double GCOMIaq::klein_nishina_integral(const double& v, const double& a)
1004 {
1005  // Compute integral (hope this is okay :o)
1006  double w = v * (1.0/a + 1.0) * (1.0/a + 1.0) +
1007  (2.0/(a*a) + 2.0/a - 1.0) * std::log(1.0 - v) -
1008  v*v/2.0 +
1009  1.0/(a*a) * 1.0/(1.0 - v);
1010 
1011  // Debug
1012  #if defined(G_DEBUG_KLEIN_NISHINA_INTEGRAL)
1013  double w_test;
1014  w_test = v*( 1.0/a + 1.0 )*( 1.0/a + 1.0 );
1015  w_test = w_test + ( 2.0/(a*a) + 2.0/a - 1.0 )* std::log(1.0 - v);
1016  w_test = w_test - (v*v)/2.0;
1017  w_test = w_test + 1.0/(a*a) * 1.0/(1.0 - v);
1018  if (w != w_test) {
1019  std::cout << "w=" << w << " w_test=" << w_test << std::endl;
1020  }
1021  #endif
1022 
1023  // Return integral
1024  return w;
1025 }
1026 
1027 
1028 /***********************************************************************//**
1029  * @brief Weight IAQ matrix
1030  *
1031  * @param[in] energy Input photon energy (MeV).
1032  *
1033  * The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F
1034  * (release 1.0, 24-FEB-93).
1035  ***************************************************************************/
1036 void GCOMIaq::weight_iaq(const double& energy)
1037 {
1038  // Compute the total D1 interaction coefficient for source on axis
1039  double d1prob = m_ict.prob_D1inter(energy);
1040 
1041  // Compute the fraction of Compton interactions within the total D1
1042  // interaction probability (COM-RP-ROL-DRG-41)
1043  double cfract = 1.0;
1044  if (energy >= 2.0) {
1045  cfract = 1.067 - 0.0295 * energy + 3.4e-4 * energy * energy;
1046  }
1047  if (cfract > 1.0) {
1048  cfract = 1.0;
1049  }
1050  else if (cfract < 0.0) {
1051  cfract = 0.0;
1052  }
1053 
1054  // Compute the transmission of material above D1
1055  double ad1trans = m_ict.trans_D1(energy);
1056 
1057  // Compute the transmission of V1 veto dome
1058  //double v1trans = m_ict.trans_V1(energy);
1059 
1060  // Compute the probability that a photon was not self vetoed assuming
1061  // a source on axis
1062  double sveto = m_ict.prob_no_selfveto(energy, 0.0);
1063 
1064  // Compute the probability that an event is not a multihit
1065  double mlhit = m_ict.prob_no_multihit(energy);
1066 
1067  // Compute the overall (shape independent) transmission coefficients
1068  //double oalltr = ad1trans * v1trans * d1prob * cfract * mlhit * sveto;
1069  double oalltr = ad1trans * d1prob * cfract * mlhit * sveto;
1070 
1071  // Debug
1072  #if defined(G_DEBUG_WEIGHT_IAQ)
1073  std::cout << "Transmission coefficients:" << std::endl;
1074  std::cout << "==========================" << std::endl;
1075  std::cout << " Above D1 transmission ..........: " << ad1trans << std::endl;
1076  std::cout << " V1 veto dome transmission ......: " << v1trans << std::endl;
1077  std::cout << " D1 interaction probability .....: " << d1prob << std::endl;
1078  std::cout << " Compton scatter fraction .......: " << cfract << std::endl;
1079  std::cout << " Compton interaction probability : " << d1prob*cfract << std::endl;
1080  std::cout << " Multihit transmission ..........: " << mlhit << std::endl;
1081  std::cout << " Self-vetoing transmission ......: " << sveto << std::endl;
1082  std::cout << " Overall shape-independent trans.: " << oalltr << std::endl;
1083  #endif
1084 
1085  // Get IAQ dimensions
1086  int n_phigeo = m_iaq.naxes(0);
1087  int n_phibar = m_iaq.naxes(1);
1088 
1089  // Loop over phigeo
1090  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1091 
1092  // Get geometrical scatter angle (deg)
1093  double phigeo = (double(i_phigeo) + 0.5) * m_phigeo_bin_size;
1094 
1095  // Compute the transmission between D1 and D2
1096  double ad2trans = m_ict.trans_D2(energy, phigeo);
1097 
1098  // Compute the transmission of V2+V3 veto domes
1099  //double v23trans = m_ict.trans_V23(energy, phigeo);
1100 
1101  // Compute the D2 interaction coefficient
1102  double d2prob = m_ict.prob_D2inter(energy, phigeo);
1103 
1104  // Compute multi-scatter transmission inside D1 module
1105  double mscatt = m_ict.multi_scatter(energy, phigeo);
1106 
1107  // Optionally compute PSD correction for the default corrected PSD
1108  // channel selection 0-110
1109  #if defined(G_APPLY_PSD_CORRECTION_IN_RESPONSE_KERNEL)
1110  double psdtrn = 1.0;
1111  #else
1112  double psdtrn = (m_psd_correct) ? m_ict.psd_correction(energy, phigeo) : 1.0;
1113  #endif
1114 
1115  // Compute the overall phigeo dependent correction
1116  //double oallpg = ad2trans * v23trans * d2prob * mscatt * psdtrn;
1117  double oallpg = ad2trans * d2prob * mscatt * psdtrn;
1118 
1119  // Compute the overall correction
1120  double weight = oalltr * oallpg;
1121 
1122  // Debug
1123  #if defined(G_DEBUG_WEIGHT_IAQ)
1124  std::cout << " Phigeo .........................: " << phigeo << std::endl;
1125  std::cout << " Between D1 & D2 transmission ..: " << ad2trans << std::endl;
1126  std::cout << " V2+V3 veto dome transmission ..: " << v23trans << std::endl;
1127  std::cout << " D2 interaction probability ....: " << d2prob << std::endl;
1128  std::cout << " D1 multi-scatter transmission .: " << mscatt << std::endl;
1129  std::cout << " PSD correction ................: " << psdtrn << std::endl;
1130  std::cout << " Overall shape-dependent trans. : " << oallpg << std::endl;
1131  std::cout << " Overall transmission ..........: " << weight << std::endl;
1132  #endif
1133 
1134  // Apply weight to all phibar layers
1135  for (int i_phibar = 0; i_phibar < n_phibar; ++i_phibar) {
1136  m_iaq(i_phigeo, i_phibar) *= weight;
1137  }
1138 
1139  } // endfor: looped over phigeo
1140 
1141  // Return
1142  return;
1143 }
1144 
1145 
1146 /***********************************************************************//**
1147  * @brief Perform location smearing
1148  *
1149  * @param[in] zenith Zenith angle of source (deg).
1150  *
1151  * The code implementation is based on the COMPASS RESPSIT2 function LOCSPR.F
1152  * (release 1.0, 05-JAN-93).
1153  ***************************************************************************/
1154 void GCOMIaq::location_smearing(const double& zenith)
1155 {
1156  // Compute location spread for all phigeo angles
1157  std::vector<double> sigmas = location_spread(zenith);
1158 
1159  // Get IAQ dimensions
1160  int n_phigeo = m_iaq.naxes(0);
1161  int n_phibar = m_iaq.naxes(1);
1162 
1163  // Setup phigeo node array for interpolated
1164  GNodeArray phigeos;
1165  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1166  phigeos.append((double(i_phigeo) + 0.5) * m_phigeo_bin_size);
1167  }
1168 
1169  // Loop over phibar
1170  for (int i_phibar = 0; i_phibar < n_phibar; ++i_phibar) {
1171 
1172  // Copy phigeo vector
1173  std::vector<double> values;
1174  double sum_before = 0.0;
1175  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1176  values.push_back(m_iaq(i_phigeo, i_phibar));
1177  sum_before += m_iaq(i_phigeo, i_phibar);
1178  }
1179 
1180  // Compute convolution integral for each phigeo pixel
1181  std::vector<double> convolved_values;
1182  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1183 
1184  // Setup integration kernel
1185  smearing_kernel integrand(phigeos[i_phigeo],
1186  sigmas[i_phigeo],
1187  phigeos,
1188  values);
1189 
1190  // Setup integral
1191  GIntegral integral(&integrand);
1192 
1193  // Set precision
1194  integral.eps(1.0e-4);
1195 
1196  // No warnings
1197  #if defined(G_LOCATION_SMEARING_NO_WARNINGS)
1198  integral.silent(true);
1199  #endif
1200 
1201  // Get integration boundaries
1202  double phigeo_min = phigeos[i_phigeo] - 3.0 * sigmas[i_phigeo];
1203  double phigeo_max = phigeos[i_phigeo] + 3.0 * sigmas[i_phigeo];
1204 
1205  // Perform integration
1206  double value = integral.romberg(phigeo_min, phigeo_max);
1207 
1208  // Store result
1209  convolved_values.push_back(value);
1210 
1211  } // endfor: convolution integral
1212 
1213  // Restore phigeo vector
1214  double sum_after = 0.0;
1215  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1216  m_iaq(i_phigeo, i_phibar) = convolved_values[i_phigeo];
1217  sum_after += convolved_values[i_phigeo];
1218  }
1219 
1220  // Debug
1221  #if defined(G_DEBUG_LOCATION_SMEARING)
1222  std::cout << "phibar=" << (double(i_phibar) + 0.5) * m_phibar_bin_size;
1223  std::cout << " before=" << sum_before;
1224  std::cout << " after=" << sum_after;
1225  if (sum_before > 0.0) {
1226  std::cout << " (" << sum_after/sum_before << ")";
1227  }
1228  std::cout << std::endl;
1229  #endif
1230 
1231  } // endfor: looped over phibar
1232 
1233  // Return
1234  return;
1235 }
1236 
1237 
1238 /***********************************************************************//**
1239  * @brief Compute location spread vector
1240  *
1241  * @param[in] zenith Zenith angle of source (deg).
1242  * @return Vector of location spreads in Gaussian sigma
1243  *
1244  * The code implementation is based on the COMPASS RESPSIT2 function LOCSPR.F
1245  * (release 1.0, 05-JAN-93).
1246  ***************************************************************************/
1247 std::vector<double> GCOMIaq::location_spread(const double& zenith) const
1248 {
1249  // Set distance between the detector midplanes of D1 and D2 (cm)
1250  const double zdist = 158.0;
1251 
1252  // Set horizontal (x,y) location spread in D1 resp. D2 (cm)
1253  const double sighd1 = 2.30;
1254  const double sighd2 = 1.96;
1255 
1256  // Set vertical (z) location spread in D1 resp. D2 based on homogeneous
1257  // event distributions in z-direction (cm)
1258  const double sigvd1 = 2.45;
1259  const double sigvd2 = 2.17;
1260 
1261  // Derive constants
1262  const double zdist2 = zdist * zdist;
1263  const double sight2 = sighd1*sighd1 + sighd2*sighd2;
1264  const double sigvt2 = sigvd1*sigvd1 + sigvd2*sigvd2;
1265 
1266  // Initialise result
1267  std::vector<double> sigmas;
1268 
1269  // Compute terms
1270  double a1 = std::sin(zenith * gammalib::deg2rad);
1271  double a1sq = a1 * a1;
1272  double a2 = std::cos(zenith * gammalib::deg2rad);
1273  double a2sq = a2 * a2;
1274 
1275  // Get IAQ dimensions
1276  int n_phigeo = m_iaq.naxes(0);
1277  //int n_phibar = m_iaq.naxes(1);
1278 
1279  // Loop over phigeo
1280  for (int i_phigeo = 0; i_phigeo < n_phigeo; ++i_phigeo) {
1281 
1282  // Get geometrical scatter angle (deg)
1283  double phigeo = (double(i_phigeo) + 0.5) * m_phigeo_bin_size;
1284 
1285  // Compute spread due to (x,y) uncertainty
1286  double cphig = std::cos(phigeo * gammalib::deg2rad);
1287  double sphig = std::sin(phigeo * gammalib::deg2rad);
1288  double cphig2 = cphig*cphig;
1289  double sphig2 = sphig*sphig;
1290  double sigxy2 = (a2sq * cphig2 *
1291  (4.0 * a1sq * sphig2 + cphig2 - a1sq) +
1292  0.5 * a1sq *
1293  (a2sq * cphig2 + a1sq * sphig2 * (1.0 - 0.75 * cphig2))) *
1294  sight2 / zdist2;
1295 
1296  // Compute spread due to (z) uncertainty
1297  double sigz2 = (a2sq * a2sq * sphig2 * cphig2 +
1298  a2sq * a2sq * (0.5 - 3.0 * cphig2 * sphig2) +
1299  3.0 * a1sq * a1sq * sphig2 * cphig2 / 8.0) *
1300  sigvt2 / zdist2;
1301 
1302  // Calculate sigma in degrees
1303  double siggeo = std::sqrt(sigxy2 + sigz2) * gammalib::rad2deg;
1304 
1305  // Append sigma
1306  sigmas.push_back(siggeo);
1307 
1308  // Debug
1309  #if defined(G_DEBUG_LOCATION_SPREAD)
1310  std::cout << "phigeo=" << phigeo << " sigma=" << siggeo << std::endl;
1311  #endif
1312 
1313  // We should now create a warning if siggeo+zenith is larger than
1314  // 90.0 deg, at least the original code does that
1315 
1316  } // endfor: looped over phigeo
1317 
1318  // Return sigmas
1319  return sigmas;
1320 }
1321 
1322 
1323 /***********************************************************************//**
1324  * @brief Computes product of D1 and D2 responses at energy E1 and phibar
1325  *
1326  * @param[in] energy1 D1 energy deposit (MeV).
1327  *
1328  * The method computes
1329  *
1330  * \f[
1331  * R(E_1,E_2|\hat{E_1},\hat{E_2}) =
1332  * R_1(E_1|\hat{E_1}) \times R_2(E_2|\hat{E_2}) \times J
1333  * \f]
1334  *
1335  * where
1336  * \f$R_1(E_1|\hat{E_1})\f$ is the D1 module response function,
1337  * \f$R_2(E_2|\hat{E_2})\f$ is the D2 module response function,
1338  * \f$E_1\f$ and \f$E_2\f$ are the measured D1 and D2 energy deposits,
1339  * respectively,
1340  * \f$\hat{E_1}\f$ and \f$\hat{E_2}\f$ are the true D1 and D2 energy
1341  * deposits, respectively, and \f$J\f$ is the Jacobian.
1342  *
1343  * The measured D2 energy deposit is computed using
1344  *
1345  * \f[
1346  * E_2 = \frac{1}{2} \times \left(
1347  * \sqrt{\frac{4 m_e c^2 E_1}{1 - \cos \bar{\varphi}} + E_1^2} - E_1
1348  * \right)
1349  * \f]
1350  *
1351  * where
1352  * \f$m_e\,c^2\f$ the electron rest mass, and
1353  * \f$\bar{\varphi}\f$ the Compton scatter angle.
1354  *
1355  * The Jacobian is given by
1356  *
1357  * \f[
1358  * J = \frac{m_e c^2 \sqrt{E_1} \sin \bar{\varphi}}
1359  * {(1 - \cos \bar{\varphi})^2
1360  * \sqrt{\frac{4 m_e c^2}{1 - \cos \bar{\varphi}} + E_1}}
1361  * \f]
1362  *
1363  * The code implementation is based on the COMPASS RESPSIT2 function FUNC.F
1364  * (release 4.0, 21-Oct-91).
1365  ***************************************************************************/
1366 double GCOMIaq::response_kernel::eval(const double& energy1)
1367 {
1368  // Initialise result
1369  double value = 0.0;
1370 
1371  // Set E1 and 1-cos(phibar). If the values are too small, limit to small
1372  // positive numbers.
1373  #if defined(G_RESPONSE_KERNEL_LIMITS)
1374  double e1 = (energy1 < 1.0e-20) ? 1.0e-20 : energy1;
1375  double cp = 1.0 - m_cos_phibar;
1376  if (std::abs(cp) < 1.0e-20) {
1377  cp = 1.0e-20;
1378  }
1379  #else
1380  double e1 = energy1;
1381  double cp = 1.0 - m_cos_phibar;
1382  #endif
1383 
1384  // Compute E2 and Jacobian
1385  double term1 = 4.0 * gammalib::mec2 / cp;
1386  double e2 = 0.5 * (std::sqrt(term1 * e1 + e1*e1) - e1);
1387  double jc = gammalib::mec2 * std::sqrt(e1) * m_sin_phibar /
1388  (cp * cp * std::sqrt(term1 + e1));
1389 
1390  // Debug
1391  #if defined(G_DEBUG_RESPONSE_KERNEL)
1392  std::cout << " e1=" << e1;
1393  std::cout << " e2=" << e2;
1394  std::cout << " jc=" << jc << std::endl;
1395  #endif
1396 
1397  // Allow to break now at any time
1398  do {
1399 
1400  // Break if the D1 energy deposit is outside the energy range
1401  if ((e1 < m_e1min) || (e1 > m_e1max)) {
1402  break;
1403  }
1404 
1405  // Break if the D2 energy deposit is outside the energy range
1406  if ((e2 < m_e2min) || (e2 > m_e2max)) {
1407  break;
1408  }
1409 
1410  // Break if the total energy deposit is outside the energy range
1411  double etot = e1 + e2;
1412  if ((etot < m_etmin) || (etot > m_etmax)) {
1413  break;
1414  }
1415 
1416  // Break if the Jacobian precision is too poor
1417  #if defined(G_RESPONSE_KERNEL_LIMITS)
1418  if (jc < 1.0e-20) {
1419  break;
1420  }
1421  #endif
1422 
1423  // Get D1 response. Break if it is too small.
1424  double d1 = m_response_d1(m_etrue1, e1);
1425  #if defined(G_RESPONSE_KERNEL_LIMITS)
1426  if (d1 < 1.0e-20) {
1427  break;
1428  }
1429  #endif
1430 
1431  // Optionally apply PSD correction
1432  #if defined(G_APPLY_PSD_CORRECTION_IN_RESPONSE_KERNEL)
1433  const double a1 = 1727.9;
1434  const double a2 = 2.530;
1435  d1 *= 1.0 - (1.0 / (a1 * std::pow(e1, a2) + 1.0));
1436  #endif
1437 
1438  // Get D2 response. Break if it is too small.
1439  double d2 = m_response_d2(m_etrue2, e2);
1440  #if defined(G_RESPONSE_KERNEL_LIMITS)
1441  if (d2 < 1.0e-20) {
1442  break;
1443  }
1444  #endif
1445 
1446  // Compute value
1447  value = jc * d1 * d2;
1448 
1449  } while(false);
1450 
1451  // Return value
1452  return value;
1453 }
1454 
1455 
1456 /***********************************************************************//**
1457  * @brief Computes product of D1 and D2 responses at energy E1 and phibar
1458  *
1459  * @param[in] phigeo Geometrical scatter angle.
1460  * @return Bla ...
1461  *
1462  ***************************************************************************/
1463 double GCOMIaq::smearing_kernel::eval(const double& phigeo)
1464 {
1465  // Get interpolated IAQ value for geometrical scatter angle. Make sure
1466  // that the value is positive
1467  double value = m_phigeos.interpolate(phigeo, m_values);
1468  if (value < 0.0) {
1469  value = 0.0;
1470  }
1471 
1472  // Multiply with Gaussian
1473  double arg = (m_phigeo - phigeo) * m_wgt;
1474  value *= m_norm * std::exp(-0.5 * arg * arg);
1475 
1476  // Return value
1477  return value;
1478 }
const double mec2
Definition: GTools.hpp:52
double m_phigeo_bin_size
Bin size in geometrical scatter angle (deg)
Definition: GCOMIaq.hpp:167
void clear(void)
Clear instance.
void compton_kinematics(const double &energy)
Generate IAQ matrix based on Compton kinematics.
Definition: GCOMIaq.cpp:681
Photon flux normalized power law spectral model class.
const int & naxis(void) const
Return dimension of image.
Definition: GFitsImage.hpp:156
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
Node array class.
Definition: GNodeArray.hpp:60
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
double prob_no_multihit(const double &energy) const
Return probability that no multihit occured.
double m_phibar_bin_size
Bin size in Compton scatter angle (deg)
Definition: GCOMIaq.hpp:168
Abstract spectral model base class.
double eval(const double &phigeo)
Computes product of D1 and D2 responses at energy E1 and phibar.
Definition: GCOMIaq.cpp:1463
Interface for the COMPTEL instrument response representation class.
Definition: GCOMIaq.hpp:52
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.
double emax(const double &etrue) const
Return maximum energy.
void init_response(void)
Initialise COMPTEL response function and instrument characteristics.
Definition: GCOMIaq.cpp:621
double m_e1min
Minimum D1 energy (MeV)
Definition: GCOMIaq.hpp:132
double trans_D1(const double &energy) const
Return transmission above D1.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
double m_etrue1
True D1 energy (MeV)
Definition: GCOMIaq.hpp:126
double prob_D1inter(const double &energy) const
Return D1 interaction probability.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
void save(const GFilename &filename, const bool &clobber=false) const
Save COMPTEL instrument response representation into FITS file.
Definition: GCOMIaq.cpp:421
double index(void) const
Return power law index.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double com_energy2(const double &energy, const double &phigeo)
Return D2 energy deposit.
double m_e1max
Maximum D1 energy (MeV)
Definition: GCOMIaq.hpp:133
bool has_extname(void) const
Signal if filename has an extension name.
Definition: GFilename.hpp:255
double m_phigeo_max
Maximum geometrical scatter angle (deg)
Definition: GCOMIaq.hpp:165
void set(const GEnergy &energy, const GEbounds &ebounds)
Set mono-energetic IAQ.
Definition: GCOMIaq.cpp:253
double m_e1max
Maximum D1 energy (MeV)
Definition: GCOMIaq.hpp:171
double m_sin_phibar
sin(phibar)
Definition: GCOMIaq.hpp:129
GCOMIaq(void)
Void constructor.
Definition: GCOMIaq.cpp:78
std::vector< double > location_spread(const double &zenith) const
Compute location spread vector.
Definition: GCOMIaq.cpp:1247
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
Interface for the COMPTEL Instrument Characteristics class.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
double eval(const double &energy1)
Computes product of D1 and D2 responses at energy E1 and phibar.
Definition: GCOMIaq.cpp:1366
Gammalib tools definition.
double m_phibar_resolution
Bin size for oversampling (deg)
Definition: GCOMIaq.hpp:169
GIntegral class interface definition.
Definition: GIntegral.hpp:45
FITS file class.
Definition: GFits.hpp:63
double m_phibar_max
Maximum Compton scatter angle (deg)
Definition: GCOMIaq.hpp:166
const GFitsHeader & header(void) const
Return extension header.
Definition: GFitsHDU.hpp:205
FITS file class interface definition.
void klein_nishina(const double &energy)
Multiply IAQ matrix by the Klein-Nishina formula.
Definition: GCOMIaq.cpp:866
Energy flux normalized power law spectral model class interface definition.
Power law spectral model class.
virtual std::string classname(void) const =0
Return class name.
Implementation of support function used by COMPTEL classes.
COMPTEL instrument response representation class interface definition.
void location_smearing(const double &zenith)
Perform location smearing.
Definition: GCOMIaq.cpp:1154
double m_e2max
Maximum D2 energy (MeV)
Definition: GCOMIaq.hpp:135
double psd_correction(const double &energy, const double &phigeo) const
Return PSD correction.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
double m_etmin
Minimum total energy (MeV)
Definition: GCOMIaq.hpp:130
double m_e2max
Maximum D2 energy (MeV)
Definition: GCOMIaq.hpp:173
void clear(void)
Clear instance.
Interface for FITS header class.
Definition: GFitsHeader.hpp:49
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
Energy boundaries container class.
Definition: GEbounds.hpp:60
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
Interface for the COMPTEL D2 module response class.
Interface for the COMPTEL D1 module response class.
double compute_iaq_bin(const double &etrue1, const double &etrue2, const double &phibar)
Computes the IAQ for one bin.
Definition: GCOMIaq.cpp:791
Filename class.
Definition: GFilename.hpp:62
double klein_nishina_bin(const double &energy, const double &phigeo)
Computes Klein-Nishina probability for one bin.
Definition: GCOMIaq.cpp:935
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
GEbounds m_ebounds
Energy boundaries.
Definition: GCOMIaq.hpp:161
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double sigma(const double &etrue) const
Return photo peak standard deviation.
double ewidth(const double &etrue) const
Return energy threshold width.
GCOMD1Response m_response_d1
D1 module response.
Definition: GCOMIaq.hpp:162
bool m_psd_correct
PSD correction usage flag.
Definition: GCOMIaq.hpp:175
Power law spectral model class interface definition.
GChatter
Definition: GTypemaps.hpp:33
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1140
int m_num_energies
Number of energies for continuum IAQ.
Definition: GCOMIaq.hpp:174
GCOMD2Response m_response_d2
D2 module response.
Definition: GCOMIaq.hpp:163
double m_etrue2
True D2 energy (MeV)
Definition: GCOMIaq.hpp:127
double m_e1min
Minimum D1 energy (MeV)
Definition: GCOMIaq.hpp:170
Single precision FITS image class.
GCOMIaq & operator=(const GCOMIaq &iaq)
Assignment operator.
Definition: GCOMIaq.cpp:183
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
Flux normalized power law spectral model class interface definition.
GCOMIaq * clone(void) const
Clone instance.
Definition: GCOMIaq.cpp:236
double m_cos_phibar
cos(phibar)
Definition: GCOMIaq.hpp:128
void silent(const bool &silent)
Set silence flag.
Definition: GIntegral.hpp:243
Abstract spectral model base class interface definition.
int naxes(const int &axis) const
Return dimension of an image axis.
Definition: GFitsImage.cpp:344
GCOMInstChars m_ict
Instrument characteristics.
Definition: GCOMIaq.hpp:164
double prob_no_selfveto(const double &energy, const double &zenith) const
Return probability that photon was not self vetoed.
double m_e2min
Minimum D2 energy (MeV)
Definition: GCOMIaq.hpp:134
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:206
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
const int & npix(void) const
Return size of pixel array.
Definition: GFitsImage.hpp:136
double m_zenith
Zenith angle for location smearing (deg)
Definition: GCOMIaq.hpp:176
double position(const double &etrue) const
Return photo peak position.
double m_e2min
Minimum D2 energy (MeV)
Definition: GCOMIaq.hpp:172
void free_members(void)
Delete class members.
Definition: GCOMIaq.cpp:611
const GCOMD1Response & m_response_d1
Reference to D1 module response.
Definition: GCOMIaq.hpp:124
double index(void) const
Return power law index.
void weight_iaq(const double &energy)
Weight IAQ matrix.
Definition: GCOMIaq.cpp:1036
GFitsImageFloat m_iaq
Response.
Definition: GCOMIaq.hpp:160
double emax(const double &etrue) const
Return maximum energy.
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL instrument response representation information.
Definition: GCOMIaq.cpp:464
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
Energy boundaries class interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:665
void init_members(void)
Initialise class members.
Definition: GCOMIaq.cpp:547
double index(void) const
Return power law index.
void clear(void)
Clear instance.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
void remove_cards(void)
Remove any extra header cards.
Definition: GCOMIaq.cpp:639
void clear(void)
Clear instance.
double emin(const double &etrue) const
Return minimum energy.
~GCOMIaq(void)
Destructor.
Definition: GCOMIaq.cpp:161
const GCOMD2Response & m_response_d2
Reference to D2 module response.
Definition: GCOMIaq.hpp:125
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
virtual double flux(const GEnergy &emin, const GEnergy &emax) const =0
const double rad2deg
Definition: GMath.hpp:44
double m_etmax
Maximum total energy (MeV)
Definition: GCOMIaq.hpp:131
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void copy_members(const GCOMIaq &iaq)
Copy class members.
Definition: GCOMIaq.cpp:580
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
Integration class interface definition.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Energy flux normalized power law spectral model class.
void clear(void)
Clear instance.
Definition: GCOMIaq.cpp:217
double klein_nishina_integral(const double &v, const double &a)
Computes Klein-Nishina integral.
Definition: GCOMIaq.cpp:1003
Filename class interface definition.
double multi_scatter(const double &energy, const double &phigeo) const
Return multi-scatter correction.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413