Simon Bonnefoy / gammalib

Commit 0e4dbb59315dd048b1cfbcc4da89421d7eae94b9

Authored by Jürgen Knödlseder
1 parent e3e3cb87

Implement GCTAEdisp2D::prob_erecobin() method (#2422)

Inline Side-by-side
Showing 14 changed files with 467 additions and 58 deletions
inst/cta/include/GCTAEdisp.hpp
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdisp.hpp - Abstract CTA energy dispersion base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,6 +35,8 @@
35 35  
36 36 /* __ Forward declarations _______________________________________________ */
37 37 class GRan;
  38 +class GEnergy;
  39 +
38 40  
39 41 /***********************************************************************//**
40 42 * @class GCTAEdisp
... ... @@ -85,6 +87,10 @@ public:
85 87 const double& phi = 0.0,
86 88 const double& zenith = 0.0,
87 89 const double& azimuth = 0.0) const = 0;
  90 + virtual double prob_erecobin(const GEnergy& ereco_min,
  91 + const GEnergy& ereco_max,
  92 + const GEnergy& etrue,
  93 + const double& theta) const = 0;
88 94 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
89 95  
90 96 protected:
... ...
inst/cta/include/GCTAEdisp2D.hpp
  View file @0e4dbb5
... ... @@ -93,6 +93,10 @@ public:
93 93 const double& phi = 0.0,
94 94 const double& zenith = 0.0,
95 95 const double& azimuth = 0.0) const;
  96 + double prob_erecobin(const GEnergy& ereco_min,
  97 + const GEnergy& ereco_max,
  98 + const GEnergy& etrue,
  99 + const double& theta) const;
96 100 std::string print(const GChatter& chatter = NORMAL) const;
97 101  
98 102 // Methods
... ... @@ -106,24 +110,46 @@ public:
106 110  
107 111 private:
108 112 // Methods
109   - void init_members(void);
110   - void copy_members(const GCTAEdisp2D& edisp);
111   - void free_members(void);
112   - void update(const double& logEobs,
113   - const double& logEsrc,
114   - const double& theta) const;
115   - void compute_ebounds_obs(const double& theta = 0.0,
116   - const double& phi = 0.0,
117   - const double& zenith = 0.0,
118   - const double& azimuth = 0.0) const;
119   - void compute_ebounds_src(const double& theta = 0.0,
120   - const double& phi = 0.0,
121   - const double& zenith = 0.0,
122   - const double& azimuth = 0.0) const;
123   - void set_table(void);
124   - void set_boundaries(void);
125   - void set_max_edisp(void) const;
126   - void normalize_table(void);
  113 + void init_members(void);
  114 + void copy_members(const GCTAEdisp2D& edisp);
  115 + void free_members(void);
  116 + void update(const double& logEobs,
  117 + const double& logEsrc,
  118 + const double& theta) const;
  119 + void compute_ebounds_obs(const double& theta = 0.0,
  120 + const double& phi = 0.0,
  121 + const double& zenith = 0.0,
  122 + const double& azimuth = 0.0) const;
  123 + void compute_ebounds_src(const double& theta = 0.0,
  124 + const double& phi = 0.0,
  125 + const double& zenith = 0.0,
  126 + const double& azimuth = 0.0) const;
  127 + void set_table(void);
  128 + void set_boundaries(void);
  129 + void set_max_edisp(void) const;
  130 + void normalize_table(void);
  131 + int table_index(const int& ietrue,
  132 + const int& imigra,
  133 + const int& itheta) const;
  134 + int table_stride(const int& axis) const;
  135 + double table_value(const int& base_ll,
  136 + const int& base_lr,
  137 + const int& base_rl,
  138 + const int& base_rr,
  139 + const double& wgt_el,
  140 + const double& wgt_er,
  141 + const double& wgt_tl,
  142 + const double& wgt_tr,
  143 + const int& offset) const;
  144 + double table_value(const int& base_ll,
  145 + const int& base_lr,
  146 + const int& base_rl,
  147 + const int& base_rr,
  148 + const double& wgt_el,
  149 + const double& wgt_er,
  150 + const double& wgt_tl,
  151 + const double& wgt_tr,
  152 + const double& migra) const;
127 153  
128 154 // Kludge
129 155 void smooth_table(void);
... ...
inst/cta/include/GCTAEdispPerfTable.hpp
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispPerfTable.hpp - CTA performance table energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -88,6 +88,10 @@ public:
88 88 const double& phi = 0.0,
89 89 const double& zenith = 0.0,
90 90 const double& azimuth = 0.0) const;
  91 + double prob_erecobin(const GEnergy& ereco_min,
  92 + const GEnergy& ereco_max,
  93 + const GEnergy& etrue,
  94 + const double& theta) const;
91 95 std::string print(const GChatter& chatter = NORMAL) const;
92 96  
93 97 private:
... ...
inst/cta/include/GCTAEdispRmf.hpp
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispRmf.hpp - CTA RMF energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -87,6 +87,10 @@ public:
87 87 const double& phi = 0.0,
88 88 const double& zenith = 0.0,
89 89 const double& azimuth = 0.0) const;
  90 + double prob_erecobin(const GEnergy& ereco_min,
  91 + const GEnergy& ereco_max,
  92 + const GEnergy& etrue,
  93 + const double& theta) const;
90 94 std::string print(const GChatter& chatter = NORMAL) const;
91 95  
92 96 // Other methods
... ...
inst/cta/pyext/GCTAEdisp.i
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdisp.i - Abstract CTA energy dispersion base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -72,6 +72,10 @@ public:
72 72 const double& phi = 0.0,
73 73 const double& zenith = 0.0,
74 74 const double& azimuth = 0.0) const = 0;
  75 + virtual double prob_erecobin(const GEnergy& ereco_min,
  76 + const GEnergy& ereco_max,
  77 + const GEnergy& etrue,
  78 + const double& theta) const = 0;
75 79 };
76 80  
77 81  
... ...
inst/cta/pyext/GCTAEdisp2D.i
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdisp2D.i - CTA 2D energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015-2016 by Florent Forest *
  4 + * copyright (C) 2015-2018 by Florent Forest *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -73,8 +73,13 @@ public:
73 73 const double& phi = 0.0,
74 74 const double& zenith = 0.0,
75 75 const double& azimuth = 0.0) const;
  76 + double prob_erecobin(const GEnergy& ereco_min,
  77 + const GEnergy& ereco_max,
  78 + const GEnergy& etrue,
  79 + const double& theta) const;
76 80  
77 81 // Methods
  82 + void fetch(void) const;
78 83 const GCTAResponseTable& table(void) const;
79 84 void table(const GCTAResponseTable& table);
80 85 void read(const GFitsTable& table);
... ...
inst/cta/pyext/GCTAEdispPerfTable.i
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispPerfTable.i - CTA performance table energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -73,6 +73,10 @@ public:
73 73 const double& phi = 0.0,
74 74 const double& zenith = 0.0,
75 75 const double& azimuth = 0.0) const;
  76 + double prob_erecobin(const GEnergy& ereco_min,
  77 + const GEnergy& ereco_max,
  78 + const GEnergy& etrue,
  79 + const double& theta) const;
76 80 };
77 81  
78 82  
... ...
inst/cta/pyext/GCTAEdispRmf.i
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispRmf.i - CTA RMF energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -73,6 +73,14 @@ public:
73 73 const double& phi = 0.0,
74 74 const double& zenith = 0.0,
75 75 const double& azimuth = 0.0) const;
  76 + double prob_erecobin(const GEnergy& ereco_min,
  77 + const GEnergy& ereco_max,
  78 + const GEnergy& etrue,
  79 + const double& theta) const;
  80 +
  81 + // Other methods
  82 + int size(void) const;
  83 + const GRmf& rmf(void) const;
76 84 };
77 85  
78 86  
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @0e4dbb5
... ... @@ -779,6 +779,146 @@ void GCTAEdisp2D::fetch(void) const
779 779  
780 780  
781 781 /***********************************************************************//**
  782 + * @brief Integrate energy dispersion probability over reconstructed energy
  783 + * interval
  784 + *
  785 + * @param[in] ereco_min Minimum of reconstructed energy interval.
  786 + * @param[in] ereco_max Maximum of reconstructed energy interval.
  787 + * @param[in] etrue True energy.
  788 + * @param[in] theta Offset angle (rad).
  789 + * @return Integrated energy dispersion probability.
  790 + *
  791 + * Integrates the energy dispersion probability over an interval in
  792 + * reconstructed energy, defined by @p ereco_min and @p ereco_max, for a
  793 + * given true energy @p etrue and offset angle @p theta.
  794 + *
  795 + * The method takes into account that the energy dispersion data are stored
  796 + * in a matrix and uses the trapezoidal rule for integration of the tabulated
  797 + * data. This assures the fastest possible integration.
  798 + ***************************************************************************/
  799 +double GCTAEdisp2D::prob_erecobin(const GEnergy& ereco_min,
  800 + const GEnergy& ereco_max,
  801 + const GEnergy& etrue,
  802 + const double& theta) const
  803 +{
  804 + // Make sure that energy dispersion is online
  805 + fetch();
  806 +
  807 + // Initalize probability
  808 + double prob = 0.0;
  809 +
  810 + // Get log10 of true energy
  811 + double logEsrc = etrue.log10TeV();
  812 +
  813 + // Get migration limits
  814 + double migra_min = ereco_min / etrue;
  815 + double migra_max = ereco_max / etrue;
  816 +
  817 + // Continue only if logEsrc and theta are in validity range and migration
  818 + // interval overlaps with response table
  819 + if ((logEsrc >= m_logEsrc_min) && (logEsrc <= m_logEsrc_max) &&
  820 + (theta >= m_theta_min) && (theta <= m_theta_max) &&
  821 + (migra_max > m_migra_min) && (migra_min < m_migra_max)) {
  822 +
  823 + // Constrain migration limits to response table
  824 + if (migra_min < m_migra_min) {
  825 + migra_min = m_migra_min;
  826 + }
  827 + if (migra_max > m_migra_max) {
  828 + migra_max = m_migra_max;
  829 + }
  830 +
  831 + // Retrieve references to node arrays
  832 + const GNodeArray& etrue_nodes = m_edisp.axis_nodes(m_inx_etrue);
  833 + const GNodeArray& theta_nodes = m_edisp.axis_nodes(m_inx_theta);
  834 + const GNodeArray& migra_nodes = m_edisp.axis_nodes(m_inx_migra);
  835 +
  836 + // Set logEsrc and theta for node array interpolation
  837 + etrue_nodes.set_value(logEsrc);
  838 + theta_nodes.set_value(theta);
  839 +
  840 + // Compute base indices
  841 + int base_ll = table_index(etrue_nodes.inx_left(), 0, theta_nodes.inx_left());
  842 + int base_lr = table_index(etrue_nodes.inx_left(), 0, theta_nodes.inx_right());
  843 + int base_rl = table_index(etrue_nodes.inx_right(), 0, theta_nodes.inx_left());
  844 + int base_rr = table_index(etrue_nodes.inx_right(), 0, theta_nodes.inx_right());
  845 +
  846 + // Get migration stride
  847 + int stride = table_stride(m_inx_migra);
  848 +
  849 + // Initialise first and second node
  850 + double x1 = migra_min;
  851 + double f1 = table_value(base_ll, base_lr, base_rl, base_rr,
  852 + etrue_nodes.wgt_left(),
  853 + etrue_nodes.wgt_right(),
  854 + theta_nodes.wgt_left(),
  855 + theta_nodes.wgt_right(), x1);
  856 + double x2 = 0.0;
  857 + double f2 = 0.0;
  858 +
  859 + // Loop over all migration nodes
  860 + for (int i = 0, offset = 0; i < migra_nodes.size(); ++i, offset += stride) {
  861 +
  862 + // If migration value is below migra_min then skip node
  863 + if (migra_nodes[i] <= migra_min) {
  864 + continue;
  865 + }
  866 +
  867 + // If migration value is above maximum migration value then use
  868 + // migra_max as migration value
  869 + if (migra_nodes[i] > migra_max) {
  870 + x2 = migra_max;
  871 + f2 = table_value(base_ll, base_lr, base_rl, base_rr,
  872 + etrue_nodes.wgt_left(),
  873 + etrue_nodes.wgt_right(),
  874 + theta_nodes.wgt_left(),
  875 + theta_nodes.wgt_right(), x2);
  876 + }
  877 +
  878 + // ... otherwise use migration value
  879 + else {
  880 + x2 = migra_nodes[i];
  881 + f2 = table_value(base_ll, base_lr, base_rl, base_rr,
  882 + etrue_nodes.wgt_left(),
  883 + etrue_nodes.wgt_right(),
  884 + theta_nodes.wgt_left(),
  885 + theta_nodes.wgt_right(), offset);
  886 + }
  887 +
  888 + // Compute integral
  889 + prob += 0.5 * (f1 + f2) * std::log10(x2/x1);
  890 +
  891 + // Set second node as first node
  892 + x1 = x2;
  893 + f1 = f2;
  894 +
  895 + // If node energy is above migra_max then break now
  896 + if (migra_nodes[i] > migra_max) {
  897 + break;
  898 + }
  899 +
  900 + } // endfor: looped over all nodes
  901 +
  902 + // If last node energy is below migra_max then compute last part of
  903 + // integral up to emax
  904 + if (x1 < migra_max) {
  905 + x2 = migra_max;
  906 + f2 = table_value(base_ll, base_lr, base_rl, base_rr,
  907 + etrue_nodes.wgt_left(),
  908 + etrue_nodes.wgt_right(),
  909 + theta_nodes.wgt_left(),
  910 + theta_nodes.wgt_right(), x2);
  911 + prob += 0.5 * (f1 + f2) * std::log10(x2/x1);
  912 + }
  913 +
  914 + } // endif: logEsrc, theta and migration range were valid
  915 +
  916 + // Return probability
  917 + return prob;
  918 +}
  919 +
  920 +
  921 +/***********************************************************************//**
782 922 * @brief Print energy dispersion information
783 923 *
784 924 * @param[in] chatter Chattiness.
... ... @@ -1345,8 +1485,9 @@ void GCTAEdisp2D::normalize_table(void)
1345 1485 for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue, ++inx) {
1346 1486 double sum = sums[inx];
1347 1487 if (sum > 0.0) {
1348   - int offset = i_etrue + (i_theta*migra_size) * etrue_size;
1349   - for (int k = 0, i = offset; k < migra_size; ++k, i += etrue_size) {
  1488 + int offset = table_index(i_etrue, 0, i_theta);
  1489 + int stride = table_stride(m_inx_migra);
  1490 + for (int k = 0, i = offset; k < migra_size; ++k, i += stride) {
1350 1491 m_edisp(m_inx_matrix,i) /= sum;
1351 1492 }
1352 1493 }
... ... @@ -1360,15 +1501,160 @@ void GCTAEdisp2D::normalize_table(void)
1360 1501 }
1361 1502  
1362 1503  
  1504 +/***********************************************************************//**
  1505 + * @brief Return index of response table element
  1506 + *
  1507 + * @param[in] ietrue True energy index.
  1508 + * @param[in] imigra Migration index.
  1509 + * @param[in] itheta Offset index.
  1510 + * @return Index of response table element.
  1511 + ***************************************************************************/
  1512 +int GCTAEdisp2D::table_index(const int& ietrue,
  1513 + const int& imigra,
  1514 + const int& itheta) const
  1515 +{
  1516 + // Set index vector
  1517 + int inx[3];
  1518 + inx[m_inx_etrue] = ietrue;
  1519 + inx[m_inx_migra] = imigra;
  1520 + inx[m_inx_theta] = itheta;
  1521 +
  1522 + // Compute index
  1523 + int index = inx[0] + (inx[1] + inx[2] * m_edisp.axis_bins(1)) *
  1524 + m_edisp.axis_bins(0);
  1525 +
  1526 + // Return index
  1527 + return index;
  1528 +}
  1529 +
  1530 +
  1531 +/***********************************************************************//**
  1532 + * @brief Return stride of response table axis
  1533 + *
  1534 + * @param[in] axis Response table axis.
  1535 + * @return Stride of response table axis.
  1536 + ***************************************************************************/
  1537 +int GCTAEdisp2D::table_stride(const int& axis) const
  1538 +{
  1539 + // Initialise stride
  1540 + int stride = 1;
  1541 +
  1542 + // Multiply in higher
  1543 + if (axis > 0) {
  1544 + stride *= m_edisp.axis_bins(0);
  1545 + }
  1546 + if (axis > 1) {
  1547 + stride *= m_edisp.axis_bins(1);
  1548 + }
  1549 +
  1550 + // Return stride
  1551 + return stride;
  1552 +}
  1553 +
  1554 +
  1555 +/***********************************************************************//**
  1556 + * @brief Return bi-linearly interpolate table value for given migration bin
  1557 + *
  1558 + * @param[in] base_ll Base index for left true energy and left offset angle.
  1559 + * @param[in] base_lr Base index for left true energy and right offset angle.
  1560 + * @param[in] base_rl Base index for right true energy and left offset angle.
  1561 + * @param[in] base_rr Base index for right true energy and right offset angle.
  1562 + * @param[in] wgt_el Weighting for left true energy.
  1563 + * @param[in] wgt_er Weighting for right true energy.
  1564 + * @param[in] wgt_tl Weighting for left offset angle.
  1565 + * @param[in] wgt_tr Weighting for right offset angle.
  1566 + * @param[in] offset Offset of migration bin with respect to base indices.
  1567 + * @return Bi-linearly interpolate table value.
  1568 + ***************************************************************************/
  1569 +double GCTAEdisp2D::table_value(const int& base_ll,
  1570 + const int& base_lr,
  1571 + const int& base_rl,
  1572 + const int& base_rr,
  1573 + const double& wgt_el,
  1574 + const double& wgt_er,
  1575 + const double& wgt_tl,
  1576 + const double& wgt_tr,
  1577 + const int& offset) const
  1578 +{
  1579 + // Compute table element indices
  1580 + int inx_ll = base_ll + offset;
  1581 + int inx_lr = base_lr + offset;
  1582 + int inx_rl = base_rl + offset;
  1583 + int inx_rr = base_rr + offset;
  1584 +
  1585 + // Get
  1586 + double value = wgt_el * (m_edisp(m_inx_matrix, inx_ll) * wgt_tl +
  1587 + m_edisp(m_inx_matrix, inx_lr) * wgt_tr) +
  1588 + wgt_er * (m_edisp(m_inx_matrix, inx_rl) * wgt_tl +
  1589 + m_edisp(m_inx_matrix, inx_rr) * wgt_tr);
  1590 +
  1591 + // Return value
  1592 + return value;
  1593 +}
  1594 +
  1595 +
  1596 +/***********************************************************************//**
  1597 + * @brief Return bi-linearly interpolate table value for given migration value
  1598 + *
  1599 + * @param[in] base_ll Base index for left true energy and left offset angle.
  1600 + * @param[in] base_lr Base index for left true energy and right offset angle.
  1601 + * @param[in] base_rl Base index for right true energy and left offset angle.
  1602 + * @param[in] base_rr Base index for right true energy and right offset angle.
  1603 + * @param[in] wgt_el Weighting for left true energy.
  1604 + * @param[in] wgt_er Weighting for right true energy.
  1605 + * @param[in] wgt_tl Weighting for left offset angle.
  1606 + * @param[in] wgt_tr Weighting for right offset angle.
  1607 + * @param[in] migra Migration value.
  1608 + * @return Bi-linearly interpolate table value.
  1609 + ***************************************************************************/
  1610 +double GCTAEdisp2D::table_value(const int& base_ll,
  1611 + const int& base_lr,
  1612 + const int& base_rl,
  1613 + const int& base_rr,
  1614 + const double& wgt_el,
  1615 + const double& wgt_er,
  1616 + const double& wgt_tl,
  1617 + const double& wgt_tr,
  1618 + const double& migra) const
  1619 +{
  1620 + // Retrieve references to migration node array
  1621 + const GNodeArray& migra_nodes = m_edisp.axis_nodes(m_inx_migra);
  1622 +
  1623 + // Set migration value for node array interpolation
  1624 + migra_nodes.set_value(migra);
  1625 +
  1626 + // Get migration stride
  1627 + int stride = table_stride(m_inx_migra);
  1628 +
  1629 + // Get left value
  1630 + double value_left = table_value(base_ll, base_lr, base_rl, base_rr,
  1631 + wgt_el, wgt_er, wgt_tl, wgt_tr,
  1632 + migra_nodes.inx_left() * stride);
  1633 +
  1634 + // Get right value
  1635 + double value_right = table_value(base_ll, base_lr, base_rl, base_rr,
  1636 + wgt_el, wgt_er, wgt_tl, wgt_tr,
  1637 + migra_nodes.inx_right() * stride);
  1638 +
  1639 + // Interpolate result
  1640 + double value = migra_nodes.wgt_left() * value_left +
  1641 + migra_nodes.wgt_right() * value_right;
  1642 +
  1643 + // Make sure that interpolation result is not negative
  1644 + if (value < 0.0) {
  1645 + value = 0.0;
  1646 + }
  1647 +
  1648 + // Return
  1649 + return value;
  1650 +}
  1651 +
1363 1652  
1364 1653 /***********************************************************************//**
1365 1654 * @brief Smoothed energy dispersion table
1366 1655 ***************************************************************************/
1367 1656 void GCTAEdisp2D::smooth_table(void)
1368 1657 {
1369   - // Log entrance
1370   - //std::cout << "GCTAEdisp2D::smooth_table in" << std::endl;
1371   -
1372 1658 // Get axes dimensions
1373 1659 int netrue = m_edisp.axis_bins(m_inx_etrue);
1374 1660 int nmigra = m_edisp.axis_bins(m_inx_migra);
... ... @@ -1432,9 +1718,6 @@ void GCTAEdisp2D::smooth_table(void)
1432 1718  
1433 1719 } // endfor: looped over offset angles
1434 1720  
1435   - // Log exit
1436   - //std::cout << "GCTAEdisp2D::smooth_table out" << std::endl;
1437   -
1438 1721 // Return
1439 1722 return;
1440 1723 }
... ...
inst/cta/src/GCTAEdispPerfTable.cpp
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispPerfTable.cpp - CTA performance table energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -415,9 +415,30 @@ GEbounds GCTAEdispPerfTable::ebounds_src(const double&amp; logEobs,
415 415  
416 416  
417 417 /***********************************************************************//**
  418 + * @brief Integrated energy dispersion probability over reconstructed energy
  419 + *
  420 + * @param[in] ereco_min Minimum reconstructed energy.
  421 + * @param[in] ereco_max Maximum reconstructed energy.
  422 + * @param[in] etrue True energy.
  423 + * @param[in] theta Offset angle (degrees).
  424 + * @return Integrated energy dispersion probability.
  425 + ***************************************************************************/
  426 +double GCTAEdispPerfTable::prob_erecobin(const GEnergy& ereco_min,
  427 + const GEnergy& ereco_max,
  428 + const GEnergy& etrue,
  429 + const double& theta) const
  430 +{
  431 + // TODO
  432 +
  433 + // Return
  434 + return 0.0;
  435 +}
  436 +
  437 +
  438 +/***********************************************************************//**
418 439 * @brief Print energy dispersion information
419 440 *
420   - * @param[in] chatter Chattiness (defaults to NORMAL).
  441 + * @param[in] chatter Chattiness.
421 442 * @return String containing energy dispersion information.
422 443 ***************************************************************************/
423 444 std::string GCTAEdispPerfTable::print(const GChatter& chatter) const
... ...
inst/cta/src/GCTAEdispRmf.cpp
  View file @0e4dbb5
1 1 /***************************************************************************
2 2 * GCTAEdispRmf.cpp - CTA RMF energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2018 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -427,9 +427,30 @@ GEbounds GCTAEdispRmf::ebounds_src(const double&amp; logEobs,
427 427  
428 428  
429 429 /***********************************************************************//**
  430 + * @brief Integrated energy dispersion probability over reconstructed energy
  431 + *
  432 + * @param[in] ereco_min Minimum reconstructed energy.
  433 + * @param[in] ereco_max Maximum reconstructed energy.
  434 + * @param[in] etrue True energy.
  435 + * @param[in] theta Offset angle (degrees).
  436 + * @return Integrated energy dispersion probability.
  437 + ***************************************************************************/
  438 +double GCTAEdispRmf::prob_erecobin(const GEnergy& ereco_min,
  439 + const GEnergy& ereco_max,
  440 + const GEnergy& etrue,
  441 + const double& theta) const
  442 +{
  443 + // TODO
  444 +
  445 + // Return
  446 + return 0.0;
  447 +}
  448 +
  449 +
  450 +/***********************************************************************//**
430 451 * @brief Print RMF information
431 452 *
432   - * @param[in] chatter Chattiness (defaults to NORMAL).
  453 + * @param[in] chatter Chattiness.
433 454 * @return String containing energy dispersion information.
434 455 ***************************************************************************/
435 456 std::string GCTAEdispRmf::print(const GChatter& chatter) const
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @0e4dbb5
... ... @@ -96,6 +96,7 @@ const double minerr = 1.0e-100; //!&lt; Minimum statistical error
96 96 /* __ Macros _____________________________________________________________ */
97 97  
98 98 /* __ Coding definitions _________________________________________________ */
  99 +//#define G_RMF_INTEGRATION //!< Use numerical integration for RMF computation
99 100  
100 101 /* __ Debug definitions __________________________________________________ */
101 102 //#define G_LIKELIHOOD_DEBUG //!< Debug likelihood computation
... ... @@ -1428,6 +1429,7 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
1428 1429 logEtrue);
1429 1430  
1430 1431 // Setup energy dispersion integral
  1432 + #if defined(G_RMF_INTEGRATION)
