gammalib / gammalib

Commit 9261f32030e3b886359ac0df9e55a051f86d4ccf

Authored by Jürgen Knödlseder
1 parent b0b38ecb

Implement phibar-dependent rate weighting for DRW

Inline Side-by-side
Showing 3 changed files with 638 additions and 239 deletions
inst/com/include/GCOMDris.hpp
  View file @9261f32
... ... @@ -34,6 +34,9 @@
34 34 #include "GContainer.hpp"
35 35  
36 36 /* __ Forward declarations _______________________________________________ */
  37 +class GCOMObservation;
  38 +class GCOMEventList;
  39 +class GCOMSelection;
37 40  
38 41 /* __ Constants __________________________________________________________ */
39 42  
... ... @@ -75,7 +78,9 @@ public:
75 78 void extend(const GCOMDris& dris);
76 79 void compute_drws(const GCOMObservation& obs,
77 80 const GCOMSelection& select = GCOMSelection(),
78   - const double& zetamin = 5.0);
  81 + const double& zetamin = 5.0,
  82 + const double& timebin = 300.0,
  83 + const std::string& method = "phibar");
79 84 std::string print(const GChatter& chatter = NORMAL) const;
80 85  
81 86 protected:
... ... @@ -83,6 +88,16 @@ protected:
83 88 void init_members(void);
84 89 void copy_members(const GCOMDris& dris);
85 90 void free_members(void);
  91 + void compute_drws_energy(const GCOMObservation& obs,
  92 + const GCOMEventList* events,
  93 + const GCOMSelection& select = GCOMSelection(),
  94 + const double& zetamin = 5.0,
  95 + const double& timebin = 300.0);
  96 + void compute_drws_phibar(const GCOMObservation& obs,
  97 + const GCOMEventList* events,
  98 + const GCOMSelection& select = GCOMSelection(),
  99 + const double& zetamin = 5.0,
  100 + const double& timebin = 300.0);
86 101  
87 102 // Protected data members
88 103 std::vector<GCOMDri> m_dris; //!< Data space instances
... ...
inst/com/pyext/GCOMDris.i
  View file @9261f32
... ... @@ -55,7 +55,9 @@ public:
55 55 void extend(const GCOMDris& dris);
56 56 void compute_drws(const GCOMObservation& obs,
57 57 const GCOMSelection& select = GCOMSelection(),
58   - const double& zetamin = 5.0);
  58 + const double& zetamin = 5.0,
  59 + const double& timebin = 300.0,
  60 + const std::string& method = "phibar");
59 61 };
60 62  
61 63  
... ...
inst/com/src/GCOMDris.cpp
  View file @9261f32
... ... @@ -39,14 +39,16 @@
39 39 #include "GCOMTim.hpp"
40 40 #include "GCOMStatus.hpp"
41 41 #include "GCOMObservation.hpp"
  42 +#include "GCOMEventList.hpp"
  43 +#include "GCOMEventAtom.hpp"
42 44 #include "GCOMSelection.hpp"
43 45  
44 46 /* __ Method name definitions ____________________________________________ */
45 47 #define G_AT "GCOMDris::at(int&)"
46 48 #define G_INSERT "GCOMDris::insert(int&, GCOMDri&)"
47 49 #define G_REMOVE "GCOMDris::remove(int&)"
48   -#define G_COMPUTE_DRWS "GCOMDris::compute_drws(GCOMObservation&, "\
49   - "GCOMSelection&, double&)"
  50 +#define G_COMPUTE_DRWS "GCOMDris::compute_drws(GCOMObservation&,"\
  51 + " GCOMSelection&, double&, double&, std::string&)"
50 52  
51 53 /* __ Macros _____________________________________________________________ */
52 54  
... ... @@ -55,6 +57,9 @@
55 57 /* __ Debug definitions __________________________________________________ */
56 58 #define G_DEBUG_COMPUTE_DRWS
57 59  
  60 +/* __ Constants __________________________________________________________ */
  61 +const double superpacket_duration = 16.384; // Duration of superpacket (s)
  62 +
