gammalib / gammalib

Commit 1fc91b4a17a5008f63b178ab79e19719bf286370

Authored by Jürgen Knödlseder
1 parent a423b3ab

Add vetorate method for DRW computation (#4209)

Inline Side-by-side
Showing 4 changed files with 1045 additions and 48 deletions
inst/com/include/GCOMDris.hpp
  View file @1fc91b4
... ... @@ -30,10 +30,13 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32 32 #include <vector>
33   -#include "GCOMDri.hpp"
34 33 #include "GContainer.hpp"
  34 +#include "GNdarray.hpp"
  35 +#include "GOptimizerFunction.hpp"
  36 +#include "GCOMDri.hpp"
35 37  
36 38 /* __ Forward declarations _______________________________________________ */
  39 +class GOptimizer;
37 40 class GCOMObservation;
38 41 class GCOMEventList;
39 42 class GCOMSelection;
... ... @@ -80,9 +83,38 @@ public:
80 83 const GCOMSelection& select = GCOMSelection(),
81 84 const double& zetamin = 5.0,
82 85 const double& timebin = 300.0,
83   - const std::string& method = "phibar");
  86 + const std::string& method = "vetorate");
84 87 std::string print(const GChatter& chatter = NORMAL) const;
85 88  
  89 + // Likelihood function
  90 + class likelihood : public GOptimizerFunction {
  91 + public:
  92 + // Constructors and destructors
  93 + likelihood(GCOMDris *dris, const int& ieng, const double& norm);
  94 +
  95 + // Implemented pure virtual base class methods
  96 + virtual void eval(const GOptimizerPars& pars);
  97 + virtual double value(void) const;
  98 + virtual GVector* gradient(void);
  99 + virtual GMatrixSparse* curvature(void);
  100 +
  101 + protected:
  102 + int m_ieng; //!< DRW energy bin
  103 + double m_norm; //!< Normalisation value
  104 + int m_noads; //!< Number of superpackets
  105 + int m_nphibar; //!< Number of phibar layers
  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
  109 + GNdarray m_vetorate_sum; //!< Time integrated vetorate array
  110 + GNdarray m_constrate_sum; //!< Time integrated constant rate array
  111 + GNdarray m_diffrate_sum; //!< Time integrated difference rate array
  112 + double m_value; //!< Function value
  113 + GVector m_gradient; //!< Gradient vector
  114 + GMatrixSparse m_curvature; //!< Curvature matrix
  115 + GCOMDris* m_this; //!< Pointer to GCOMDris object
  116 + };
  117 +
86 118 protected:
87 119 // Protected methods
88 120 void init_members(void);
... ... @@ -98,9 +130,31 @@ protected:
98 130 const GCOMSelection& select,
99 131 const double& zetamin,
100 132 const double& timebin);
  133 + void compute_drws_vetorate(const GCOMObservation& obs,
  134 + const GCOMEventList* events,
  135 + const GCOMSelection& select,
  136 + const double& zetamin);
  137 + void vetorate_setup(const GCOMObservation& obs,
  138 + const GCOMEventList* events,
  139 + const GCOMSelection& select,
  140 + const double& zetamin);
  141 + void vetorate_fit(void);
  142 + void vetorate_generate(const GCOMObservation& obs,
  143 + const GCOMSelection& select,
  144 + const double& zetamin);
  145 + void vetorate_save(const GFilename& filename) const;
  146 + void vetorate_load(const GFilename& filename);
101 147  
102 148 // Protected data members
103 149 std::vector<GCOMDri> m_dris; //!< Data space instances
  150 +
  151 + // Working arrays for vetorate computation
  152 + GNdarray m_wrk_counts; //!< 3D event cube array
  153 + GNdarray m_wrk_ehacutcorr; //!< 2D geometry response array
  154 + GNdarray m_wrk_vetorate; //!< 1D vetorates array
  155 + GNdarray m_wrk_constrate; //!< 1D constant rate array
  156 + GNdarray m_wrk_rate; //!< 2D rate array
  157 + std::vector<bool> m_wrk_use_sp; //!< 1D superpacket usage array
104 158 };
105 159  
106 160  
... ... @@ -187,4 +241,46 @@ void GCOMDris::reserve(const int&amp; num)
187 241 return;
188 242 }
189 243  
  244 +
  245 +/***********************************************************************//**
  246 + * @brief Return log-likelihood value of optimizer function
  247 + *
  248 + * @return Log-likelihood value of optimizer function.
  249 + *
  250 + * Returns the log-likelihood value of optimizer function.
  251 + ***************************************************************************/
  252 +inline
  253 +double GCOMDris::likelihood::value(void) const
  254 +{
  255 + return m_value;
  256 +}
  257 +
  258 +
  259 +/***********************************************************************//**
  260 + * @brief Return pointer to gradient vector
  261 + *
  262 + * @return Pointer to gradient vector.
  263 + *
  264 + * Returns a pointer to the parameter gradient vector.
  265 + ***************************************************************************/
  266 +inline
  267 +GVector* GCOMDris::likelihood::gradient(void)
  268 +{
  269 + return &(m_gradient);
  270 +}
  271 +
  272 +
  273 +/***********************************************************************//**
  274 + * @brief Return pointer to curvature matrix
  275 + *
  276 + * @return Pointer to curvature matrix.
  277 + *
  278 + * Returns a pointer to the parameter curvature matrix.
  279 + ***************************************************************************/
  280 +inline
  281 +GMatrixSparse* GCOMDris::likelihood::curvature(void)
  282 +{
  283 + return &(m_curvature);
  284 +}
  285 +
190 286 #endif /* GCOMDRIS_HPP */
... ...
inst/com/pyext/GCOMDris.i
  View file @1fc91b4
... ... @@ -57,7 +57,7 @@ public:
57 57 const GCOMSelection& select = GCOMSelection(),
58 58 const double& zetamin = 5.0,
59 59 const double& timebin = 300.0,
60   - const std::string& method = "phibar");
  60 + const std::string& method = "vetorate");
61 61 };
62 62  
63 63  
... ...
inst/com/src/GCOMDris.cpp
  View file @1fc91b4
