gammalib / gammalib

Commit 2791cc64b148a2b43cd71d8795de405ba7f5ee27

Authored by Jürgen Knödlseder
1 parent e7137bed

Implement iterative adjustment of energy-dependent activation rate (#4209)

Inline Side-by-side
Showing 2 changed files with 407 additions and 111 deletions
inst/com/include/GCOMDris.hpp
  View file @2791cc6
... ... @@ -83,7 +83,7 @@ public:
83 83 const GCOMSelection& select = GCOMSelection(),
84 84 const double& zetamin = 5.0,
85 85 const double& timebin = 300.0,
86   - const std::string& method = "vetorate");
  86 + const std::string& method = "phibar");
87 87 std::string print(const GChatter& chatter = NORMAL) const;
88 88  
89 89 // Likelihood function
... ... @@ -101,13 +101,13 @@ public:
101 101 protected:
102 102 int m_ieng; //!< DRW energy bin
103 103 double m_norm; //!< Normalisation value
104   - int m_noads; //!< Number of superpackets
  104 + int m_nsp; //!< Number of superpackets
105 105 int m_nphibar; //!< Number of phibar layers
106 106 GNdarray m_vetorate; //!< Vetorate array multiplied by EHA cut correction
107   - GNdarray m_constrate; //!< Constant rate array multiplied by EHA cut correction
108   - GNdarray m_diffrate; //!< Vetorate - Constant rate array
  107 + GNdarray m_activrate; //!< Activation rate array multiplied by EHA cut correction
  108 + GNdarray m_diffrate; //!< Vetorate - activation rate array
109 109 GNdarray m_vetorate_sum; //!< Time integrated vetorate array
110   - GNdarray m_constrate_sum; //!< Time integrated constant rate array
  110 + GNdarray m_activrate_sum; //!< Time integrated activation rate array
111 111 GNdarray m_diffrate_sum; //!< Time integrated difference rate array
112 112 double m_value; //!< Function value
113 113 GVector m_gradient; //!< Gradient vector
... ... @@ -139,9 +139,11 @@ protected:
139 139 const GCOMSelection& select,
140 140 const double& zetamin);
141 141 void vetorate_fit(void);
  142 + void vetorate_update_activ(void);
142 143 void vetorate_generate(const GCOMObservation& obs,
143 144 const GCOMSelection& select,
144 145 const double& zetamin);
  146 + void vetorate_finish(const GCOMObservation& obs);
145 147 void vetorate_save(const GFilename& filename) const;
146 148 void vetorate_load(const GFilename& filename);
147 149  
... ... @@ -152,7 +154,7 @@ protected:
152 154 GNdarray m_wrk_counts; //!< 3D event cube array
153 155 GNdarray m_wrk_ehacutcorr; //!< 2D geometry response array
154 156 GNdarray m_wrk_vetorate; //!< 1D vetorates array
155   - GNdarray m_wrk_constrate; //!< 1D constant rate array
  157 + GNdarray m_wrk_activrate; //!< 2D activation rate array
156 158 GNdarray m_wrk_rate; //!< 2D rate array
157 159 std::vector<bool> m_wrk_use_sp; //!< 1D superpacket usage array
158 160 };
... ...
inst/com/src/GCOMDris.cpp
  View file @2791cc6
... ... @@ -35,9 +35,12 @@
35 35 #include "GFits.hpp"
36 36 #include "GFitsImageDouble.hpp"
37 37 #include "GFitsImageShort.hpp"
  38 +#include "GFitsTableLongCol.hpp"
  39 +#include "GFitsTableDoubleCol.hpp"
38 40 #include "GOptimizerPars.hpp"
39 41 #include "GOptimizerPar.hpp"
40 42 #include "GOptimizerLM.hpp"
  43 +#include "GFft.hpp"
41 44 #include "GCOMTools.hpp"
42 45 #include "GCOMSupport.hpp"
43 46 #include "GCOMDri.hpp"
... ... @@ -64,8 +67,8 @@
64 67 #define G_COMPUTE_DRWS_EHACORR //!< Correct Earth horizon cut
65 68  
66 69 /* __ Debug definitions __________________________________________________ */
67   -#define G_DEBUG_COMPUTE_DRWS
68   -#define G_DEBUG_SAVE_WORKING_ARRAYS
  70 +//#define G_DEBUG_COMPUTE_DRWS
  71 +//#define G_DEBUG_SAVE_WORKING_ARRAYS
69 72  
70 73 /* __ Constants __________________________________________________________ */
71 74  
... ... @@ -531,7 +534,7 @@ void GCOMDris::init_members(void)
531 534 m_wrk_counts.clear();
532 535 m_wrk_ehacutcorr.clear();
533 536 m_wrk_vetorate.clear();
534   - m_wrk_constrate.clear();
  537 + m_wrk_activrate.clear();
535 538 m_wrk_rate.clear();
536 539 m_wrk_use_sp.clear();
537 540  
... ... @@ -554,7 +557,7 @@ void GCOMDris::copy_members(const GCOMDris&amp; dris)
554 557 m_wrk_counts = dris.m_wrk_counts;
555 558 m_wrk_ehacutcorr = dris.m_wrk_ehacutcorr;
556 559 m_wrk_vetorate = dris.m_wrk_vetorate;
557   - m_wrk_constrate = dris.m_wrk_constrate;
  560 + m_wrk_activrate = dris.m_wrk_activrate;