58 63  
59 64 /*==========================================================================
60 65 = =
... ... @@ -338,11 +343,13 @@ void GCOMDris::extend(const GCOMDris&amp; dris)
338 343  
339 344  
340 345 /***********************************************************************//**
341   - * @brief Compute backgroud weighting cubes
  346 + * @brief Compute background weighting cubes
342 347 *
343 348 * @param[in] obs COMPTEL observation.
344 349 * @param[in] select Selection set.
345 350 * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  351 + * @param[in] timebin Time binning for rate determination (sec).
  352 + * @param[in] method Method to compute DRWs ("energy" or "phibar").
346 353 *
347 354 * @exception GException::invalid_value
348 355 * DRW with mismatching definition encountered in container.
... ... @@ -350,34 +357,18 @@ void GCOMDris::extend(const GCOMDris&amp; dris)
350 357 * No event list found in COMPTEL observation.
351 358 *
352 359 * Compute DRW cubes for a COMPTEL observation. DRW cubes are event-rate
353   - * weighted geometry cubes, where event rate are time and energy-dependent.
354   - *
355   - * Three hard-wired energy ranges are defined for which the event rate is
356   - * determined on a 5 minutes basis:
357   - *
358   - * E < 1.8 MeV
359   - * 1.8 < E < 4.3 MeV
360   - * E > 4.3 MeV
361   - *
362   - * For each superpacket, these event rates are linearly interpolated to the
363   - * relevant superpacket time.
  360 + * weighted geometry cubes which were multiplied by the solid angle of
  361 + * the (Chi,Psi) pixels.
364 362 *
365 363 * The DRW cubes are normalised so that the sum of each Phibar layer equals
366 364 * to unity.
367 365 ***************************************************************************/
368 366 void GCOMDris::compute_drws(const GCOMObservation& obs,
369 367 const GCOMSelection& select,
370   - const double& zetamin)
  368 + const double& zetamin,
  369 + const double& timebin,
  370 + const std::string& method)