... ... @@ -32,6 +32,12 @@
32 32 #include <iostream>
33 33 #include "GException.hpp"
34 34 #include "GMath.hpp"
  35 +#include "GFits.hpp"
  36 +#include "GFitsImageDouble.hpp"
  37 +#include "GFitsImageShort.hpp"
  38 +#include "GOptimizerPars.hpp"
  39 +#include "GOptimizerPar.hpp"
  40 +#include "GOptimizerLM.hpp"
35 41 #include "GCOMTools.hpp"
36 42 #include "GCOMSupport.hpp"
37 43 #include "GCOMDri.hpp"
... ... @@ -58,7 +64,8 @@
58 64 #define G_COMPUTE_DRWS_EHACORR //!< Correct Earth horizon cut
59 65  
60 66 /* __ Debug definitions __________________________________________________ */
61   -//#define G_DEBUG_COMPUTE_DRWS
  67 +#define G_DEBUG_COMPUTE_DRWS
  68 +#define G_DEBUG_SAVE_WORKING_ARRAYS
62 69  
63 70 /* __ Constants __________________________________________________________ */
64 71  
... ... @@ -351,7 +358,7 @@ void GCOMDris::extend(const GCOMDris&amp; dris)
351 358 * @param[in] select Selection set.
352 359 * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
353 360 * @param[in] timebin Time binning for rate determination (sec).
354   - * @param[in] method Method to compute DRWs ("energy" or "phibar").
  361 + * @param[in] method Method to compute DRWs ("energy", "phibar" or "vetorate").
355 362 *
356 363 * @exception GException::invalid_value
357 364 * DRW with mismatching definition encountered in container.
... ... @@ -407,6 +414,9 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
407 414 // Set all DRG bins to zero
408 415 m_dris[i].init_cube();
409 416  
  417 + // Set DRW header method parameter
  418 + m_dris[i].m_drw_method = method;
  419 +
410 420 // Debug
411 421 #if defined(G_DEBUG_COMPUTE_DRWS)
412 422 std::cout << m_dris[i].print() << std::endl;
... ... @@ -433,6 +443,11 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
433 443 compute_drws_phibar(obs, events, select, zetamin, timebin);
434 444 }
435 445  
  446 + // Method C: Veto rate weighting
  447 + else if (method == "vetorate") {
  448 + compute_drws_vetorate(obs, events, select, zetamin);
  449 + }
  450 +
436 451 // Normalise Phibar layers of all DRWs to unity
437 452 for (int i = 0; i < size(); ++i) {
438 453  
... ... @@ -512,6 +527,14 @@ void GCOMDris::init_members(void)
512 527 // Initialise members
513 528 m_dris.clear();
514 529  
  530 + // Initialise working arrays
  531 + m_wrk_counts.clear();
  532 + m_wrk_ehacutcorr.clear();
  533 + m_wrk_vetorate.clear();
  534 + m_wrk_constrate.clear();
  535 + m_wrk_rate.clear();
  536 + m_wrk_use_sp.clear();
  537 +
515 538 // Return
516 539 return;
517 540 }
... ... @@ -527,6 +550,14 @@ void GCOMDris::copy_members(const GCOMDris&amp; dris)
527 550 // Copy members
528 551 m_dris = dris.m_dris;
529 552  
  553 + // Copy working arrays
  554 + m_wrk_counts = dris.m_wrk_counts;
  555 + m_wrk_ehacutcorr = dris.m_wrk_ehacutcorr;
  556 + m_wrk_vetorate = dris.m_wrk_vetorate;
  557 + m_wrk_constrate = dris.m_wrk_constrate;
  558 + m_wrk_rate = dris.m_wrk_rate;
  559 + m_wrk_use_sp = dris.m_wrk_use_sp;
  560 +
530 561 // Return
531 562 return;
532 563 }
... ... @@ -546,7 +577,7 @@ void GCOMDris::free_members(void)
546 577 * @brief Compute background weighting cubes using energy dependent rates
547 578 *
548 579 * @param[in] obs COMPTEL observation.
549   - * @param[in] event COMPTEL event list.
  580 + * @param[in] events COMPTEL event list.