558 561 m_wrk_rate = dris.m_wrk_rate;
559 562 m_wrk_use_sp = dris.m_wrk_use_sp;
560 563  
... ... @@ -1264,12 +1267,24 @@ void GCOMDris::compute_drws_vetorate(const GCOMObservation&amp; obs,
1264 1267 vetorate_load("debug_gcomdris_compute_drws_vetorate.fits");
1265 1268 #endif
1266 1269  
  1270 + // Initialise fprompt for all energies
  1271 + for (int i = 0; i < size(); ++i) {
  1272 + m_dris[i].m_drw_fprompt = 0.5; //!< Initial fprompt value
  1273 + }
  1274 +
1267 1275 // Fit working arrays to determine f_prompt parameter
1268 1276 vetorate_fit();
  1277 + for (int iter = 0; iter < 4; ++iter) { //!< Perform 4 iterations
  1278 + vetorate_update_activ();
  1279 + vetorate_fit();
  1280 + }
1269 1281  
1270 1282 // Generate DRWs
1271 1283 vetorate_generate(obs, select, zetamin);
1272 1284  
  1285 + // Finish DRWs
  1286 + vetorate_finish(obs);
  1287 +
1273 1288 // Return
1274 1289 return;
1275 1290 }
... ... @@ -1283,7 +1298,7 @@ void GCOMDris::compute_drws_vetorate(const GCOMObservation&amp; obs,
1283 1298 * @param[in] select Selection set.
1284 1299 * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
1285 1300 *
1286   - * Setup working arrays for vetorate DRW computation. There are four working
  1301 + * Setup working arrays for vetorate DRW computation. There are five working
1287 1302 * arrays:
1288 1303 *
1289 1304 * @c m_wrk_counts is a 3D working array spanned by superpacket index,
... ... @@ -1298,6 +1313,9 @@ void GCOMDris::compute_drws_vetorate(const GCOMObservation&amp; obs,
1298 1313 * @c m_wrk_vetorate is a 1D working array, containing the veto rate as
1299 1314 * function of superpacket index.
1300 1315 *
  1316 + * @c m_wrk_activrate is a 2D working array, containing the activation rate
  1317 + * as function of superpacket index and energy range.
  1318 + *
1301 1319 * @c m_wrk_use_sp is a 1D working array, signalling the usage of a given
1302 1320 * superpacket.
1303 1321 *
... ... @@ -1329,7 +1347,7 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1329 1347  
1330 1348 // Debug
1331 1349 #if defined(G_DEBUG_COMPUTE_DRWS)
1332   - std::cout << "Number of superpackets: " << noads << std::endl;
  1350 + std::cout << "Number of OAD superpackets: " << noads << std::endl;
1333 1351 std::cout << "Number of pixels: " << npix << std::endl;
1334 1352 std::cout << "Number of Phibar layers: " << nphibar << std::endl;
1335 1353 std::cout << "Number of energies: " << neng << std::endl;
... ... @@ -1339,8 +1357,14 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1339 1357 m_wrk_counts = GNdarray(noads, nphibar, neng);
1340 1358 m_wrk_ehacutcorr = GNdarray(noads, nphibar);
1341 1359 m_wrk_vetorate = GNdarray(noads);
  1360 + m_wrk_activrate = GNdarray(noads, neng);
1342 1361 m_wrk_use_sp = std::vector<bool>(noads);
1343 1362  
  1363 + // Initialise superpacket usage array
  1364 + for (int i_oad = 0; i_oad < noads; ++i_oad) {
  1365 + m_wrk_use_sp[i_oad] = false;
  1366 + }
  1367 +
1344 1368 // Allocate further working arrays
1345 1369 GNdarray geo(npix);
1346 1370 GNdarray eha(npix);
... ... @@ -1384,8 +1408,9 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1384 1408 tim.reduce(select.pulsar().validity());
1385 1409 }
1386 1410  
1387   - // Initialise event counter
  1411 + // Initialise event and superpacket counters
1388 1412 int i_evt = 0;
  1413 + int i_sp = 0;
1389 1414  
1390 1415 // Signal that loop should be terminated
1391 1416 bool terminate = false;
... ... @@ -1398,10 +1423,9 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1398 1423  
1399 1424 // Check superpacket usage for all DRWs so that their statistics
1400 1425 // get correctly updated
1401   - m_wrk_use_sp[i_oad] = true;
1402 1426 for (int i = 0; i < size(); ++i) {
1403   - if (!m_dris[i].use_superpacket(oad, tim, select)) {
1404   - m_wrk_use_sp[i_oad] = false;
  1427 + if (m_dris[i].use_superpacket(oad, tim, select)) {
  1428 + m_wrk_use_sp[i_oad] = true;
1405 1429 }
1406 1430 }
1407 1431 if (!m_wrk_use_sp[i_oad]) {
... ... @@ -1459,15 +1483,20 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1459 1483 accepted_geo += geo(index);
1460 1484 }
1461 1485 }
1462   - m_wrk_ehacutcorr(i_oad, iphibar) = accepted_geo / total_geo;
  1486 + m_wrk_ehacutcorr(i_sp, iphibar) = accepted_geo / total_geo;
1463 1487 }
1464 1488 }
1465 1489  
1466   - // Compute veto rate for superpacket time if rates are available
  1490 + // Compute veto rate and activation rates for superpacket time if rates
  1491 + // are available. The activation rates are initialised with a constant
  1492 + // rate.
1467 1493 if (!veto_times.is_empty()) {
1468 1494 double value = veto_times.interpolate(time, veto_rates);
1469 1495 if (value > 0.0) {
1470   - m_wrk_vetorate(i_oad) = value;
  1496 + m_wrk_vetorate(i_sp) = value;
  1497 + for (int ieng = 0; ieng < neng; ++ieng) {
  1498 + m_wrk_activrate(i_sp, ieng) = 1.0;
  1499 + }
1471 1500 }
1472 1501 }
1473 1502  
... ... @@ -1544,19 +1573,23 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1544 1573 }
1545 1574  
1546 1575 // Now fill the event in the 3D counts array
1547   - m_wrk_counts(i_oad, iphibar, ienergy) += 1.0;
  1576 + m_wrk_counts(i_sp, iphibar, ienergy) += 1.0;
1548 1577  
1549 1578 } // endfor: collected events
1550 1579  
1551 1580 // Debug
1552 1581 #if defined(G_DEBUG_COMPUTE_DRWS)
1553 1582 std::cout << "Superpacket " << i_oad;
  1583 + std::cout << " i_sp=" << i_sp;