371 371 {
372   - // Debug
373   - #if defined(G_DEBUG_COMPUTE_DRWS)
374   - std::cout << "GCOMDris::compute_drws" << std::endl;
375   - std::cout << "======================" << std::endl;
376   - #endif
377   -
378   - // Set constants
379   - const double dtime = 5.0 * 60.0; // Event rate integration time in seconds
380   -
381 372 // Continue only if there are DRWs in container
382 373 if (!is_empty()) {
383 374  
... ... @@ -387,9 +378,6 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
387 378 std::cout << "--------------" << std::endl;
388 379 #endif
389 380  
390   - // Initialise D1 & D2 module status
391   - GCOMStatus status;
392   -
393 381 // Get pointer to first DRW in container (for convenient access to
394 382 // members of GCOMDri; we won't use it for modifications, hence we
395 383 // can declare it as constant)
... ... @@ -433,219 +421,16 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
433 421 throw GException::invalid_argument(G_COMPUTE_DRWS, msg);
434 422 }
435 423  
436   - // Debug
437   - #if defined(G_DEBUG_COMPUTE_DRWS)
438   - std::cout << "Determine energy-dependent event rates" << std::endl;
439   - std::cout << "--------------------------------------" << std::endl;
440   - #endif
441   -
442   - // Initialise current time and vectors for computation of node function
443   - // and vectors for energy-dependent rate interpolation
444   - double time;
445   - double rate0 = 0.0;
446   - double rate1 = 0.0;
447   - double rate2 = 0.0;
448   - GNodeArray times;
449   - std::vector<double> rates0;
450   - std::vector<double> rates1;
451   - std::vector<double> rates2;
452   -
453   - // Set time of first event
454   - if (events->size() > 0) {
455   - time = (*events)[0]->time().secs();
  424 + // Method A: energy-dependent weighting
  425 + if (method == "energy") {
  426 + compute_drws_energy(obs, events, select, zetamin, timebin);
456 427 }
457 428  
458   - // Determine energy-dependent event rates. The standard event selection
459   - // criteria and earth horizon cut is applied so that the event rate
460   - // corresponds to the current selection, yet we do not consider any
461   - // Good Time Intervals or OAD time information as we later only access
462   - // relevant times (YET THERE COULD BE AN ISSUE WITH A GLITCH). Rates
463   - // are computed for three hard-coded energy intervals.
464   - for (int i_evt = 0; i_evt < events->size(); ++i_evt) {
465   -
466   - // Get pointer to event
467   - const GCOMEventAtom* event = (*events)[i_evt];
468   -
469   - // Apply event selection
470   - if (!select.use_event(*event)) {
471   - continue;
472   - }
473   -
474   - // Apply Earth horizon cut
475   - if ((event->eha() - event->phibar()) < zetamin) {
476   - continue;
477   - }
478   -
479   - // If time exceeds time bin then add rates
480   - while (event->time().secs() > (time + dtime)) {
481   -
482   - // If rates are not all zero then append rates to vectors
483   - if ((rate0 != 0.0) || (rate1 != 0.0) || (rate2 != 0.0)) {
484   -
485   - // Append time at mid-point of interval
486   - times.append(time + 0.5 * dtime);
487   -
488   - // Append rates
489   - rates0.push_back(rate0);
490   - rates1.push_back(rate1);
491   - rates2.push_back(rate2);
492   -
493   - // Debug
494   - #if defined(G_DEBUG_COMPUTE_DRWS)
495   - std::cout << "t=" << time + 0.5 * dtime;
496   - std::cout << " r0=" << rate0;
497   - std::cout << " r1=" << rate1;
498   - std::cout << " r2=" << rate2;
499   - std::cout << std::endl;
500   - #endif
501   -
502   - // Initialise rates for new accumulation
503   - rate0 = 0.0;
504   - rate1 = 0.0;
505   - rate2 = 0.0;
506   -
507   - } // endif: there were rates to append
508   -
509   - // Increment time
510   - time += dtime;
511   -
512   - } // endwhile: event time was after the current integration time stop
513   -
514   - // Accumulate rates
515   - double energy = event->energy().MeV();
516   - if (energy < 1.8) {
517   - rate0 += 1.0;
518   - }
519   - else if (energy < 4.3) {
520   - rate1 += 1.0;
521   - }
522   - else {
523   - rate2 += 1.0;
524   - }
525   -
526   - } // endfor: looped over events
527   -
528   - // Debug
529   - #if defined(G_DEBUG_COMPUTE_DRWS)
530   - std::cout << "Compute DRWs" << std::endl;
531   - std::cout << "------------" << std::endl;
532   - #endif
533   -
534   - // Get Good Time Intervals. If a pulsar selection is specified then
535   - // reduce the Good Time Intervals to the validity intervals of the
536   - // pulsar emphemerides.
537   - GCOMTim tim = obs.tim();
538   - if (select.has_pulsar()) {
539   - tim.reduce(select.pulsar().validity());
  429 + // Method B: Phibar-dependent weighting
  430 + else if (method == "phibar") {
  431 + compute_drws_phibar(obs, events, select, zetamin, timebin);
540 432 }
541 433  
542   - // Loop over Orbit Aspect Data
543   - for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
544   -
545   - // Get reference to Orbit Aspect Data of superpacket
546   - const GCOMOad &oad = obs.oads()[i_oad];
547   -
548   - // Check superpacket usage for all DRWs so that their statistics
549   - // get correctly updated
550   - bool skip = false;
551   - for (int i = 0; i < size(); ++i) {
552   - if (!m_dris[i].use_superpacket(oad, tim, select)) {
553   - skip = true;
554   - }
555   - }
556   - if (skip) {
557   - continue;
558   - }
559   -
560   - // Get Orbit Aspect Data mid-time in seconds
561   - double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
562   -
563   - // Compute energy-dependent rates
564   - std::vector<double> rates;
565   - for (int i = 0; i < size(); ++i) {
566   - if (m_dris[i].m_ebounds.emin().MeV() > 4.3) {
567   - rates.push_back(times.interpolate(time, rates2));
568   - }
569   - else if (m_dris[i].m_ebounds.emin().MeV() > 1.8) {
570   - rates.push_back(times.interpolate(time, rates1));
571   - }
572   - else {
573   - rates.push_back(times.interpolate(time, rates0));
574   - }
575   - }
576   -
577   - // Debug
578   - #if defined(G_DEBUG_COMPUTE_DRWS)
579   - std::cout << "time=" << time;
580   - for (int i = 0; i < size(); ++i) {
581   - std::cout << " r[" << i << "]=" << rates[i];
582   - }
583   - std::cout << std::endl;
584   - #endif
585   -
586   - // Prepare Earth horizon angle computation. The celestial system
587   - // is reinterpreted as the COMPTEL coordinate system, where the
588   - // 90 degrees - zenith angle becomes the declination and the azimuth
589   - // angle becomes Right Ascension. This allows us later to use the
590   - // GSkyDir::dist method to compute the distance between the geocentre
591   - // and the telescope Z-axis.
592   - double theta_geocentre = double(oad.gcel());
593   - double phi_geocentre = double(oad.gcaz());
594   - GSkyDir geocentre_comptel;
595   - geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
596   -
597   - // Compute geometrical factors for all (Chi, Psi). We assume here
598   - // that their size and definition is identical, which was checked
599   - // before
600   - for (int index = 0; index < dri->m_dri.npix(); ++index) {
601   -
602   - // Get sky direction and solid angle for (Chi, Psi)
603   - GSkyDir sky = dri->m_dri.inx2dir(index);
604   - double omega = dri->m_dri.solidangle(index);
605   -
606   - // Convert sky direction to COMPTEL coordinates
607   - double theta = oad.theta(sky);
608   - double phi = oad.phi(sky);
609   -
610   - // Compute geometric factor summed over all D1, D2 times
611   - // the solid angle of the weighting cube pixel
612   - double geometry = dri->compute_geometry(oad.tjd(), theta, phi, select, status) *
613   - omega;
614   -
615   - // Compute Earth horizon angle as the distance between the Earth
616   - // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
617   - // COMPTEL coordinates minus the radius of the Earth.
618   - GSkyDir chipsi_comptel;
619   - chipsi_comptel.radec_deg(phi, 90.0-theta);
620   - double eha = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
621   -
622   - // Loop over Phibar
623   - for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
624   -
625   - // Compute minimum Earth horizon angle for Phibar layer
626   - double ehamin = double(iphibar) * dri->m_phibin + zetamin;
627   -
628   - // Add up weighted geometry if the Earth horizon angle is
629   - // equal or larger than the minimum
630   - if (eha >= ehamin) {
631   -
632   - // Get data space index
633   - int inx = index + iphibar * dri->m_dri.npix();
634   -
635   - // Loop over DRWs and multiply geometry with appopriate
636   - // rates
637   - for (int i = 0; i < size(); ++i) {
638   - m_dris[i][inx] += geometry * rates[i];
639   - }
640   -
641   - } // endif: Earth horizon cut passed
642   -
643   - } // endfor: looped over Phibar
644   -
645   - } // endfor: looped over (Chi, Psi)
646   -
647   - } // endfor: looped over Orbit Aspect Data
648   -
649 434 // Normalise Phibar layers of all DRWs to unity
650 435 for (int i = 0; i < size(); ++i) {
651 436  
... ... @@ -676,7 +461,7 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
676 461  
677 462 } // endfor: looped over DRWs
678 463  
679   - } // endif: container was not empty
  464 + } // endif: there were DRWs in container