550 581 * @param[in] select Selection set.
551 582 * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
552 583 * @param[in] timebin Time binning for rate determination (sec).
... ... @@ -869,7 +900,7 @@ void GCOMDris::compute_drws_energy(const GCOMObservation&amp; obs,
869 900 * @brief Compute background weighting cubes using Phibar dependent rates
870 901 *
871 902 * @param[in] obs COMPTEL observation.
872   - * @param[in] event COMPTEL event list.
  903 + * @param[in] events COMPTEL event list.
873 904 * @param[in] select Selection set.
874 905 * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
875 906 * @param[in] timebin Time binning for rate determination (sec).
... ... @@ -1199,3 +1230,881 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1199 1230 // Return
1200 1231 return;
1201 1232 }
  1233 +
  1234 +
  1235 +/***********************************************************************//**
  1236 + * @brief Compute background weighting cubes using veto rates
  1237 + *
  1238 + * @param[in] obs COMPTEL observation.
  1239 + * @param[in] events COMPTEL event list.
  1240 + * @param[in] select Selection set.
  1241 + * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  1242 + *
  1243 + * Compute DRW cubes for a COMPTEL observation. DRW cubes are event-rate
  1244 + * weighted geometry cubes which were multiplied by the solid angle of
  1245 + * the (Chi,Psi) pixels.
  1246 + ***************************************************************************/
  1247 +void GCOMDris::compute_drws_vetorate(const GCOMObservation& obs,
  1248 + const GCOMEventList* events,
  1249 + const GCOMSelection& select,
  1250 + const double& zetamin)
  1251 +{
  1252 + // Debug
  1253 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1254 + std::cout << "GCOMDris::compute_drws_vetorate" << std::endl;
  1255 + std::cout << "===============================" << std::endl;
  1256 + #endif
  1257 +
  1258 + // Setup working arrays
  1259 + vetorate_setup(obs, events, select, zetamin);
  1260 +
  1261 + // Debug: Save working arrays in FITS file
  1262 + #if defined(G_DEBUG_SAVE_WORKING_ARRAYS)
  1263 + vetorate_save("debug_gcomdris_compute_drws_vetorate.fits");
  1264 + vetorate_load("debug_gcomdris_compute_drws_vetorate.fits");
  1265 + #endif
  1266 +
  1267 + // Fit working arrays to determine f_prompt parameter
  1268 + vetorate_fit();
  1269 +
  1270 + // Generate DRWs
  1271 + vetorate_generate(obs, select, zetamin);
  1272 +
  1273 + // Return
  1274 + return;
  1275 +}
  1276 +
  1277 +
  1278 +/***********************************************************************//**
  1279 + * @brief Setup working arrays for vetorate computation
  1280 + *
  1281 + * @param[in] obs COMPTEL observation.
  1282 + * @param[in] events COMPTEL event list.
  1283 + * @param[in] select Selection set.
  1284 + * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  1285 + *
  1286 + * Setup working arrays for vetorate DRW computation. There are four working
  1287 + * arrays:
  1288 + *
  1289 + * @c m_wrk_counts is a 3D working array spanned by superpacket index,
  1290 + * phibar layer and energy bin, containing the number of events passing the
  1291 + * event selections for the specified bin.
  1292 + *
  1293 + * @c m_wrk_ehacutcorr is a 2D working array spanned by superpacket index
  1294 + * and phibar layer, containing the fraction of the solid-angle weighted
  1295 + * geometry that passes the Earth horizon cut for the specified bin. This
  1296 + * array represents the response to a given incident event rate.
  1297 + *
  1298 + * @c m_wrk_vetorate is a 1D working array, containing the veto rate as
  1299 + * function of superpacket index.
  1300 + *
  1301 + * @c m_wrk_use_sp is a 1D working array, signalling the usage of a given
  1302 + * superpacket.
  1303 + *
  1304 + * Note that the method also sets the superpacket usage and event selection
  1305 + * statistics of all DRWs.
  1306 + ***************************************************************************/
  1307 +void GCOMDris::vetorate_setup(const GCOMObservation& obs,
  1308 + const GCOMEventList* events,
  1309 + const GCOMSelection& select,
  1310 + const double& zetamin)
  1311 +{
  1312 + // Debug
  1313 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1314 + std::cout << "GCOMDris::vetorate_setup" << std::endl;
  1315 + std::cout << "========================" << std::endl;
  1316 + #endif
  1317 +
  1318 + // Initialise D1 & D2 module status
  1319 + GCOMStatus status;
  1320 +
  1321 + // Get pointer to first DRW (assume their definition is identical)
  1322 + const GCOMDri* dri0 = &(m_dris[0]);
  1323 +
  1324 + // Extract dimensions
  1325 + int noads = obs.oads().size();
  1326 + int npix = dri0->m_dri.npix();
  1327 + int nphibar = dri0->nphibar();
  1328 + int neng = size();
  1329 +
  1330 + // Debug
  1331 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1332 + std::cout << "Number of superpackets: " << noads << std::endl;
  1333 + std::cout << "Number of pixels: " << npix << std::endl;
  1334 + std::cout << "Number of Phibar layers: " << nphibar << std::endl;
  1335 + std::cout << "Number of energies: " << neng << std::endl;
  1336 + #endif
  1337 +
  1338 + // Define working arrays
  1339 + m_wrk_counts = GNdarray(noads, nphibar, neng);
  1340 + m_wrk_ehacutcorr = GNdarray(noads, nphibar);
  1341 + m_wrk_vetorate = GNdarray(noads);
  1342 + m_wrk_use_sp = std::vector<bool>(noads);
  1343 +
  1344 + // Allocate further working arrays
  1345 + GNdarray geo(npix);
  1346 + GNdarray eha(npix);
  1347 +
  1348 + // Setup node array and rate vector for veto rate interpolation
  1349 + const GCOMHkd& hdk = obs.hkds()["SCV2M"];
  1350 + GNodeArray veto_times;
  1351 + std::vector<double> veto_rates;
  1352 + veto_times.reserve(hdk.size());
  1353 + veto_rates.reserve(hdk.size());
  1354 + for (int i = 0; i < hdk.size(); ++i) {
  1355 + double rate = hdk.value(i);
  1356 + if (rate > 0.0) {
  1357 + veto_times.append(hdk.time(i).secs());
  1358 + veto_rates.push_back(rate);
  1359 + }
  1360 + }
  1361 +
  1362 + // Setup sky directions and solid angles for all (Chi, Psi)
  1363 + std::vector<GSkyDir> skys;
  1364 + std::vector<double> omegas;
  1365 + for (int index = 0; index < dri0->m_dri.npix(); ++index) {
  1366 + GSkyDir sky = dri0->m_dri.inx2dir(index);
  1367 + double omega = dri0->m_dri.solidangle(index);
  1368 + skys.push_back(sky);
  1369 + omegas.push_back(omega);
  1370 + }
  1371 +
  1372 + // Debug
  1373 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1374 + std::cout << "Setup node array for veto rate interpolation: ";
  1375 + std::cout << veto_times.size() << " times, ";
  1376 + std::cout << veto_rates.size() << " rates" << std::endl;
  1377 + #endif
  1378 +
  1379 + // Get Good Time Intervals. If a pulsar selection is specified then
  1380 + // reduce the Good Time Intervals to the validity intervals of the
  1381 + // pulsar emphemerides.
  1382 + GCOMTim tim = obs.tim();
  1383 + if (select.has_pulsar()) {
  1384 + tim.reduce(select.pulsar().validity());
  1385 + }
  1386 +
  1387 + // Initialise event counter
  1388 + int i_evt = 0;
  1389 +
  1390 + // Signal that loop should be terminated
  1391 + bool terminate = false;
  1392 +
  1393 + // Loop over Orbit Aspect Data
  1394 + for (int i_oad = 0; i_oad < noads; ++i_oad) {
  1395 +
  1396 + // Get reference to Orbit Aspect Data of superpacket
  1397 + const GCOMOad &oad = obs.oads()[i_oad];
  1398 +
  1399 + // Check superpacket usage for all DRWs so that their statistics
  1400 + // get correctly updated
  1401 + m_wrk_use_sp[i_oad] = true;
  1402 + 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;
  1405 + }
  1406 + }
  1407 + if (!m_wrk_use_sp[i_oad]) {
  1408 + continue;
  1409 + }
  1410 +
  1411 + // Get Orbit Aspect Data mid-time in seconds
  1412 + double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
  1413 +
  1414 + // Prepare Earth horizon angle computation. The celestial system
  1415 + // is reinterpreted as the COMPTEL coordinate system, where the
  1416 + // 90 degrees - zenith angle becomes the declination and the azimuth
  1417 + // angle becomes Right Ascension. This allows us later to use the
  1418 + // GSkyDir::dist method to compute the distance between the geocentre
  1419 + // and the telescope Z-axis.
  1420 + double theta_geocentre = double(oad.gcel());
  1421 + double phi_geocentre = double(oad.gcaz());
  1422 + GSkyDir geocentre_comptel;
  1423 + geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
  1424 +
  1425 + // Compute solid-angle-corrected geometry factors and Earth horizon
  1426 + // angles for all (Chi,Psi) pixels
  1427 + double total_geo = 0.0;
  1428 + for (int index = 0; index < npix; ++index) {
  1429 +
  1430 + // Convert sky direction to COMPTEL coordinates
  1431 + double theta = oad.theta(skys[index]);
  1432 + double phi = oad.phi(skys[index]);
  1433 +
  1434 + // Compute geometric factor summed over all D1, D2 times
  1435 + // the solid angle of the weighting cube pixel
  1436 + double geometry = dri0->compute_geometry(oad.tjd(), theta, phi,
  1437 + select, status) * omegas[index];
  1438 +
  1439 + // Store and sum up geometric factor
  1440 + geo(index) = geometry;
  1441 + total_geo += geometry;
  1442 +
  1443 + // Compute Earth horizon angle as the distance between the Earth
  1444 + // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
  1445 + // COMPTEL coordinates minus the radius of the Earth.
  1446 + GSkyDir chipsi_comptel;
  1447 + chipsi_comptel.radec_deg(phi, 90.0-theta);
  1448 + eha(index) = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
  1449 +
  1450 + } // endfor: looped over (Chi, Psi)
  1451 +
  1452 + // Compute Earth horizon cut correction factor for all Phibar layers
  1453 + if (total_geo > 0.0) {
  1454 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1455 + double ehamin = double(iphibar) * dri0->m_phibin + zetamin;
  1456 + double accepted_geo = 0.0;
  1457 + for (int index = 0; index < npix; ++index) {
  1458 + if (eha(index) >= ehamin) {
  1459 + accepted_geo += geo(index);
  1460 + }
  1461 + }
  1462 + m_wrk_ehacutcorr(i_oad, iphibar) = accepted_geo / total_geo;
  1463 + }
  1464 + }
  1465 +
  1466 + // Compute veto rate for superpacket time if rates are available
  1467 + if (!veto_times.is_empty()) {
  1468 + double value = veto_times.interpolate(time, veto_rates);
  1469 + if (value > 0.0) {
  1470 + m_wrk_vetorate(i_oad) = value;
  1471 + }
  1472 + }
  1473 +
  1474 + // Collect all events within superpacket and store result in 3D array.
  1475 + // Break if the end of the event list was reached.
  1476 + for (; i_evt < events->size(); ++i_evt) {
  1477 +
  1478 + // Get pointer to event
  1479 + const GCOMEventAtom* event = (*events)[i_evt];
  1480 +
  1481 + // Break loop if the end of the superpacket was reached
  1482 + if (event->time() > oad.tstop()) {
  1483 + break;
  1484 + }
  1485 +
  1486 + // Determine energy bin
  1487 + int ienergy = -1;
  1488 + for (int i = 0; i < size(); ++i) {
  1489 + if ((event->energy() >= m_dris[i].m_ebounds.emin()) &&
  1490 + (event->energy() <= m_dris[i].m_ebounds.emax())) {
  1491 + ienergy = i;
  1492 + break;
  1493 + }
  1494 + }
  1495 +
  1496 + // Skip event if it is not contained in any energy bin
  1497 + if (ienergy == -1) {
  1498 + continue;
  1499 + }
  1500 +
  1501 + // Get pointer to appropriate DRW
  1502 + const GCOMDri* dri = &(m_dris[ienergy]);
  1503 +
  1504 + // Check GTIs if the DRW has GTIs
  1505 + if (dri->gti().size() > 0) {
  1506 +
  1507 + // Skip event if it lies before the GTI start
  1508 + if (event->time() < dri->gti().tstart()) {
  1509 + continue;
  1510 + }
  1511 +
  1512 + // Break Orbit Aspect Data loop if event lies after the GTI stop
  1513 + else if (event->time() > dri->gti().tstop()) {
  1514 + terminate = true;
  1515 + break;
  1516 + }
  1517 +
  1518 + } // endif: DRW had GTIs
  1519 +
  1520 + // Skip event if it lies before the superpacket start
  1521 + if (event->time() < oad.tstart()) {
  1522 + continue;
  1523 + }
  1524 +
  1525 + // Apply event selection
  1526 + if (!dri->m_selection.use_event(*event)) {
  1527 + continue;
  1528 + }
  1529 +
  1530 + // Compute Compton scatter angle index. Skip if it's invalid.
  1531 + int iphibar = int((event->phibar() - dri->m_phimin) / dri->m_phibin);
  1532 + if (iphibar < 0) {
  1533 + continue;
  1534 + }
  1535 + else if (iphibar >= dri->nphibar()) {
  1536 + continue;
  1537 + }
  1538 +
  1539 + // Check for Earth horizon angle. There is a constant EHA limit
  1540 + // over a Phibar layer to be compliant with the DRG.
  1541 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  1542 + if (event->eha() < ehamin) {
  1543 + continue;
  1544 + }
  1545 +
  1546 + // Now fill the event in the 3D counts array
  1547 + m_wrk_counts(i_oad, iphibar, ienergy) += 1.0;
  1548 +
  1549 + } // endfor: collected events
  1550 +
  1551 + // Debug
  1552 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1553 + std::cout << "Superpacket " << i_oad;
  1554 + std::cout << " time=" << time;
  1555 + std::cout << " vetorate=" << m_wrk_vetorate(i_oad);
  1556 + std::cout << " total_geo=" << total_geo;
  1557 + std::cout << std::endl;
  1558 + #endif
  1559 +
  1560 + // Break Orbit Aspect Data loop if termination was signalled or if there
  1561 + // are no more events
  1562 + if (terminate || i_evt >= events->size()) {
  1563 + break;
  1564 + }
  1565 +
  1566 + } // endfor: looped over Orbit Aspect Data
  1567 +
  1568 + // Return
  1569 + return;
  1570 +}
  1571 +
  1572 +
  1573 +/***********************************************************************//**
  1574 + * @brief Fit working arrays for vetorate computation
  1575 + *
  1576 + * Determines the f_prompt parameter for each energy bin using a maximum
  1577 + * log-likelihood fit using the data in the working arrays for each energy
  1578 + * bin. Based on the fitted f_prompt parameter the method also computes the
  1579 + * expected background rate time variation for each energy bin.
  1580 + ***************************************************************************/
  1581 +void GCOMDris::vetorate_fit(void)
  1582 +{
  1583 + // Set constants
  1584 + const double fprompt_init = 0.5; // Initial f_prompt value
  1585 + const double norm = 1.0; // Model normalisation
  1586 +
  1587 + // Debug
  1588 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1589 + std::cout << "GCOMDris::vetorate_fit" << std::endl;
  1590 + std::cout << "======================" << std::endl;
  1591 + #endif
  1592 +
  1593 + // Extract dimension from working arrays
  1594 + int noads = m_wrk_counts.shape()[0];
  1595 + int nphibar = m_wrk_counts.shape()[1];
  1596 + int neng = m_wrk_counts.shape()[2];
  1597 +
  1598 + // Debug
  1599 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1600 + std::cout << "Number of superpackets: " << noads << std::endl;
  1601 + std::cout << "Number of Phibar layers: " << nphibar << std::endl;
  1602 + std::cout << "Number of energies: " << neng << std::endl;
  1603 + #endif
  1604 +
  1605 + // Setup rate result array
  1606 + m_wrk_rate = GNdarray(noads, neng);
  1607 +
  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);
  1617 + }
  1618 + }
  1619 + if (nvetorate > 0.0) {
  1620 + for (int ioad = 0; ioad < noads; ++ioad) {
  1621 + m_wrk_vetorate(ioad) /= nvetorate;
  1622 + }
  1623 + }
  1624 + if (nconstrate > 0.0) {
  1625 + for (int ioad = 0; ioad < noads; ++ioad) {
  1626 + m_wrk_constrate(ioad) /= nconstrate;
  1627 + }
  1628 + }
  1629 +
  1630 + // Loop over energy bins
  1631 + for (int ieng = 0; ieng < neng; ++ieng) {
  1632 +
  1633 + // Set fprompt parameter
  1634 + GOptimizerPars pars;
  1635 + GOptimizerPar par("fprompt", fprompt_init);
  1636 + par.range(0.0, 1.0);
  1637 + par.factor_gradient(1.0); // To avoid zero gradient log message
  1638 + pars.append(par);
  1639 +
  1640 + // Set fit function
  1641 + GCOMDris::likelihood fct(this, ieng, norm);
  1642 +
  1643 + // Optimize function and compute errors
  1644 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1645 + GLog log;
  1646 + log.cout(true);
  1647 + GOptimizerLM opt(&log);
  1648 + #else
  1649 + GOptimizerLM opt;
  1650 + #endif
  1651 + opt.eps(0.1);
  1652 + opt.max_iter(500);
  1653 + opt.optimize(fct, pars);
  1654 + opt.errors(fct, pars);
  1655 +
  1656 + // Retrieve logL
  1657 + double logL = opt.value();
  1658 +
  1659 + // Recover fprompt parameter and errors
  1660 + double fprompt = pars[0]->value();
  1661 + double e_fprompt = pars[0]->error();
  1662 + double fconst = 1.0 - fprompt;
  1663 +
  1664 + // 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);
  1668 + }
  1669 +
  1670 + // Set DRW header parameters
  1671 + m_dris[ieng].m_drw_status = opt.status_string();
  1672 + m_dris[ieng].m_drw_fprompt = fprompt;
  1673 + m_dris[ieng].m_drw_e_fprompt = e_fprompt;
  1674 + m_dris[ieng].m_drw_iter = opt.iter();
  1675 +
  1676 + // Debug: Print fit results
  1677 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1678 + std::cout << ieng;
  1679 + std::cout << " logL=" << logL;
  1680 + std::cout << " fprompt=" << fprompt << " +/- " << e_fprompt;
  1681 + std::cout << std::endl;
  1682 + #endif
  1683 +
  1684 + } // endfor: looped over energy bins
  1685 +
  1686 + // Return
  1687 + return;
  1688 +}
  1689 +
  1690 +
  1691 +/***********************************************************************//**
  1692 + * @brief Generate DRWs
  1693 + *
  1694 + * @param[in] obs COMPTEL observation.
  1695 + * @param[in] select Selection set.
  1696 + * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  1697 + ***************************************************************************/
  1698 +void GCOMDris::vetorate_generate(const GCOMObservation& obs,
  1699 + const GCOMSelection& select,
  1700 + const double& zetamin)
  1701 +{
  1702 + // Debug
  1703 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1704 + std::cout << "GCOMDris::vetorate_generate" << std::endl;
  1705 + std::cout << "===========================" << std::endl;
  1706 + #endif
  1707 +
  1708 + // Initialise D1 & D2 module status
  1709 + GCOMStatus status;
  1710 +
  1711 + // Get pointer to first DRW (assume their definition is identical)
  1712 + const GCOMDri* dri0 = &(m_dris[0]);
  1713 +
  1714 + // Extract dimensions
  1715 + int noads = obs.oads().size();
  1716 + int npix = dri0->m_dri.npix();
  1717 + int nphibar = dri0->nphibar();
  1718 + int neng = size();
  1719 +
  1720 + // Debug
  1721 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1722 + std::cout << "Number of superpackets: " << noads << std::endl;
  1723 + std::cout << "Number of pixels: " << npix << std::endl;
  1724 + std::cout << "Number of Phibar layers: " << nphibar << std::endl;
  1725 + std::cout << "Number of energies: " << neng << std::endl;
  1726 + #endif
  1727 +
  1728 + // Allocate geometry factor and Earth horizon angle working arrays
  1729 + GNdarray geometry(npix);
  1730 + GNdarray eha(npix);
  1731 +
  1732 + // Setup sky directions and solid angles for all (Chi, Psi)
  1733 + std::vector<GSkyDir> skys;
  1734 + std::vector<double> omegas;
  1735 + for (int index = 0; index < npix; ++index) {
  1736 + GSkyDir sky = dri0->m_dri.inx2dir(index);
  1737 + double omega = dri0->m_dri.solidangle(index);
  1738 + skys.push_back(sky);
  1739 + omegas.push_back(omega);
  1740 + }
  1741 +
  1742 + // Get Good Time Intervals. If a pulsar selection is specified then
  1743 + // reduce the Good Time Intervals to the validity intervals of the
  1744 + // pulsar emphemerides.
  1745 + GCOMTim tim = obs.tim();
  1746 + if (select.has_pulsar()) {
  1747 + tim.reduce(select.pulsar().validity());
  1748 + }
  1749 +
  1750 + // Loop over Orbit Aspect Data
  1751 + for (int i_oad = 0; i_oad < noads; ++i_oad) {
  1752 +
  1753 + // Get reference to Orbit Aspect Data of superpacket
  1754 + const GCOMOad &oad = obs.oads()[i_oad];
  1755 +
  1756 + // Skip not-to-be-used superpackets
  1757 + if (!m_wrk_use_sp[i_oad]) {
  1758 + continue;
  1759 + }
  1760 +
  1761 + // Prepare Earth horizon angle computation. The celestial system
  1762 + // is reinterpreted as the COMPTEL coordinate system, where the
  1763 + // 90 degrees - zenith angle becomes the declination and the azimuth
  1764 + // angle becomes Right Ascension. This allows us later to use the
  1765 + // GSkyDir::dist method to compute the distance between the geocentre
  1766 + // and the telescope Z-axis.
  1767 + double theta_geocentre = double(oad.gcel());
  1768 + double phi_geocentre = double(oad.gcaz());
  1769 + GSkyDir geocentre_comptel;
  1770 + geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
  1771 +
  1772 + // Compute solid-angle-corrected geometry factors and Earth horizon
  1773 + // angles for all (Chi,Psi) pixels
  1774 + for (int index = 0; index < npix; ++index) {
  1775 +
  1776 + // Convert sky direction to COMPTEL coordinates
  1777 + double theta = oad.theta(skys[index]);
  1778 + double phi = oad.phi(skys[index]);
  1779 +
  1780 + // Compute geometric factor summed over all D1, D2 times
  1781 + // the solid angle of the weighting cube pixel
  1782 + geometry(index) = dri0->compute_geometry(oad.tjd(), theta, phi,
  1783 + select, status) * omegas[index];
  1784 +
  1785 + // Compute Earth horizon angle as the distance between the Earth
  1786 + // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
  1787 + // COMPTEL coordinates minus the radius of the Earth.
  1788 + GSkyDir chipsi_comptel;
  1789 + chipsi_comptel.radec_deg(phi, 90.0-theta);
  1790 + eha(index) = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
  1791 +
  1792 + } // endfor: computed solid-angle-corrected geometry factors
  1793 +
  1794 + // Loop over all DRWs
  1795 + for (int i = 0; i < size(); ++i) {
  1796 +
  1797 + // Get pointer to appropriate DRW
  1798 + const GCOMDri* dri = &(m_dris[i]);
  1799 +
  1800 + // Loop over all Phibar layers
  1801 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1802 +
  1803 + // Compute minimum Earth horizon angle for Phibar layer
  1804 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  1805 +
  1806 + // Loop over all (Chi,Psi) pixels
  1807 + for (int index = 0; index < npix; ++index) {
  1808 +
  1809 + // If Earth horizon angle is equal or larger than the minimum
  1810 + // requirement then add up the rate weighted geometry factor
  1811 + if (eha(index) >= ehamin) {
  1812 +
  1813 + // Get data space index
  1814 + int inx = index + iphibar * npix;
  1815 +
  1816 + // Store geometry weighted with rate as DRW
  1817 + m_dris[i][inx] += geometry(index) * m_wrk_rate(i_oad, i);
  1818 +
  1819 + } // endif: Earth horizon cut passed
  1820 +
  1821 + } // endfor: looped over (Chi, Psi)
  1822 +
  1823 + } // endfor: looped over Phibar
  1824 +
  1825 + } // endfor: looped over all DRWs
  1826 +
  1827 + // Debug
  1828 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1829 + std::cout << "Superpacket " << i_oad << std::endl;
  1830 + #endif
  1831 +
  1832 + } // endfor: looped over Orbit Aspect Data
  1833 +
  1834 + // Return
  1835 + return;
  1836 +}
  1837 +
  1838 +
  1839 +/***********************************************************************//**
  1840 + * @brief Save working arrays for vetorate computation
  1841 + *
  1842 + * @param[in] filename FITS filename.
  1843 + *
  1844 + * Save working arrays for vetorate DRW computation in FITS file.
  1845 + ***************************************************************************/
  1846 +void GCOMDris::vetorate_save(const GFilename& filename) const
  1847 +{
  1848 + // Debug
  1849 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1850 + std::cout << "GCOMDris::vetorate_save" << std::endl;
  1851 + std::cout << "=======================" << std::endl;
  1852 + #endif
  1853 +
  1854 + // Extract dimension from working arrays
  1855 + int noads = m_wrk_counts.shape()[0];
  1856 + int nphibar = m_wrk_counts.shape()[1];
  1857 + int neng = m_wrk_counts.shape()[2];
  1858 +
  1859 + // Allocate FITS images
  1860 + GFitsImageDouble image_counts(noads, nphibar, neng);
  1861 + GFitsImageDouble image_ehacutcorr(noads, nphibar);
  1862 + GFitsImageDouble image_vetorate(noads);
  1863 + GFitsImageShort image_use_sp(noads);
  1864 +
  1865 + // Fill FITS images from working arrays
  1866 + for (int ioad = 0; ioad < noads; ++ioad) {
  1867 + 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);
  1870 + }
  1871 + image_ehacutcorr(ioad, iphibar) = m_wrk_ehacutcorr(ioad, iphibar);
  1872 + }
  1873 + image_vetorate(ioad) = m_wrk_vetorate(ioad);
  1874 + image_use_sp(ioad) = (m_wrk_use_sp[ioad]) ? 1 : 0;
  1875 + }
  1876 +
  1877 + // Set image attributes
  1878 + image_counts.extname("COUNTS");
  1879 + image_ehacutcorr.extname("EHACUTCORR");
  1880 + image_vetorate.extname("VETORATE");
  1881 + image_use_sp.extname("SPUSAGE");
  1882 +
  1883 + // Create FITS file and append images
  1884 + GFits fits;
  1885 + fits.append(image_counts);
  1886 + fits.append(image_ehacutcorr);
  1887 + fits.append(image_vetorate);
  1888 + fits.append(image_use_sp);
  1889 +
  1890 + // Save FITS file
  1891 + fits.saveto(filename, true);
  1892 +
  1893 + // Return
  1894 + return;
  1895 +}
  1896 +
  1897 +
  1898 +/***********************************************************************//**
  1899 + * @brief Load working arrays for vetorate computation
  1900 + *
  1901 + * @param[in] filename FITS filename.
  1902 + *
  1903 + * Load working arrays for vetorate DRW computation from FITS file.
  1904 + ***************************************************************************/
  1905 +void GCOMDris::vetorate_load(const GFilename& filename)
  1906 +{
  1907 + // Debug
  1908 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1909 + std::cout << "GCOMDris::vetorate_load" << std::endl;
  1910 + std::cout << "=======================" << std::endl;
  1911 + #endif
  1912 +
  1913 + // Open FITS file
  1914 + GFits fits(filename);
  1915 +
  1916 + // Get images
  1917 + const GFitsImage* image_counts = fits.image("COUNTS");
  1918 + const GFitsImage* image_ehacutcorr = fits.image("EHACUTCORR");
  1919 + const GFitsImage* image_vetorate = fits.image("VETORATE");
  1920 + const GFitsImage* image_use_sp = fits.image("SPUSAGE");
  1921 +
  1922 + // Extract dimensions
  1923 + int noads = image_counts->naxes(0);
  1924 + int nphibar = image_counts->naxes(1);
  1925 + int neng = image_counts->naxes(2);
  1926 +
  1927 + // Debug
  1928 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1929 + std::cout << "Number of superpackets: " << noads << std::endl;
  1930 + std::cout << "Number of Phibar layers: " << nphibar << std::endl;
  1931 + std::cout << "Number of energies: " << neng << std::endl;
  1932 + #endif
  1933 +
  1934 + // Allocate working array
  1935 + m_wrk_counts = GNdarray(noads, nphibar, neng);
  1936 + m_wrk_ehacutcorr = GNdarray(noads, nphibar);
  1937 + m_wrk_vetorate = GNdarray(noads);
  1938 + m_wrk_use_sp = std::vector<bool>(noads);
  1939 +
  1940 + // Extract data from images
  1941 + for (int ioad = 0; ioad < noads; ++ioad) {
  1942 + 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);
  1945 + }
  1946 + m_wrk_ehacutcorr(ioad, iphibar) = image_ehacutcorr->pixel(ioad, iphibar);
  1947 + }
  1948 + m_wrk_vetorate(ioad) = image_vetorate->pixel(ioad);
  1949 + m_wrk_use_sp[ioad] = (image_use_sp->pixel(ioad) == 1);
  1950 + }
  1951 +
  1952 + // Close FITS file
  1953 + fits.close();
  1954 +
  1955 + // Return
  1956 + return;
  1957 +}
  1958 +
  1959 +
  1960 +/***********************************************************************//**
  1961 + * @brief Log-likelihood function constructor
  1962 + *
  1963 + * @param[in] dris Parent calling likelihood function.
  1964 + * @param[in] ieng DRW energy bin.
  1965 + * @param[in] norm Normalisation constant.
  1966 + *
  1967 + * Constructs the log-likelihood function that is used to determine the
  1968 + * value of f_prompt.
  1969 + ***************************************************************************/
  1970 +GCOMDris::likelihood::likelihood(GCOMDris *dris,
  1971 + const int& ieng,
  1972 + const double& norm) : GOptimizerFunction()
  1973 +{
  1974 + // Store arguments
  1975 + m_this = dris;
  1976 + m_ieng = ieng;
  1977 + m_norm = norm;
  1978 +
  1979 + // Set gradient and curvature members
  1980 + m_gradient = GVector(1);
  1981 + m_curvature = GMatrixSparse(1,1);
  1982 +
  1983 + // Extract dimension from working array
  1984 + m_noads = m_this->m_wrk_counts.shape()[0];
  1985 + m_nphibar = m_this->m_wrk_counts.shape()[1];
  1986 +
  1987 + // Multiply veto and constant rate by EHA cut correction. Rates are zero
  1988 + // 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);
  1994 + if (vetorate > 0.0) {
  1995 + double constrate = m_this->m_wrk_constrate(ioad);
  1996 + 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);
  2000 + }
  2001 + }
  2002 + }
  2003 +
  2004 + // Compute rate sums
  2005 + m_vetorate_sum = GNdarray(m_nphibar);
  2006 + m_constrate_sum = GNdarray(m_nphibar);
  2007 + m_diffrate_sum = GNdarray(m_nphibar);
  2008 + 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);
  2012 + }
  2013 + m_diffrate_sum(iphibar) = m_vetorate_sum(iphibar) - m_constrate_sum(iphibar);
  2014 + }
  2015 +
  2016 + // Return
  2017 + return;
  2018 +}
  2019 +
  2020 +
  2021 +/***********************************************************************//**
  2022 + * @brief Log-likelihood function evaluation
  2023 + *
  2024 + * @param[in] pars Function parameters.
  2025 + *
  2026 + * Computes the log-likelihood function, its gradient and its curvature at
  2027 + * the specified function parameters.
  2028 + ***************************************************************************/
  2029 +void GCOMDris::likelihood::eval(const GOptimizerPars& pars)
  2030 +{
  2031 + // Recover parameters
  2032 + 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];
  2038 +
  2039 + // Allocate phibar normalisation arrays
  2040 + GNdarray nevents(nphibar);
  2041 + GNdarray nmodel(nphibar);
  2042 + GNdarray norm(nphibar);
  2043 +
  2044 + // Compute model
  2045 + GNdarray model = fprompt * m_vetorate + fconst * m_constrate;
  2046 +
  2047 + // Compute Phibar normalised model
  2048 + GNdarray model_norm = model;
  2049 + 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);
  2054 + }
  2055 + }
  2056 + if (nmodel(iphibar) > 0.0) {
  2057 + norm(iphibar) = m_norm * nevents(iphibar) / nmodel(iphibar);
  2058 + for (int ioad = 0; ioad < m_noads; ++ioad) {
  2059 + model_norm(ioad, iphibar) *= norm(iphibar);
  2060 + }
  2061 + }
  2062 + }
  2063 +
  2064 + // Compute LogL
  2065 + m_value = 0.0;
  2066 + for (int ioad = 0; ioad < m_noads; ++ioad) {
  2067 + 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);
  2071 + }
  2072 + }
  2073 + }
  2074 +
  2075 + // Evaluate gradient and curvature
  2076 + for (int iphibar = 0; iphibar < m_nphibar; ++iphibar) {
  2077 +
  2078 + // Precompute normalisation gradient
  2079 + double nmodel2 = nmodel(iphibar) * nmodel(iphibar);
  2080 + double g_norm = -m_norm * nevents(iphibar) / nmodel2 * m_diffrate_sum(iphibar);
  2081 +
  2082 + // Loop over superpackets
  2083 + for (int ioad = 0; ioad < m_noads; ++ioad) {
  2084 +
  2085 + // Continue only if model is positive
  2086 + if (model_norm(ioad, iphibar) > 0.0) {
  2087 +
  2088 + // Pre computation
  2089 + double fb = m_this->m_wrk_counts(ioad, iphibar, m_ieng) / model_norm(ioad, iphibar);
  2090 + double fc = (1.0 - fb);
  2091 + double fa = fb / model_norm(ioad, iphibar);
  2092 +
  2093 + // Compute parameter gradient
  2094 + double g = norm(iphibar) * m_diffrate(ioad, iphibar) + g_norm * model(ioad, iphibar);
  2095 +
  2096 + // Update gradient
  2097 + m_gradient[0] += fc * g;
  2098 +
  2099 + // Update curvature
  2100 + m_curvature(0,0) += fa * g * g;
  2101 +
  2102 + } // endif: model was positive
  2103 +
  2104 + } // endfor: looped over superpackets
  2105 +
  2106 + } // endfor: looped over Phibar
  2107 +
  2108 + // Return
  2109 + return;
  2110 +}