1554 1584 std::cout << " time=" << time;
1555 1585 std::cout << " vetorate=" << m_wrk_vetorate(i_oad);
1556 1586 std::cout << " total_geo=" << total_geo;
1557 1587 std::cout << std::endl;
1558 1588 #endif
1559 1589  
  1590 + // Increment superpacket counter
  1591 + i_sp++;
  1592 +
1560 1593 // Break Orbit Aspect Data loop if termination was signalled or if there
1561 1594 // are no more events
1562 1595 if (terminate || i_evt >= events->size()) {
... ... @@ -1565,6 +1598,28 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1565 1598  
1566 1599 } // endfor: looped over Orbit Aspect Data
1567 1600  
  1601 + // Copy results in shortened Ndarrays
  1602 + GNdarray wrk_counts(i_sp, nphibar, neng);
  1603 + GNdarray wrk_ehacutcorr(i_sp, nphibar);
  1604 + GNdarray wrk_vetorate(i_sp);
  1605 + GNdarray wrk_activrate(i_sp, neng);
  1606 + for (int i = 0; i < i_sp; ++i) {
  1607 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1608 + for (int ieng = 0; ieng < neng; ++ieng) {
  1609 + wrk_counts(i, iphibar, ieng) = m_wrk_counts(i, iphibar, ieng);
  1610 + }
  1611 + wrk_ehacutcorr(i, iphibar) = m_wrk_ehacutcorr(i, iphibar);
  1612 + }
  1613 + for (int ieng = 0; ieng < neng; ++ieng) {
  1614 + wrk_activrate(i, ieng) = m_wrk_activrate(i, ieng);
  1615 + }
  1616 + wrk_vetorate(i) = m_wrk_vetorate(i);
  1617 + }
  1618 + m_wrk_counts = wrk_counts;
  1619 + m_wrk_ehacutcorr = wrk_ehacutcorr;
  1620 + m_wrk_vetorate = wrk_vetorate;
  1621 + m_wrk_activrate = wrk_activrate;
  1622 +