680 465  
681 466 // Return
682 467 return;
... ... @@ -753,3 +538,600 @@ void GCOMDris::free_members(void)
753 538 // Return
754 539 return;
755 540 }
  541 +
  542 +
  543 +/***********************************************************************//**
  544 + * @brief Compute background weighting cubes using energy dependent rates
  545 + *
  546 + * @param[in] obs COMPTEL observation.
  547 + * @param[in] event COMPTEL event list.
  548 + * @param[in] select Selection set.
  549 + * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  550 + * @param[in] timebin Time binning for rate determination (sec).
  551 + *
  552 + * Compute DRW cubes for a COMPTEL observation. DRW cubes are event-rate
  553 + * weighted geometry cubes which were multiplied by the solid angle of
  554 + * the (Chi,Psi) pixels.
  555 + *
  556 + * For this method the event rates depend on time and energy. Three
  557 + * hard-wired energy ranges are defined for which the event rate is
  558 + * determined:
  559 + *
  560 + * E < 1.8 MeV
  561 + * 1.8 < E < 4.3 MeV
  562 + * E > 4.3 MeV
  563 + *
  564 + * For each superpacket, the event rates are linearly interpolated to the
  565 + * relevant superpacket time.
  566 + *
  567 + * The DRW cubes are normalised so that the sum of each Phibar layer equals
  568 + * to unity.
  569 + ***************************************************************************/
  570 +void GCOMDris::compute_drws_energy(const GCOMObservation& obs,
  571 + const GCOMEventList* events,
  572 + const GCOMSelection& select,
  573 + const double& zetamin,
  574 + const double& timebin)
  575 +{
  576 + // Debug
  577 + #if defined(G_DEBUG_COMPUTE_DRWS)
  578 + std::cout << "GCOMDris::compute_drws_energy" << std::endl;
  579 + std::cout << "=============================" << std::endl;
  580 + #endif
  581 +
  582 + // Initialise D1 & D2 module status
  583 + GCOMStatus status;
  584 +
  585 + // Get pointer to first DRW in container (for convenient access to members
  586 + // of GCOMDri; we won't use it for modifications, hence we can declare it
  587 + // as constant)
  588 + const GCOMDri* dri = &(m_dris[0]);
  589 +
  590 + // Debug
  591 + #if defined(G_DEBUG_COMPUTE_DRWS)
  592 + std::cout << "Determine energy-dependent event rates" << std::endl;
  593 + std::cout << "--------------------------------------" << std::endl;
  594 + #endif
  595 +
  596 + // Initialise current time and vectors for computation of node function
  597 + // and vectors for energy-dependent rate interpolation
  598 + double time;
  599 + double rate0 = 0.0;
  600 + double rate1 = 0.0;
  601 + double rate2 = 0.0;
  602 + GNodeArray times;
  603 + std::vector<double> rates0;
  604 + std::vector<double> rates1;
  605 + std::vector<double> rates2;
  606 + bool set = false;
  607 +
  608 + // Set time of first event
  609 + if (events->size() > 0) {
  610 + time = (*events)[0]->time().secs();
  611 + }
  612 +
  613 + // Determine energy-dependent event rates. The standard event selection
  614 + // criteria and earth horizon cut is applied so that the event rate
  615 + // corresponds to the current selection, yet we do not consider any
  616 + // Good Time Intervals or OAD time information as we later only access
  617 + // relevant times (YET THERE COULD BE AN ISSUE WITH A GLITCH). Rates
  618 + // are computed for three hard-coded energy intervals.
  619 + for (int i_evt = 0; i_evt < events->size(); ++i_evt) {
  620 +
  621 + // Get pointer to event
  622 + const GCOMEventAtom* event = (*events)[i_evt];
  623 +
  624 + // Apply event selection
  625 + if (!select.use_event(*event)) {
  626 + continue;
  627 + }
  628 +
  629 + // Apply Earth horizon cut. Note that we need to correct later
  630 + // for this cut to remove the additional time variation coming
  631 + // from the cut.
  632 + if ((event->eha() - event->phibar()) < zetamin) {
  633 + continue;
  634 + }
  635 +
  636 + // If time exceeds time bin then add rates
  637 + while (event->time().secs() > (time + timebin)) {
  638 +
  639 + // If rates are not all zero then append rates to vectors
  640 + if (set) {
  641 +
  642 + // Append time at mid-point of interval
  643 + times.append(time + 0.5 * timebin);
  644 +
  645 + // Append rates
  646 + rates0.push_back(rate0);
  647 + rates1.push_back(rate1);
  648 + rates2.push_back(rate2);
  649 +
  650 + // Debug
  651 + #if defined(G_DEBUG_COMPUTE_DRWS)
  652 + std::cout << "t=" << time + 0.5 * timebin;
  653 + std::cout << " r0=" << rate0;
  654 + std::cout << " r1=" << rate1;
  655 + std::cout << " r2=" << rate2;
  656 + std::cout << std::endl;
  657 + #endif
  658 +
  659 + // Initialise rates for new accumulation
  660 + rate0 = 0.0;
  661 + rate1 = 0.0;
  662 + rate2 = 0.0;
  663 + set = false;
  664 +
  665 + } // endif: there were rates to append
  666 +
  667 + // Increment time
  668 + time += timebin;
  669 +
  670 + } // endwhile: event time was after the current integration time stop
  671 +
  672 + // Accumulate rates
  673 + double energy = event->energy().MeV();
  674 + if (energy < 1.8) {
  675 + rate0 += 1.0;
  676 + set = true;
  677 + }
  678 + else if (energy < 4.3) {
  679 + rate1 += 1.0;
  680 + set = true;
  681 + }
  682 + else {
  683 + rate2 += 1.0;
  684 + set = true;
  685 + }
  686 +
  687 + } // endfor: looped over events
  688 +
  689 + // Append remaining rates
  690 + if (set) {
  691 +
  692 + // Append time at mid-point of interval
  693 + times.append(time + 0.5 * timebin);
  694 +
  695 + // Append rates
  696 + rates0.push_back(rate0);
  697 + rates1.push_back(rate1);
  698 + rates2.push_back(rate2);
  699 +
  700 + } // endif: appended remaining rates
  701 +
  702 + // Debug
  703 + #if defined(G_DEBUG_COMPUTE_DRWS)
  704 + std::cout << "Compute DRWs" << std::endl;
  705 + std::cout << "------------" << std::endl;
  706 + #endif
  707 +
  708 + // Allocate geometry factor and Earth horizon angle working arrays
  709 + std::vector<double> geometry(dri->m_dri.npix());
  710 + std::vector<double> eha(dri->m_dri.npix());
  711 +
  712 + // Get Good Time Intervals. If a pulsar selection is specified then
  713 + // reduce the Good Time Intervals to the validity intervals of the
  714 + // pulsar emphemerides.
  715 + GCOMTim tim = obs.tim();
  716 + if (select.has_pulsar()) {
  717 + tim.reduce(select.pulsar().validity());
  718 + }
  719 +
  720 + // Loop over Orbit Aspect Data
  721 + for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
  722 +
  723 + // Get reference to Orbit Aspect Data of superpacket
  724 + const GCOMOad &oad = obs.oads()[i_oad];
  725 +
  726 + // Check superpacket usage for all DRWs so that their statistics
  727 + // get correctly updated
  728 + bool skip = false;
  729 + for (int i = 0; i < size(); ++i) {
  730 + if (!m_dris[i].use_superpacket(oad, tim, select)) {
  731 + skip = true;
  732 + }
  733 + }
  734 + if (skip) {
  735 + continue;
  736 + }
  737 +
  738 + // Get Orbit Aspect Data mid-time in seconds
  739 + double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
  740 +
  741 + // Compute energy-dependent rates
  742 + std::vector<double> rates;
  743 + for (int i = 0; i < size(); ++i) {
  744 + if (m_dris[i].m_ebounds.emin().MeV() > 4.3) {
  745 + rates.push_back(times.interpolate(time, rates2));
  746 + }
  747 + else if (m_dris[i].m_ebounds.emin().MeV() > 1.8) {
  748 + rates.push_back(times.interpolate(time, rates1));
  749 + }
  750 + else {
  751 + rates.push_back(times.interpolate(time, rates0));
  752 + }
  753 + }
  754 +
  755 + // Debug
  756 + #if defined(G_DEBUG_COMPUTE_DRWS)
  757 + std::cout << "time=" << time;
  758 + for (int i = 0; i < size(); ++i) {
  759 + std::cout << " r[" << i << "]=" << rates[i];
  760 + }
  761 + std::cout << std::endl;
  762 + #endif
  763 +
  764 + // Prepare Earth horizon angle computation. The celestial system
  765 + // is reinterpreted as the COMPTEL coordinate system, where the
  766 + // 90 degrees - zenith angle becomes the declination and the azimuth
  767 + // angle becomes Right Ascension. This allows us later to use the
  768 + // GSkyDir::dist method to compute the distance between the geocentre
  769 + // and the telescope Z-axis.
  770 + double theta_geocentre = double(oad.gcel());
  771 + double phi_geocentre = double(oad.gcaz());
  772 + GSkyDir geocentre_comptel;
  773 + geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
  774 +
  775 + // Compute solid-angle-corrected geometry factors and Earth horizon
  776 + // angles for all (Chi,Psi) pixels
  777 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  778 +
  779 + // Get sky direction and solid angle for (Chi, Psi)
  780 + GSkyDir sky = dri->m_dri.inx2dir(index);
  781 + double omega = dri->m_dri.solidangle(index);
  782 +
  783 + // Convert sky direction to COMPTEL coordinates
  784 + double theta = oad.theta(sky);
  785 + double phi = oad.phi(sky);
  786 +
  787 + // Compute geometric factor summed over all D1, D2 times
  788 + // the solid angle of the weighting cube pixel
  789 + geometry[index] = dri->compute_geometry(oad.tjd(), theta, phi,
  790 + select, status) * omega;
  791 +
  792 + // Compute Earth horizon angle as the distance between the Earth
  793 + // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
  794 + // COMPTEL coordinates minus the radius of the Earth.
  795 + GSkyDir chipsi_comptel;
  796 + chipsi_comptel.radec_deg(phi, 90.0-theta);
  797 + eha[index] = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
  798 +
  799 + } // endfor: computed solid-angle-corrected geometry factors
  800 +
  801 + // Loop over all DRWs
  802 + for (int i = 0; i < size(); ++i) {
  803 +
  804 + // Compute Earth horizon cut correction factor
  805 + double rates_sum_all = 0.0;
  806 + double rates_sum_cut = 0.0;
  807 + double ehacorr = 1.0;
  808 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  809 +
  810 + // Sum product of geometry and rates for all pixels
  811 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  812 + rates_sum_all += geometry[index] * rates[i];
  813 + }
  814 +
  815 + // Sum product of geometry and rates for pixels passing
  816 + // the Earth horizon cut
  817 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  818 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  819 + if (eha[index] >= ehamin) {
  820 + rates_sum_cut += geometry[index] * rates[i];
  821 + }
  822 + }
  823 +
  824 + } // endfor: looped over Phibar
  825 +
  826 + // If something passes the Earth horizon cut then compute the
  827 + // correction factor now
  828 + if (rates_sum_cut != 0.0) {
  829 + ehacorr = rates_sum_all / rates_sum_cut;
  830 + }
  831 +
  832 + // Loop over all Phibar layers
  833 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  834 +
  835 + // Compute minimum Earth horizon angle for Phibar layer
  836 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  837 +
  838 + // Loop over all (Chi,Psi) pixels
  839 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  840 +
  841 + // If Earth horizon angle is equal or larger than the minimum
  842 + // requirement then add up the rate weighted geometry factor
  843 + if (eha[index] >= ehamin) {
  844 +
  845 + // Get data space index
  846 + int inx = index + iphibar * dri->m_dri.npix();
  847 +
  848 + // Store product as DRW value
  849 + m_dris[i][inx] += geometry[index] * rates[i] * ehacorr;
  850 +
  851 + } // endif: Earth horizon cut passed
  852 +
  853 + } // endfor: looped over (Chi, Psi)
  854 +
  855 + } // endfor: looped over Phibar
  856 +
  857 + } // endfor: looped over all DRWs
  858 +
  859 + } // endfor: looped over Orbit Aspect Data
  860 +
  861 + // Return
  862 + return;
  863 +}
  864 +
  865 +
  866 +/***********************************************************************//**
  867 + * @brief Compute background weighting cubes using Phibar dependent rates
  868 + *
  869 + * @param[in] obs COMPTEL observation.
  870 + * @param[in] event COMPTEL event list.
  871 + * @param[in] select Selection set.
  872 + * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
  873 + * @param[in] timebin Time binning for rate determination (sec).
  874 + *
  875 + * Compute DRW cubes for a COMPTEL observation. DRW cubes are event-rate
  876 + * weighted geometry cubes which were multiplied by the solid angle of
  877 + * the (Chi,Psi) pixels.
  878 + *
  879 + * For this method the event rates depend on time and Phibar angle. For
  880 + * each superpacket, the event rates are linearly interpolated to the
  881 + * relevant superpacket time.
  882 + *
  883 + * The DRW cubes are normalised so that the sum of each Phibar layer equals
  884 + * to unity.
  885 + ***************************************************************************/
  886 +void GCOMDris::compute_drws_phibar(const GCOMObservation& obs,
  887 + const GCOMEventList* events,
  888 + const GCOMSelection& select,
  889 + const double& zetamin,
  890 + const double& timebin)
  891 +{
  892 + // Debug
  893 + #if defined(G_DEBUG_COMPUTE_DRWS)
  894 + std::cout << "GCOMDris::compute_drws_phibar" << std::endl;
  895 + std::cout << "=============================" << std::endl;
  896 + #endif
  897 +
  898 + // Initialise D1 & D2 module status
  899 + GCOMStatus status;
  900 +
  901 + // Get pointer to first DRW in container (for convenient access to members
  902 + // of GCOMDri; we won't use it for modifications, hence we can declare it
  903 + // as constant)
  904 + const GCOMDri* dri = &(m_dris[0]);
  905 +
  906 + // Debug
  907 + #if defined(G_DEBUG_COMPUTE_DRWS)
  908 + std::cout << "Determine Phibar-dependent event rates" << std::endl;
  909 + std::cout << "--------------------------------------" << std::endl;
  910 + #endif
  911 +
  912 + // Initialise current time and vectors for computation of node function
  913 + // and vectors for energy-dependent rate interpolation
  914 + double time;
  915 + std::vector<double> rate(dri->nphibar());
  916 + GNodeArray times;
  917 + std::vector<std::vector<double> > rates(dri->nphibar());
  918 + bool set = false;
  919 +
  920 + // Set time of first event
  921 + if (events->size() > 0) {
  922 + time = (*events)[0]->time().secs();
  923 + }
  924 +
  925 + // Determine energy-dependent event rates. The standard event selection
  926 + // criteria and earth horizon cut is applied so that the event rate
  927 + // corresponds to the current selection, yet we do not consider any
  928 + // Good Time Intervals or OAD time information as we later only access
  929 + // relevant times.
  930 + for (int i_evt = 0; i_evt < events->size(); ++i_evt) {
  931 +
  932 + // Get pointer to event
  933 + const GCOMEventAtom* event = (*events)[i_evt];
  934 +
  935 + // Apply event selection
  936 + if (!select.use_event(*event)) {
  937 + continue;
  938 + }
  939 +
  940 + // Apply Earth horizon cut. Note that we need to correct later
  941 + // for this cut to remove the additional time variation coming
  942 + // from the cut.
  943 + if ((event->eha() - event->phibar()) < zetamin) {
  944 + continue;
  945 + }
  946 +
  947 + // If time exceeds time bin then add rates
  948 + while (event->time().secs() > (time + timebin)) {
  949 +
  950 + // If rates are set then append rates to vectors
  951 + if (set) {
  952 +
  953 + // Append time at mid-point of interval
  954 + times.append(time + 0.5 * timebin);
  955 +
  956 + // Append rates and initialise for new accumulation
  957 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  958 + rates[iphibar].push_back(rate[iphibar]);
  959 + rate[iphibar] = 0.0;
  960 + }
  961 +
  962 + // Signal that rates were not yet set
  963 + set = false;
  964 +
  965 + } // endif: there were rates to append
  966 +
  967 + // Increment time
  968 + time += timebin;
  969 +
  970 + } // endwhile: event time was after the current integration time stop
  971 +
  972 + // Accumulate rates
  973 + int iphibar = int((event->phibar() - dri->phimin()) / dri->phibin());
  974 + if (iphibar >= 0 and iphibar < dri->nphibar()) {
  975 + rate[iphibar] += 1.0;
  976 + set = true;
  977 + }
  978 +
  979 + } // endfor: looped over events
  980 +
  981 + // Append remaining rates
  982 + if (set) {
  983 +
  984 + // Append time at mid-point of interval
  985 + times.append(time + 0.5 * timebin);
  986 +
  987 + // Append rates
  988 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  989 + rates[iphibar].push_back(rate[iphibar]);
  990 + }
  991 +
  992 + } // endif: appended remaining rates
  993 +
  994 + // Debug
  995 + #if defined(G_DEBUG_COMPUTE_DRWS)
  996 + std::cout << "Compute DRWs" << std::endl;
  997 + std::cout << "------------" << std::endl;
  998 + #endif
  999 +
  1000 + // Allocate geometry factor and Earth horizon angle working arrays
  1001 + std::vector<double> geometry(dri->m_dri.npix());
  1002 + std::vector<double> eha(dri->m_dri.npix());
  1003 +
  1004 + // Get Good Time Intervals. If a pulsar selection is specified then
  1005 + // reduce the Good Time Intervals to the validity intervals of the
  1006 + // pulsar emphemerides.
  1007 + GCOMTim tim = obs.tim();
  1008 + if (select.has_pulsar()) {
  1009 + tim.reduce(select.pulsar().validity());
  1010 + }
  1011 +
  1012 + // Loop over Orbit Aspect Data
  1013 + for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
  1014 +
  1015 + // Get reference to Orbit Aspect Data of superpacket
  1016 + const GCOMOad &oad = obs.oads()[i_oad];
  1017 +
  1018 + // Check superpacket usage for all DRWs so that their statistics
  1019 + // get correctly updated
  1020 + bool skip = false;
  1021 + for (int i = 0; i < size(); ++i) {
  1022 + if (!m_dris[i].use_superpacket(oad, tim, select)) {
  1023 + skip = true;
  1024 + }
  1025 + }
  1026 + if (skip) {
  1027 + continue;
  1028 + }
  1029 +
  1030 + // Get Orbit Aspect Data mid-time in seconds
  1031 + double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
  1032 +
  1033 + // Compute Phibar-dependent rates by interpolation of precomputed
  1034 + // rate vectors
  1035 + std::vector<double> rate_oad(dri->nphibar());
  1036 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  1037 + rate_oad[iphibar] = times.interpolate(time, rates[iphibar]);
  1038 + }
  1039 +
  1040 + // Prepare Earth horizon angle computation. The celestial system
  1041 + // is reinterpreted as the COMPTEL coordinate system, where the
  1042 + // 90 degrees - zenith angle becomes the declination and the azimuth
  1043 + // angle becomes Right Ascension. This allows us later to use the
  1044 + // GSkyDir::dist method to compute the distance between the geocentre
  1045 + // and the telescope Z-axis.
  1046 + double theta_geocentre = double(oad.gcel());
  1047 + double phi_geocentre = double(oad.gcaz());
  1048 + GSkyDir geocentre_comptel;
  1049 + geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
  1050 +
  1051 + // Compute solid-angle-corrected geometry factors and Earth horizon
  1052 + // angles for all (Chi,Psi) pixels
  1053 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  1054 +
  1055 + // Get sky direction and solid angle for (Chi, Psi)
  1056 + GSkyDir sky = dri->m_dri.inx2dir(index);
  1057 + double omega = dri->m_dri.solidangle(index);
  1058 +
  1059 + // Convert sky direction to COMPTEL coordinates
  1060 + double theta = oad.theta(sky);
  1061 + double phi = oad.phi(sky);
  1062 +
  1063 + // Compute geometric factor summed over all D1, D2 times
  1064 + // the solid angle of the weighting cube pixel
  1065 + geometry[index] = dri->compute_geometry(oad.tjd(), theta, phi,
  1066 + select, status) * omega;
  1067 +
  1068 + // Compute Earth horizon angle as the distance between the Earth
  1069 + // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
  1070 + // COMPTEL coordinates minus the radius of the Earth.
  1071 + GSkyDir chipsi_comptel;
  1072 + chipsi_comptel.radec_deg(phi, 90.0-theta);
  1073 + eha[index] = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
  1074 +
  1075 + } // endfor: computed solid-angle-corrected geometry factors
  1076 +
  1077 + // Correct rates for Earth horizon cut
  1078 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  1079 +
  1080 + // Compute minimum Earth Horizon angle for Phibar layer
  1081 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  1082 +
  1083 + // Initialise sums
  1084 + double rates_sum_all = 0.0;
  1085 + double rates_sum_cut = 0.0;
  1086 +
  1087 + // Sum geometry for all pixels and pixels passing the Earth horizon cut
  1088 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  1089 + rates_sum_all += geometry[index];
  1090 + if (eha[index] >= ehamin) {
  1091 + rates_sum_cut += geometry[index];
  1092 + }
  1093 + }
  1094 +
  1095 + // If something passes the Earth horizon cut then correct
  1096 + // the rate
  1097 + if (rates_sum_cut != 0.0) {
  1098 + rate_oad[iphibar] *= rates_sum_all / rates_sum_cut;
  1099 + }
  1100 +
  1101 + } // endfor: looped over Phibar
  1102 +
  1103 + // Loop over all DRWs
  1104 + for (int i = 0; i < size(); ++i) {
  1105 +
  1106 + // Loop over all Phibar layers
  1107 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  1108 +
  1109 + // Compute minimum Earth horizon angle for Phibar layer
  1110 + double ehamin = double(iphibar) * dri->m_phibin + zetamin;
  1111 +
  1112 + // Loop over all (Chi,Psi) pixels
  1113 + for (int index = 0; index < dri->m_dri.npix(); ++index) {
  1114 +
  1115 + // If Earth horizon angle is equal or larger than the minimum
  1116 + // requirement then add up the rate weighted geometry factor
  1117 + if (eha[index] >= ehamin) {
  1118 +
  1119 + // Get data space index
  1120 + int inx = index + iphibar * dri->m_dri.npix();
  1121 +
  1122 + // Store product as DRW value
  1123 + m_dris[i][inx] += geometry[index] * rate_oad[iphibar];
  1124 +
  1125 + } // endif: Earth horizon cut passed
  1126 +
  1127 + } // endfor: looped over (Chi, Psi)
  1128 +
  1129 + } // endfor: looped over Phibar
  1130 +
  1131 + } // endfor: looped over all DRWs
  1132 +
  1133 + } // endfor: looped over Orbit Aspect Data
  1134 +
  1135 + // Return
  1136 + return;
  1137 +}
... ...