... ...
inst/com/test/dev/test_bgdlixf.py
  View file @1fc91b4
... ... @@ -56,7 +56,6 @@ def get_veto_rates(xmlname):
56 56  
57 57 # Get 2nd Veto dome rates housekeeping data
58 58 hkd = com.hkds()['SCV2M']
59   - print(hkd)
60 59  
61 60 # Extract data
62 61 for i in range(hkd.size()):
... ... @@ -73,35 +72,37 @@ def get_veto_rates(xmlname):
73 72 # ====================================== #
74 73 # Check time ordering of Veto dome rates #
75 74 # ====================================== #
76   -def check_veto_rates(xmlname):
  75 +def check_veto_rates():
77 76 """
78 77 Check time ordering of Veto dome rates
79   -
80   - Parameters
81   - ----------
82   - xmlname : str
83   - Observation definition XML file name
84 78 """
85   - # Load observations
86   - obs = gammalib.GObservations(xmlname)
  79 + # Get all XML files
  80 + xmlnames = glob.glob('dbase/xml/vp*.xml')
  81 + xmlnames.sort()
87 82  
88   - # Loop over observations
89   - for com in obs:
  83 + # Loop over XML files
  84 + for xmlname in xmlnames:
90 85  
91   - # Get 2nd Veto dome rates housekeeping data
92   - hkd = com.hkds()['SCV2M']
  86 + # Load observations
  87 + obs = gammalib.GObservations(xmlname)