1431 1433 GCTAOnOffObservation::edisp_kern integrand(response,
1432 1434 theta,
1433 1435 phi,
... ... @@ -1436,16 +1438,24 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
1436 1438 logEtrue);
1437 1439 GIntegral integral(&integrand);
1438 1440 integral.eps(1.0e-4);
  1441 + #endif
1439 1442  
1440 1443 // Loop over reconstructed energy
1441 1444 for (int ireco = 0; ireco < nreco; ++ireco) {
1442 1445  
1443   - // Get log of reconstructed energy boundaries
  1446 + #if defined(G_RMF_INTEGRATION)
  1447 + // Get log of reconstructed energy boundaries in MeV
1444 1448 double ereco_min = std::log(ereco.emin(ireco).MeV());
1445 1449 double ereco_max = std::log(ereco.emax(ireco).MeV());
1446 1450  
1447 1451 // Do Romberg integration
1448 1452 double value = integral.romberg(ereco_min, ereco_max);
  1453 + #else
  1454 + double value = response->edisp()->prob_erecobin(ereco.emin(ireco),
  1455 + ereco.emax(ireco),
  1456 + etrue.elogmean(itrue),
  1457 + theta);
  1458 + #endif
1449 1459  
1450 1460 // Update Rmf value and weight
1451 1461 m_rmf(itrue, ireco) += value * aeff;
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @0e4dbb5
... ... @@ -1818,7 +1818,7 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
1818 1818  
1819 1819  
1820 1820 /***********************************************************************//**
1821   - * @brief Return energy dispersion (in units or MeV^-1)
  1821 + * @brief Return energy dispersion (in units of MeV^-1)
1822 1822 *
1823 1823 * @param[in] obsEng Measured event energy.
1824 1824 * @param[in] theta Radial offset angle in camera (radians).
... ...
inst/cta/src/GCTAResponseTable.cpp
  View file @0e4dbb5
... ... @@ -609,7 +609,7 @@ double GCTAResponseTable::operator()(const int&amp; table,
609 609 const double& arg2,
610 610 const double& arg3) const
611 611 {
612   - // Throw exception if table is not 2D
  612 + // Throw exception if table is not 3D
613 613 if (axes() != 3) {
614 614 std::string msg = "Invalid response table dimension "+
615 615 gammalib::str(axes())+" encountered. Response "
... ... @@ -2020,30 +2020,43 @@ void GCTAResponseTable::update(const double&amp; arg1, const double&amp; arg2,
2020 2020 // Compute offsets
2021 2021 int size1 = axis_bins(0);
2022 2022 int size2 = axis_bins(1);
  2023 + int size12 = size1 * size2;
2023 2024 int offset_left_2 = nodes2->inx_left() * size1;
2024 2025 int offset_right_2 = nodes2->inx_right() * size1;
2025   - int offset_left_3 = nodes3->inx_left() * size1 * size2;
2026   - int offset_right_3 = nodes3->inx_right() * size1 * size2;
  2026 + int offset_left_3 = nodes3->inx_left() * size12;
  2027 + int offset_right_3 = nodes3->inx_right() * size12;
  2028 +
  2029 + // Pre-compute stuff
  2030 + int inx1l2l = nodes1->inx_left() + offset_left_2;
  2031 + int inx1l2r = nodes1->inx_left() + offset_right_2;
  2032 + int inx1r2l = nodes1->inx_right() + offset_left_2;
  2033 + int inx1r2r = nodes1->inx_right() + offset_right_2;
2027 2034  
2028 2035 // Set indices for tri-linear interpolation
2029   - m_inx1 = nodes1->inx_left() + offset_left_2 + offset_left_3 ;
2030   - m_inx2 = nodes1->inx_left() + offset_left_2 + offset_right_3;
2031   - m_inx3 = nodes1->inx_left() + offset_right_2 + offset_left_3 ;
2032   - m_inx4 = nodes1->inx_left() + offset_right_2 + offset_right_3;
2033   - m_inx5 = nodes1->inx_right() + offset_left_2 + offset_left_3 ;
2034   - m_inx6 = nodes1->inx_right() + offset_left_2 + offset_right_3;
2035   - m_inx7 = nodes1->inx_right() + offset_right_2 + offset_left_3 ;
2036   - m_inx8 = nodes1->inx_right() + offset_right_2 + offset_right_3;
  2036 + m_inx1 = inx1l2l + offset_left_3;
  2037 + m_inx2 = inx1l2l + offset_right_3;
  2038 + m_inx3 = inx1l2r + offset_left_3;
  2039 + m_inx4 = inx1l2r + offset_right_3;
  2040 + m_inx5 = inx1r2l + offset_left_3;
  2041 + m_inx6 = inx1r2l + offset_right_3;
  2042 + m_inx7 = inx1r2r + offset_left_3;
  2043 + m_inx8 = inx1r2r + offset_right_3;
  2044 +
  2045 + // Pre-compute stuff
  2046 + double wgt1l2l = nodes1->wgt_left() * nodes2->wgt_left();
  2047 + double wgt1l2r = nodes1->wgt_left() * nodes2->wgt_right();
  2048 + double wgt1r2l = nodes1->wgt_right() * nodes2->wgt_left();
  2049 + double wgt1r2r = nodes1->wgt_right() * nodes2->wgt_right();
2037 2050  
2038 2051 // Set weighting factors for tri-linear interpolation
2039   - m_wgt1 = nodes1->wgt_left() * nodes2->wgt_left() * nodes3->wgt_left();
2040   - m_wgt2 = nodes1->wgt_left() * nodes2->wgt_left() * nodes3->wgt_right();
2041   - m_wgt3 = nodes1->wgt_left() * nodes2->wgt_right() * nodes3->wgt_left();
2042   - m_wgt4 = nodes1->wgt_left() * nodes2->wgt_right() * nodes3->wgt_right();
2043   - m_wgt5 = nodes1->wgt_right() * nodes2->wgt_left() * nodes3->wgt_left();
2044   - m_wgt6 = nodes1->wgt_right() * nodes2->wgt_left() * nodes3->wgt_right();
2045   - m_wgt7 = nodes1->wgt_right() * nodes2->wgt_right() * nodes3->wgt_left();
2046   - m_wgt8 = nodes1->wgt_right() * nodes2->wgt_right() * nodes3->wgt_right();
  2052 + m_wgt1 = wgt1l2l * nodes3->wgt_left();
  2053 + m_wgt2 = wgt1l2l * nodes3->wgt_right();
  2054 + m_wgt3 = wgt1l2r * nodes3->wgt_left();
  2055 + m_wgt4 = wgt1l2r * nodes3->wgt_right();
  2056 + m_wgt5 = wgt1r2l * nodes3->wgt_left();
  2057 + m_wgt6 = wgt1r2l * nodes3->wgt_right();
  2058 + m_wgt7 = wgt1r2r * nodes3->wgt_left();
  2059 + m_wgt8 = wgt1r2r * nodes3->wgt_right();
2047 2060  
2048 2061 // Return
2049 2062 return;
... ...