1568 1623 // Return
1569 1624 return;
1570 1625 }
... ... @@ -1581,8 +1636,7 @@ void GCOMDris::vetorate_setup(const GCOMObservation&amp; obs,
1581 1636 void GCOMDris::vetorate_fit(void)
1582 1637 {
1583 1638 // Set constants
1584   - const double fprompt_init = 0.5; // Initial f_prompt value
1585   - const double norm = 1.0; // Model normalisation
  1639 + const double norm = 1.0; // Model normalisation
1586 1640  
1587 1641 // Debug
1588 1642 #if defined(G_DEBUG_COMPUTE_DRWS)
... ... @@ -1591,39 +1645,45 @@ void GCOMDris::vetorate_fit(void)
1591 1645 #endif
1592 1646  
1593 1647 // Extract dimension from working arrays
1594   - int noads = m_wrk_counts.shape()[0];
  1648 + int nsp = m_wrk_counts.shape()[0];
1595 1649 int nphibar = m_wrk_counts.shape()[1];
1596 1650 int neng = m_wrk_counts.shape()[2];
1597 1651  
1598 1652 // Debug
1599 1653 #if defined(G_DEBUG_COMPUTE_DRWS)
1600   - std::cout << "Number of superpackets: " << noads << std::endl;
  1654 + std::cout << "Number of used superpackets: " << nsp << std::endl;
1601 1655 std::cout << "Number of Phibar layers: " << nphibar << std::endl;
1602 1656 std::cout << "Number of energies: " << neng << std::endl;
1603 1657 #endif
1604 1658  
1605 1659 // Setup rate result array
1606   - m_wrk_rate = GNdarray(noads, neng);
  1660 + m_wrk_rate = GNdarray(nsp, neng);
1607 1661  
1608   - // Setup normalised constant rate array and normalise veto rates
1609   - m_wrk_constrate = GNdarray(noads);
1610   - double nvetorate = 0.0;
1611   - double nconstrate = 0.0;
1612   - for (int ioad = 0; ioad < noads; ++ioad) {
1613   - if (m_wrk_vetorate(ioad) > 0.0) {
1614   - m_wrk_constrate(ioad) = 1.0;
1615   - nvetorate += m_wrk_vetorate(ioad);
1616   - nconstrate += m_wrk_constrate(ioad);
  1662 + // Normalise veto rates
  1663 + double nvetorate = 0.0;
  1664 + for (int isp = 0; isp < nsp; ++isp) {
  1665 + if (m_wrk_vetorate(isp) > 0.0) {
  1666 + nvetorate += m_wrk_vetorate(isp);
1617 1667 }
1618 1668 }
1619 1669 if (nvetorate > 0.0) {
1620   - for (int ioad = 0; ioad < noads; ++ioad) {
1621   - m_wrk_vetorate(ioad) /= nvetorate;
  1670 + for (int isp = 0; isp < nsp; ++isp) {
  1671 + m_wrk_vetorate(isp) /= nvetorate;
1622 1672 }
1623 1673 }
1624   - if (nconstrate > 0.0) {
1625   - for (int ioad = 0; ioad < noads; ++ioad) {
1626   - m_wrk_constrate(ioad) /= nconstrate;
  1674 +
  1675 + // Normalise activation rates
  1676 + for (int ieng = 0; ieng < neng; ++ieng) {
  1677 + double nactivrate = 0.0;
  1678 + for (int isp = 0; isp < nsp; ++isp) {
  1679 + if (m_wrk_activrate(isp, ieng) > 0.0) {
  1680 + nactivrate += m_wrk_activrate(isp, ieng);
  1681 + }
  1682 + }
  1683 + if (nactivrate > 0.0) {
  1684 + for (int isp = 0; isp < nsp; ++isp) {
  1685 + m_wrk_activrate(isp, ieng) /= nactivrate;
  1686 + }
1627 1687 }
1628 1688 }
1629 1689  
... ... @@ -1632,7 +1692,7 @@ void GCOMDris::vetorate_fit(void)
1632 1692  
1633 1693 // Set fprompt parameter
1634 1694 GOptimizerPars pars;
1635   - GOptimizerPar par("fprompt", fprompt_init);
  1695 + GOptimizerPar par("fprompt", m_dris[ieng].m_drw_fprompt);
1636 1696 par.range(0.0, 1.0);
1637 1697 par.factor_gradient(1.0); // To avoid zero gradient log message
1638 1698 pars.append(par);
... ... @@ -1659,15 +1719,15 @@ void GCOMDris::vetorate_fit(void)
1659 1719 // Recover fprompt parameter and errors
1660 1720 double fprompt = pars[0]->value();
1661 1721 double e_fprompt = pars[0]->error();
1662   - double fconst = 1.0 - fprompt;
  1722 + double factiv = 1.0 - fprompt;
1663 1723  
1664 1724 // Set resulting rate vector for this energy bin
1665   - for (int ioad = 0; ioad < noads; ++ioad) {
1666   - m_wrk_rate(ioad, ieng) = fprompt * m_wrk_vetorate(ioad) +
1667   - fconst * m_wrk_constrate(ioad);
  1725 + for (int isp = 0; isp < nsp; ++isp) {
  1726 + m_wrk_rate(isp, ieng) = fprompt * m_wrk_vetorate(isp) +
  1727 + factiv * m_wrk_activrate(isp, ieng);
1668 1728 }
1669 1729  
1670   - // Set DRW header parameters
  1730 + // Set DRW header members
1671 1731 m_dris[ieng].m_drw_status = opt.status_string();
1672 1732 m_dris[ieng].m_drw_fprompt = fprompt;
1673 1733 m_dris[ieng].m_drw_e_fprompt = e_fprompt;
... ... @@ -1689,6 +1749,136 @@ void GCOMDris::vetorate_fit(void)
1689 1749  
1690 1750  
1691 1751 /***********************************************************************//**
  1752 + * @brief Update activation rate
  1753 + ***************************************************************************/
  1754 +void GCOMDris::vetorate_update_activ(void)
  1755 +{
  1756 + // Set constants
  1757 + const double norm = 1.0; // Model normalisation
  1758 + const double t_hour = 1.0; // Smoothing duration (hours)
  1759 +
  1760 + // Debug
  1761 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1762 + std::cout << "GCOMDris::vetorate_update_activ" << std::endl;
  1763 + std::cout << "===============================" << std::endl;
  1764 + #endif
  1765 +
  1766 + // Extract dimension from working arrays
  1767 + int nsp = m_wrk_counts.shape()[0];
  1768 + int nphibar = m_wrk_counts.shape()[1];
  1769 + int neng = m_wrk_counts.shape()[2];
  1770 +
  1771 + // Setup Gaussian smoothing kernel
  1772 + double nsmooth = 3600.0/16.384 * t_hour;
  1773 + double weight = -0.5 / (nsmooth * nsmooth);
  1774 + GNdarray kernel(nsp);
  1775 + double sum = 0.0;
  1776 + for (int isp = 0; isp < nsp/2; ++isp) {
  1777 + double kvalue = std::exp(weight*isp*isp);
  1778 + if (kvalue <= 0.0) {
  1779 + break;
  1780 + }
  1781 + kernel(isp) = kvalue;
  1782 + sum += kvalue;
  1783 + if (isp > 0) {
  1784 + kernel(nsp-isp) = kvalue;
  1785 + sum += kvalue;
  1786 + }
  1787 + }
  1788 + if (sum > 0.0) {
  1789 + for (int isp = 0; isp < nsp; ++isp) {
  1790 + kernel(isp) /= sum;
  1791 + }
  1792 + }
  1793 + GFft fft_kernel(kernel);
  1794 +
  1795 + // Loop over energy bins
  1796 + for (int ieng = 0; ieng < neng; ++ieng) {
  1797 +
  1798 + // Compute model normalised to number of events for each Phibar
  1799 + GNdarray model(nsp, nphibar);
  1800 + GNdarray n(nphibar);
  1801 + for (int isp = 0; isp < nsp; ++isp) {
  1802 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1803 + model(isp, iphibar) = m_wrk_rate(isp, ieng) * m_wrk_ehacutcorr(isp, iphibar);
  1804 + }
  1805 + }
  1806 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1807 + double nevents = 0.0;
  1808 + double nmodel = 0.0;
  1809 + for (int isp = 0; isp < nsp; ++isp) {
  1810 + if (model(isp, iphibar) > 0.0) {
  1811 + nevents += m_wrk_counts(isp, iphibar, ieng);
  1812 + nmodel += model(isp, iphibar);
  1813 + }
  1814 + }
  1815 + if (nmodel > 0.0) {
  1816 + n(iphibar) = norm * (nevents / nmodel);
  1817 + for (int isp = 0; isp < nsp; ++isp) {
  1818 + model(isp, iphibar) *= n(iphibar);
  1819 + }
  1820 + }
  1821 + }
  1822 +
  1823 + // Compute residual and Earth horizon correction cut by summing over
  1824 + // all Phibar layers
  1825 + GNdarray residual(nsp);
  1826 + GNdarray ehacutcorr(nsp);
  1827 + for (int isp = 0; isp < nsp; ++isp) {
  1828 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1829 + residual(isp) += m_wrk_counts(isp, iphibar, ieng) - model(isp, iphibar);
  1830 + ehacutcorr(isp) += n(iphibar) * m_wrk_ehacutcorr(isp, iphibar);
  1831 + }
  1832 + }
  1833 +
  1834 + // Smooth residual
  1835 + GFft fft_residual(residual);
  1836 + fft_residual = fft_residual * fft_kernel;
  1837 + residual = fft_residual.backward();
  1838 +
  1839 + // Smooth ehacutcorr
  1840 + GFft fft_ehacutcorr(ehacutcorr);
  1841 + fft_ehacutcorr = fft_ehacutcorr * fft_kernel;
  1842 + ehacutcorr = fft_ehacutcorr.backward();
  1843 +
  1844 + // Derive indicent residual background rate
  1845 + for (int isp = 0; isp < nsp; ++isp) {
  1846 + if (ehacutcorr(isp) > 0.0) {
  1847 + residual(isp) /= ehacutcorr(isp);
  1848 + }
  1849 + }
  1850 +
  1851 + // Update activation rate
  1852 + for (int isp = 0; isp < nsp; ++isp) {
  1853 + double activ = m_wrk_activrate(isp, ieng);
  1854 + if (activ > 0.0) {
  1855 + double update = activ + residual(isp);
  1856 + if (update < 0.05 * activ) { //!< Cap reduction at 5%
  1857 + update = 0.05 * activ;
  1858 + }
  1859 + m_wrk_activrate(isp, ieng) = update;
  1860 + }
  1861 + }
  1862 +
  1863 + // Normalise activation rate
  1864 + double nactivrate = 0.0;
  1865 + for (int isp = 0; isp < nsp; ++isp) {
  1866 + nactivrate += m_wrk_activrate(isp, ieng);
  1867 + }
  1868 + if (nactivrate > 0.0) {
  1869 + for (int isp = 0; isp < nsp; ++isp) {
  1870 + m_wrk_activrate(isp, ieng) /= nactivrate;
  1871 + }
  1872 + }
  1873 +
  1874 + } // endfor: looped over all energy bins
  1875 +
  1876 + // Return
  1877 + return;
  1878 +}
  1879 +
  1880 +
  1881 +/***********************************************************************//**
1692 1882 * @brief Generate DRWs
1693 1883 *
1694 1884 * @param[in] obs COMPTEL observation.
... ... @@ -1719,7 +1909,7 @@ void GCOMDris::vetorate_generate(const GCOMObservation&amp; obs,
1719 1909  
1720 1910 // Debug
1721 1911 #if defined(G_DEBUG_COMPUTE_DRWS)
1722   - std::cout << "Number of superpackets: " << noads << std::endl;
  1912 + std::cout << "Number of OAD superpackets: " << noads << std::endl;
1723 1913 std::cout << "Number of pixels: " << npix << std::endl;
1724 1914 std::cout << "Number of Phibar layers: " << nphibar << std::endl;
1725 1915 std::cout << "Number of energies: " << neng << std::endl;
... ... @@ -1747,6 +1937,9 @@ void GCOMDris::vetorate_generate(const GCOMObservation&amp; obs,
1747 1937 tim.reduce(select.pulsar().validity());
1748 1938 }
1749 1939  
  1940 + // Initialise superpacket counter
  1941 + int i_sp = 0;
  1942 +
1750 1943 // Loop over Orbit Aspect Data
1751 1944 for (int i_oad = 0; i_oad < noads; ++i_oad) {
1752 1945  
... ... @@ -1814,7 +2007,7 @@ void GCOMDris::vetorate_generate(const GCOMObservation&amp; obs,
1814 2007 int inx = index + iphibar * npix;
1815 2008  
1816 2009 // Store geometry weighted with rate as DRW
1817   - m_dris[i][inx] += geometry(index) * m_wrk_rate(i_oad, i);
  2010 + m_dris[i][inx] += geometry(index) * m_wrk_rate(i_sp, i);
1818 2011  
1819 2012 } // endif: Earth horizon cut passed
1820 2013  
... ... @@ -1822,12 +2015,12 @@ void GCOMDris::vetorate_generate(const GCOMObservation&amp; obs,
1822 2015  
1823 2016 } // endfor: looped over Phibar
1824 2017  
  2018 + // Set binary table information
  2019 +
1825 2020 } // endfor: looped over all DRWs
1826 2021  
1827   - // Debug
1828   - #if defined(G_DEBUG_COMPUTE_DRWS)
1829   - std::cout << "Superpacket " << i_oad << std::endl;
1830   - #endif
  2022 + // Increment superpacket counter
  2023 + i_sp++;
1831 2024  
1832 2025 } // endfor: looped over Orbit Aspect Data
1833 2026  
... ... @@ -1837,6 +2030,93 @@ void GCOMDris::vetorate_generate(const GCOMObservation&amp; obs,
1837 2030  
1838 2031  
1839 2032 /***********************************************************************//**
  2033 + * @brief Finish DRWs
  2034 + *
  2035 + * @param[in] obs COMPTEL observation.
  2036 + *
  2037 + * Set DRW binary tables
  2038 + ***************************************************************************/
  2039 +void GCOMDris::vetorate_finish(const GCOMObservation& obs)
  2040 +{
  2041 + // Debug
  2042 + #if defined(G_DEBUG_COMPUTE_DRWS)
  2043 + std::cout << "GCOMDris::vetorate_finish" << std::endl;
  2044 + std::cout << "=========================" << std::endl;
  2045 + #endif
  2046 +
  2047 + // Extract dimensions
  2048 + int noads = obs.oads().size();
  2049 +
  2050 + // Compute number of to-be-used superpackets
  2051 + int nsp = 0;
  2052 + for (int i_oad = 0; i_oad < noads; ++i_oad) {
  2053 + if (m_wrk_use_sp[i_oad]) {
  2054 + nsp++;
  2055 + }
  2056 + }
  2057 +
  2058 + // Debug
  2059 + #if defined(G_DEBUG_COMPUTE_DRWS)
  2060 + std::cout << "Number of OAD superpackets: " << noads << std::endl;
  2061 + std::cout << "Number of used superpackets: " << nsp << std::endl;
  2062 + #endif
  2063 +
  2064 + // Setup DRW binary tables
  2065 + for (int i = 0; i < size(); ++i) {
  2066 +
  2067 + // Get pointer to DRW
  2068 + GCOMDri* dri = &(m_dris[i]);
  2069 +
  2070 + // Allocate columns
  2071 + GFitsTableLongCol tjd("TJD", nsp);
  2072 + GFitsTableLongCol tics("TICS", nsp);
  2073 + GFitsTableDoubleCol rate("RATE", nsp);
  2074 + GFitsTableDoubleCol vetorate("VETORATE", nsp);
  2075 + GFitsTableDoubleCol activrate("ACTIVRATE", nsp);
  2076 +
  2077 + // Initialise superpacket counter
  2078 + int isp = 0;
  2079 +
  2080 + // Loop over Orbit Aspect Data
  2081 + for (int i_oad = 0; i_oad < noads; ++i_oad) {
  2082 +
  2083 + // Get reference to Orbit Aspect Data of superpacket
  2084 + const GCOMOad &oad = obs.oads()[i_oad];
  2085 +
  2086 + // Skip not-to-be-used superpackets
  2087 + if (!m_wrk_use_sp[i_oad]) {
  2088 + continue;
  2089 + }
  2090 +
  2091 + // Fill table columns
  2092 + tjd(isp) = oad.tjd();
  2093 + tics(isp) = oad.tics();
  2094 + rate(isp) = m_wrk_rate(isp, i);
  2095 + vetorate(isp) = m_wrk_vetorate(isp);
  2096 + activrate(isp) = m_wrk_activrate(isp, i);
  2097 +
  2098 + // Increment superpacket counter
  2099 + isp++;
  2100 +
  2101 + } // endfor: looped over Orbit Aspect Data
  2102 +
  2103 + // Append columns to binary table
  2104 + dri->m_drw_table.clear();
  2105 + dri->m_drw_table.extname("RATES");
  2106 + dri->m_drw_table.append(tjd);
  2107 + dri->m_drw_table.append(tics);
  2108 + dri->m_drw_table.append(rate);
  2109 + dri->m_drw_table.append(vetorate);
  2110 + dri->m_drw_table.append(activrate);
  2111 +
  2112 + } // endfor: looped over DRWs
  2113 +
  2114 + // Return
  2115 + return;
  2116 +}
  2117 +
  2118 +
  2119 +/***********************************************************************//**
1840 2120 * @brief Save working arrays for vetorate computation
1841 2121 *
1842 2122 * @param[in] filename FITS filename.
... ... @@ -1852,32 +2132,40 @@ void GCOMDris::vetorate_save(const GFilename&amp; filename) const
1852 2132 #endif
1853 2133  
1854 2134 // Extract dimension from working arrays
1855   - int noads = m_wrk_counts.shape()[0];
  2135 + int noads = m_wrk_use_sp.size();
  2136 + int nsp = m_wrk_counts.shape()[0];
1856 2137 int nphibar = m_wrk_counts.shape()[1];
1857 2138 int neng = m_wrk_counts.shape()[2];
1858 2139  
1859 2140 // Allocate FITS images
1860   - GFitsImageDouble image_counts(noads, nphibar, neng);
1861   - GFitsImageDouble image_ehacutcorr(noads, nphibar);
1862   - GFitsImageDouble image_vetorate(noads);
  2141 + GFitsImageDouble image_counts(nsp, nphibar, neng);
  2142 + GFitsImageDouble image_ehacutcorr(nsp, nphibar);
  2143 + GFitsImageDouble image_vetorate(nsp);
  2144 + GFitsImageDouble image_activrate(nsp, neng);
1863 2145 GFitsImageShort image_use_sp(noads);
1864 2146  
1865 2147 // Fill FITS images from working arrays
1866   - for (int ioad = 0; ioad < noads; ++ioad) {
  2148 + for (int isp = 0; isp < nsp; ++isp) {
1867 2149 for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1868   - for (int ienergy = 0; ienergy < neng; ++ienergy) {
1869   - image_counts(ioad, iphibar, ienergy) = m_wrk_counts(ioad, iphibar, ienergy);
  2150 + for (int ieng = 0; ieng < neng; ++ieng) {
  2151 + image_counts(isp, iphibar, ieng) = m_wrk_counts(isp, iphibar, ieng);
1870 2152 }
1871   - image_ehacutcorr(ioad, iphibar) = m_wrk_ehacutcorr(ioad, iphibar);
  2153 + image_ehacutcorr(isp, iphibar) = m_wrk_ehacutcorr(isp, iphibar);
1872 2154 }
1873   - image_vetorate(ioad) = m_wrk_vetorate(ioad);
1874   - image_use_sp(ioad) = (m_wrk_use_sp[ioad]) ? 1 : 0;
  2155 + for (int ieng = 0; ieng < neng; ++ieng) {
  2156 + image_activrate(isp, ieng) = m_wrk_activrate(isp, ieng);
  2157 + }
  2158 + image_vetorate(isp) = m_wrk_vetorate(isp);
  2159 + }
  2160 + for (int ioad = 0; ioad < noads; ++ioad) {
  2161 + image_use_sp(ioad) = (m_wrk_use_sp[ioad]) ? 1 : 0;
1875 2162 }
1876 2163  
1877 2164 // Set image attributes
1878 2165 image_counts.extname("COUNTS");
1879 2166 image_ehacutcorr.extname("EHACUTCORR");
1880 2167 image_vetorate.extname("VETORATE");
  2168 + image_activrate.extname("ACTIVRATE");
1881 2169 image_use_sp.extname("SPUSAGE");
1882 2170  
1883 2171 // Create FITS file and append images
... ... @@ -1885,6 +2173,7 @@ void GCOMDris::vetorate_save(const GFilename&amp; filename) const
1885 2173 fits.append(image_counts);
1886 2174 fits.append(image_ehacutcorr);
1887 2175 fits.append(image_vetorate);
  2176 + fits.append(image_activrate);
1888 2177 fits.append(image_use_sp);
1889 2178  
1890 2179 // Save FITS file
... ... @@ -1917,36 +2206,45 @@ void GCOMDris::vetorate_load(const GFilename&amp; filename)
1917 2206 const GFitsImage* image_counts = fits.image("COUNTS");
1918 2207 const GFitsImage* image_ehacutcorr = fits.image("EHACUTCORR");
1919 2208 const GFitsImage* image_vetorate = fits.image("VETORATE");
  2209 + const GFitsImage* image_activrate = fits.image("ACTIVRATE");
1920 2210 const GFitsImage* image_use_sp = fits.image("SPUSAGE");
1921 2211  
1922 2212 // Extract dimensions
1923   - int noads = image_counts->naxes(0);
  2213 + int noads = image_use_sp->naxes(0);
  2214 + int nsp = image_counts->naxes(0);
1924 2215 int nphibar = image_counts->naxes(1);
1925 2216 int neng = image_counts->naxes(2);
1926 2217  
1927 2218 // Debug
1928 2219 #if defined(G_DEBUG_COMPUTE_DRWS)
1929   - std::cout << "Number of superpackets: " << noads << std::endl;
  2220 + std::cout << "Number of OAD superpackets: " << noads << std::endl;
  2221 + std::cout << "Number of used superpackets: " << nsp << std::endl;
1930 2222 std::cout << "Number of Phibar layers: " << nphibar << std::endl;
1931 2223 std::cout << "Number of energies: " << neng << std::endl;
1932 2224 #endif
1933 2225  
1934 2226 // Allocate working array
1935   - m_wrk_counts = GNdarray(noads, nphibar, neng);
1936   - m_wrk_ehacutcorr = GNdarray(noads, nphibar);
1937   - m_wrk_vetorate = GNdarray(noads);
  2227 + m_wrk_counts = GNdarray(nsp, nphibar, neng);
  2228 + m_wrk_ehacutcorr = GNdarray(nsp, nphibar);
  2229 + m_wrk_vetorate = GNdarray(nsp);
  2230 + m_wrk_activrate = GNdarray(nsp, neng);
1938 2231 m_wrk_use_sp = std::vector<bool>(noads);
1939 2232  
1940 2233 // Extract data from images
1941   - for (int ioad = 0; ioad < noads; ++ioad) {
  2234 + for (int isp = 0; isp < nsp; ++isp) {
1942 2235 for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1943   - for (int ienergy = 0; ienergy < neng; ++ienergy) {
1944   - m_wrk_counts(ioad, iphibar, ienergy) = image_counts->pixel(ioad, iphibar, ienergy);
  2236 + for (int ieng = 0; ieng < neng; ++ieng) {
  2237 + m_wrk_counts(isp, iphibar, ieng) = image_counts->pixel(isp, iphibar, ieng);
1945 2238 }
1946   - m_wrk_ehacutcorr(ioad, iphibar) = image_ehacutcorr->pixel(ioad, iphibar);
  2239 + m_wrk_ehacutcorr(isp, iphibar) = image_ehacutcorr->pixel(isp, iphibar);
  2240 + }
  2241 + for (int ieng = 0; ieng < neng; ++ieng) {
  2242 + m_wrk_activrate(isp, ieng) = image_activrate->pixel(isp, ieng);
1947 2243 }
1948   - m_wrk_vetorate(ioad) = image_vetorate->pixel(ioad);
1949   - m_wrk_use_sp[ioad] = (image_use_sp->pixel(ioad) == 1);
  2244 + m_wrk_vetorate(isp) = image_vetorate->pixel(isp);
  2245 + }
  2246 + for (int ioad = 0; ioad < noads; ++ioad) {
  2247 + m_wrk_use_sp[ioad] = (image_use_sp->pixel(ioad) == 1);
1950 2248 }
1951 2249  
1952 2250 // Close FITS file
... ... @@ -1981,36 +2279,36 @@ GCOMDris::likelihood::likelihood(GCOMDris *dris,
1981 2279 m_curvature = GMatrixSparse(1,1);
1982 2280  
1983 2281 // Extract dimension from working array
1984   - m_noads = m_this->m_wrk_counts.shape()[0];
  2282 + m_nsp = m_this->m_wrk_counts.shape()[0];
1985 2283 m_nphibar = m_this->m_wrk_counts.shape()[1];
1986 2284  
1987 2285 // Multiply veto and constant rate by EHA cut correction. Rates are zero
1988 2286 // in case that the vetorate is zero.
1989   - m_vetorate = GNdarray(m_noads, m_nphibar);
1990   - m_constrate = GNdarray(m_noads, m_nphibar);
1991   - m_diffrate = GNdarray(m_noads, m_nphibar);
1992   - for (int ioad = 0; ioad < m_noads; ++ioad) {
1993   - double vetorate = m_this->m_wrk_vetorate(ioad);
  2287 + m_vetorate = GNdarray(m_nsp, m_nphibar);
  2288 + m_activrate = GNdarray(m_nsp, m_nphibar);
  2289 + m_diffrate = GNdarray(m_nsp, m_nphibar);
  2290 + for (int isp = 0; isp < m_nsp; ++isp) {
  2291 + double vetorate = m_this->m_wrk_vetorate(isp);
1994 2292 if (vetorate > 0.0) {
1995   - double constrate = m_this->m_wrk_constrate(ioad);
  2293 + double activrate = m_this->m_wrk_activrate(isp, ieng);
1996 2294 for (int iphibar = 0; iphibar < m_nphibar; ++iphibar) {
1997   - m_vetorate(ioad, iphibar) = vetorate * m_this->m_wrk_ehacutcorr(ioad, iphibar);
1998   - m_constrate(ioad, iphibar) = constrate * m_this->m_wrk_ehacutcorr(ioad, iphibar);
1999   - m_diffrate(ioad, iphibar) = m_vetorate(ioad, iphibar) - m_constrate(ioad, iphibar);
  2295 + m_vetorate(isp, iphibar) = vetorate * m_this->m_wrk_ehacutcorr(isp, iphibar);
  2296 + m_activrate(isp, iphibar) = activrate * m_this->m_wrk_ehacutcorr(isp, iphibar);
  2297 + m_diffrate(isp, iphibar) = m_vetorate(isp, iphibar) - m_activrate(isp, iphibar);
2000 2298 }
2001 2299 }
2002 2300 }
2003 2301  
2004 2302 // Compute rate sums
2005 2303 m_vetorate_sum = GNdarray(m_nphibar);
2006   - m_constrate_sum = GNdarray(m_nphibar);
  2304 + m_activrate_sum = GNdarray(m_nphibar);
2007 2305 m_diffrate_sum = GNdarray(m_nphibar);
2008 2306 for (int iphibar = 0; iphibar < m_nphibar; ++iphibar) {
2009   - for (int ioad = 0; ioad < m_noads; ++ioad) {
2010   - m_vetorate_sum(iphibar) += m_vetorate(ioad, iphibar);
2011   - m_constrate_sum(iphibar) += m_constrate(ioad, iphibar);
  2307 + for (int isp = 0; isp < m_nsp; ++isp) {
  2308 + m_vetorate_sum(iphibar) += m_vetorate(isp, iphibar);
  2309 + m_activrate_sum(iphibar) += m_activrate(isp, iphibar);
2012 2310 }
2013   - m_diffrate_sum(iphibar) = m_vetorate_sum(iphibar) - m_constrate_sum(iphibar);
  2311 + m_diffrate_sum(iphibar) = m_vetorate_sum(iphibar) - m_activrate_sum(iphibar);
2014 2312 }
2015 2313  
2016 2314 // Return
... ... @@ -2030,44 +2328,40 @@ void GCOMDris::likelihood::eval(const GOptimizerPars&amp; pars)
2030 2328 {
2031 2329 // Recover parameters
2032 2330 double fprompt = pars[0]->value();
2033   - double fconst = 1.0 - fprompt;
2034   -
2035   - // Extract dimension from working array
2036   - int noads = m_this->m_wrk_counts.shape()[0];
2037   - int nphibar = m_this->m_wrk_counts.shape()[1];
  2331 + double factiv = 1.0 - fprompt;
2038 2332  
2039 2333 // Allocate phibar normalisation arrays
2040   - GNdarray nevents(nphibar);
2041   - GNdarray nmodel(nphibar);
2042   - GNdarray norm(nphibar);
  2334 + GNdarray nevents(m_nphibar);
  2335 + GNdarray nmodel(m_nphibar);
  2336 + GNdarray norm(m_nphibar);
2043 2337  
2044 2338 // Compute model
2045   - GNdarray model = fprompt * m_vetorate + fconst * m_constrate;
  2339 + GNdarray model = fprompt * m_vetorate + factiv * m_activrate;
2046 2340  
2047 2341 // Compute Phibar normalised model
2048 2342 GNdarray model_norm = model;
2049 2343 for (int iphibar = 0; iphibar < m_nphibar; ++iphibar) {
2050   - for (int ioad = 0; ioad < m_noads; ++ioad) {
2051   - if (model(ioad, iphibar) > 0.0) {
2052   - nevents(iphibar) += m_this->m_wrk_counts(ioad, iphibar, m_ieng);
2053   - nmodel(iphibar) += model(ioad, iphibar);
  2344 + for (int isp = 0; isp < m_nsp; ++isp) {
  2345 + if (model(isp, iphibar) > 0.0) {
  2346 + nevents(iphibar) += m_this->m_wrk_counts(isp, iphibar, m_ieng);
  2347 + nmodel(iphibar) += model(isp, iphibar);
2054 2348 }
2055 2349 }
2056 2350 if (nmodel(iphibar) > 0.0) {
2057 2351 norm(iphibar) = m_norm * nevents(iphibar) / nmodel(iphibar);
2058   - for (int ioad = 0; ioad < m_noads; ++ioad) {
2059   - model_norm(ioad, iphibar) *= norm(iphibar);
  2352 + for (int isp = 0; isp < m_nsp; ++isp) {
  2353 + model_norm(isp, iphibar) *= norm(iphibar);
2060 2354 }
2061 2355 }
2062 2356 }
2063 2357  
2064 2358 // Compute LogL
2065 2359 m_value = 0.0;
2066   - for (int ioad = 0; ioad < m_noads; ++ioad) {
  2360 + for (int isp = 0; isp < m_nsp; ++isp) {
2067 2361 for (int iphibar = 0; iphibar < m_nphibar; ++iphibar) {
2068   - if (model_norm(ioad, iphibar) > 0.0) {
2069   - m_value -= m_this->m_wrk_counts(ioad, iphibar, m_ieng) *
2070   - std::log(model_norm(ioad, iphibar)) - model_norm(ioad, iphibar);
  2362 + if (model_norm(isp, iphibar) > 0.0) {
  2363 + m_value -= m_this->m_wrk_counts(isp, iphibar, m_ieng) *
  2364 + std::log(model_norm(isp, iphibar)) - model_norm(isp, iphibar);
2071 2365 }
2072 2366 }
2073 2367 }
... ... @@ -2080,18 +2374,18 @@ void GCOMDris::likelihood::eval(const GOptimizerPars&amp; pars)
2080 2374 double g_norm = -m_norm * nevents(iphibar) / nmodel2 * m_diffrate_sum(iphibar);
2081 2375  
2082 2376 // Loop over superpackets
2083   - for (int ioad = 0; ioad < m_noads; ++ioad) {
  2377 + for (int isp = 0; isp < m_nsp; ++isp) {
2084 2378  
2085 2379 // Continue only if model is positive
2086   - if (model_norm(ioad, iphibar) > 0.0) {
  2380 + if (model_norm(isp, iphibar) > 0.0) {
2087 2381  
2088 2382 // Pre computation
2089   - double fb = m_this->m_wrk_counts(ioad, iphibar, m_ieng) / model_norm(ioad, iphibar);
  2383 + double fb = m_this->m_wrk_counts(isp, iphibar, m_ieng) / model_norm(isp, iphibar);
2090 2384 double fc = (1.0 - fb);
2091   - double fa = fb / model_norm(ioad, iphibar);
  2385 + double fa = fb / model_norm(isp, iphibar);
2092 2386  
2093 2387 // Compute parameter gradient
2094   - double g = norm(iphibar) * m_diffrate(ioad, iphibar) + g_norm * model(ioad, iphibar);
  2388 + double g = norm(iphibar) * m_diffrate(isp, iphibar) + g_norm * model(isp, iphibar);
2095 2389  
2096 2390 // Update gradient
2097 2391 m_gradient[0] += fc * g;
... ...