93 88  
94   - # Check times
95   - nbad = 0
96   - for i in range(hkd.size()):
97   - if i == 0:
98   - last_time = hkd.time(i)
99   - else:
100   - if hkd.time(i) <= last_time:
101   - nbad += 1
102   - print('%s: %.5f <= %.5f' % (com.id(), hkd.time(i).mjd(), last_time.mjd()))
103   - if nbad == 0:
104   - print('%s: ok.' % (com.id()))
  89 + # Loop over observations
  90 + for com in obs:
  91 +
  92 + # Get 2nd Veto dome rates housekeeping data
  93 + hkd = com.hkds()['SCV2M']
  94 +
  95 + # Check times
  96 + nbad = 0
  97 + for i in range(hkd.size()):
  98 + if i == 0:
  99 + last_time = hkd.time(i)
  100 + else:
  101 + if hkd.time(i) <= last_time:
  102 + nbad += 1
  103 + print('%s: %.5f <= %.5f' % (com.id(), hkd.time(i).mjd(), last_time.mjd()))
  104 + if nbad == 0:
  105 + print('%s: ok.' % (com.id()))
105 106  
106 107 # Return
107 108 return
... ... @@ -147,23 +148,14 @@ if __name__ == &#39;__main__&#39;:
147 148 print('* Test BGDLIXF rate fitting *')
148 149 print('*****************************')
149 150  
150   - # Set file names
151   - #hkdname = '/project-data/comptel/data/phase01/vp0001_0/m20074_hkd.fits'
152   - #hkdname = '/project-data/comptel/data/phase07/vp0701_0/m26689_hkd.fits'
153   - #xmlname = 'dbase/xml/vp0001_0.xml'
154   -
155   - # Get all XML files
156   - xmlnames = glob.glob('dbase/xml/vp*.xml')
157   - xmlnames.sort()
  151 + # Check time ordering of veto rates
  152 + #check_veto_rates()
158 153  
159   - # Loop over XML files
160   - for xmlname in xmlnames:
161   -
162   - # Check Veto rates
163   - check_veto_rates(xmlname)
  154 + # Set file names
  155 + xmlname = 'dbase/xml/vp0001_0.xml'
164 156  
165   - # Get Veto dome rates
166   - #rates = get_veto_rates(xmlname)
  157 + # Get Veto dome rates
  158 + rates = get_veto_rates(xmlname)
167 159  
168   - # Show Veto dome rates
169   - #show_veto_rates(rates)
  160 + # Show Veto dome rates
  161 + show_veto_rates(rates)
... ...