gammalib / gammalib

Commit b0b38ecbab27c7b93dcc2dd581489afa6358540b

Authored by Jürgen Knödlseder
1 parent 6f9976c0

Implement handling of DRWs

Inline Side-by-side
Showing 4 changed files with 209 additions and 42 deletions
inst/com/include/GCOMObservation.hpp
  View file @b0b38ec
1 1 /***************************************************************************
2 2 * GCOMObservation.hpp - COMPTEL observation class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2022 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,10 +70,12 @@ public:
70 70 explicit GCOMObservation(const GXmlElement& xml);
71 71 GCOMObservation(const GCOMDri& dre,
72 72 const GCOMDri& drb,
  73 + const GCOMDri& drw,
73 74 const GCOMDri& drg,
74 75 const GCOMDri& drx);
75 76 GCOMObservation(const GFilename& drename,
76 77 const GFilename& drbname,
  78 + const GFilename& drwname,
77 79 const GFilename& drgname,
78 80 const GFilename& drxname);
79 81 GCOMObservation(const GFilename& evpname,
... ... @@ -108,6 +110,7 @@ public:
108 110 bool is_binned(void) const;
109 111 void load(const GFilename& drename,
110 112 const GFilename& drbname,
  113 + const GFilename& drwname,
111 114 const GFilename& drgname,
112 115 const GFilename& drxname);
113 116 void load(const GFilename& evpname,
... ... @@ -124,6 +127,7 @@ public:
124 127 const double& obs_id(void) const;
125 128 const double& ewidth(void) const;
126 129 const GCOMDri& drb(void) const;
  130 + const GCOMDri& drw(void) const;
127 131 const GCOMDri& drg(void) const;
128 132 const GCOMDri& drx(void) const;
129 133 GCOMDri drm(const GModels& models) const;
... ... @@ -135,6 +139,7 @@ public:
135 139 void bvcs(const GCOMBvcs& bvcs);
136 140 const GFilename& drename(void) const;
137 141 const GFilename& drbname(void) const;
  142 + const GFilename& drwname(void) const;
138 143 const GFilename& drgname(void) const;
139 144 const GFilename& drxname(void) const;
140 145 const GFilename& rspname(void) const;
... ... @@ -142,6 +147,7 @@ public:
142 147 const int& phi_last(void) const;
143 148 void drename(const GFilename& drename);
144 149 void drbname(const GFilename& drbname);
  150 + void drwname(const GFilename& drwname);
145 151 void drgname(const GFilename& drgname);
146 152 void drxname(const GFilename& drxname);
147 153 void rspname(const GFilename& rspname);
... ... @@ -161,6 +167,7 @@ protected:
161 167 void free_members(void);
162 168 void load_dre(const GFilename& drename);
163 169 void load_drb(const GFilename& drbname);
  170 + void load_drw(const GFilename& drwname);
164 171 void load_drg(const GFilename& drgname);
165 172 void load_drx(const GFilename& drxname);
166 173 bool check_dri(const GCOMDri& map) const;
... ... @@ -200,10 +207,12 @@ protected:
200 207 // Protected members for binned observation
201 208 GFilename m_drename; //!< DRE filename
202 209 GFilename m_drbname; //!< DRB filename
  210 + GFilename m_drwname; //!< DRW filename
203 211 GFilename m_drgname; //!< DRG filename
204 212 GFilename m_drxname; //!< DRX filename
205 213 GFilename m_rspname; //!< Response cache filename
206 214 GCOMDri m_drb; //!< Background model
  215 + GCOMDri m_drw; //!< Weighting cube
207 216 GCOMDri m_drg; //!< Geometry factors
208 217 GCOMDri m_drx; //!< Exposure map
209 218 double m_ewidth; //!< Energy width (MeV)
... ... @@ -418,6 +427,19 @@ const GCOMDri&amp; GCOMObservation::drb(void) const
418 427  
419 428  
420 429 /***********************************************************************//**
  430 + * @brief Return weighting cube
  431 + *
  432 + * @return Weighting cube.
  433 + ***************************************************************************/
  434 +inline
  435 +const GCOMDri& GCOMObservation::drw(void) const
  436 +{
  437 + // Return weighting cube
  438 + return (m_drw);
  439 +}
  440 +
  441 +
  442 +/***********************************************************************//**
421 443 * @brief Return geometry factors
422 444 *
423 445 * @return Geometry factors.
... ... @@ -572,6 +594,19 @@ const GFilename&amp; GCOMObservation::drbname(void) const
572 594  
573 595  
574 596 /***********************************************************************//**
  597 + * @brief Return DRW filename
  598 + *
  599 + * @return DRW filename.
  600 + ***************************************************************************/
  601 +inline
  602 +const GFilename& GCOMObservation::drwname(void) const
  603 +{
  604 + // Return DRW filename
  605 + return (m_drwname);
  606 +}
  607 +
  608 +
  609 +/***********************************************************************//**
575 610 * @brief Return DRG filename
576 611 *
577 612 * @return DRG filename.
... ... @@ -663,6 +698,19 @@ void GCOMObservation::drbname(const GFilename&amp; drbname)
663 698  
664 699  
665 700 /***********************************************************************//**
  701 + * @brief Set DRW filename
  702 + *
  703 + * @param[in] drwname DRW filename.
  704 + ***************************************************************************/
  705 +inline
  706 +void GCOMObservation::drwname(const GFilename& drwname)
  707 +{
  708 + m_drwname = drwname;
  709 + return;
  710 +}
  711 +
  712 +
  713 +/***********************************************************************//**
666 714 * @brief Set DRG filename
667 715 *
668 716 * @param[in] drgname DRG filename.
... ...
inst/com/pyext/GCOMObservation.i
  View file @b0b38ec
1 1 /***************************************************************************
2 2 * GCOMObservation.i - COMPTEL observation class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2022 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -41,10 +41,12 @@ public:
41 41 explicit GCOMObservation(const GXmlElement& xml);
42 42 GCOMObservation(const GCOMDri& dre,
43 43 const GCOMDri& drb,
  44 + const GCOMDri& drw,
44 45 const GCOMDri& drg,
45 46 const GCOMDri& drx);
46 47 GCOMObservation(const GFilename& drename,
47 48 const GFilename& drbname,
  49 + const GFilename& drwname,
48 50 const GFilename& drgname,
49 51 const GFilename& drxname);
50 52 GCOMObservation(const GFilename& evpname,
... ... @@ -75,6 +77,7 @@ public:
75 77 bool is_binned(void) const;
76 78 void load(const GFilename& drename,
77 79 const GFilename& drbname,
  80 + const GFilename& drwname,
78 81 const GFilename& drgname,
79 82 const GFilename& drxname);
80 83 void load(const GFilename& evpname,
... ... @@ -91,6 +94,7 @@ public:
91 94 const double& obs_id(void) const;
92 95 const double& ewidth(void) const;
93 96 const GCOMDri& drb(void) const;
  97 + const GCOMDri& drw(void) const;
94 98 const GCOMDri& drg(void) const;
95 99 const GCOMDri& drx(void) const;
96 100 GCOMDri drm(const GModels& models) const;
... ... @@ -102,6 +106,7 @@ public:
102 106 void bvcs(const GCOMBvcs& bvcs);
103 107 const GFilename& drename(void) const;
104 108 const GFilename& drbname(void) const;
  109 + const GFilename& drwname(void) const;
105 110 const GFilename& drgname(void) const;
106 111 const GFilename& drxname(void) const;
107 112 const GFilename& rspname(void) const;
... ... @@ -109,6 +114,7 @@ public:
109 114 const int& phi_last(void) const;
110 115 void drename(const GFilename& drename);
111 116 void drbname(const GFilename& drbname);
  117 + void drwname(const GFilename& drwname);
112 118 void drgname(const GFilename& drgname);
113 119 void drxname(const GFilename& drxname);
114 120 void rspname(const GFilename& rspname);
... ...
inst/com/src/GCOMDris.cpp
  View file @b0b38ec
... ... @@ -349,7 +349,21 @@ void GCOMDris::extend(const GCOMDris&amp; dris)
349 349 * @exception GException::invalid_argument
350 350 * No event list found in COMPTEL observation.
351 351 *
352   - * Compute DRW cubes for a COMPTEL observation.
  352 + * 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.
  364 + *
  365 + * The DRW cubes are normalised so that the sum of each Phibar layer equals
  366 + * to unity.
353 367 ***************************************************************************/
354 368 void GCOMDris::compute_drws(const GCOMObservation& obs,
355 369 const GCOMSelection& select,
... ... @@ -367,6 +381,12 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
367 381 // Continue only if there are DRWs in container
368 382 if (!is_empty()) {
369 383  
  384 + // Debug
  385 + #if defined(G_DEBUG_COMPUTE_DRWS)
  386 + std::cout << "Initialisation" << std::endl;
  387 + std::cout << "--------------" << std::endl;
  388 + #endif
  389 +
370 390 // Initialise D1 & D2 module status
371 391 GCOMStatus status;
372 392  
... ... @@ -397,6 +417,11 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
397 417 // Set all DRG bins to zero
398 418 m_dris[i].init_cube();
399 419  
  420 + // Debug
  421 + #if defined(G_DEBUG_COMPUTE_DRWS)
  422 + std::cout << m_dris[i].print() << std::endl;
  423 + #endif
  424 +
400 425 } // endfor: looped over all DRWs
401 426  
402 427 // Get pointer to event list. Throw an exception if the observation
... ... @@ -408,6 +433,12 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
408 433 throw GException::invalid_argument(G_COMPUTE_DRWS, msg);
409 434 }
410 435  
  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 +
411 442 // Initialise current time and vectors for computation of node function
412 443 // and vectors for energy-dependent rate interpolation
413 444 double time;
... ... @@ -494,6 +525,12 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
494 525  
495 526 } // endfor: looped over events
496 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 +
497 534 // Get Good Time Intervals. If a pulsar selection is specified then
498 535 // reduce the Good Time Intervals to the validity intervals of the
499 536 // pulsar emphemerides.
... ... @@ -523,6 +560,29 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
523 560 // Get Orbit Aspect Data mid-time in seconds
524 561 double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
525 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 +
526 586 // Prepare Earth horizon angle computation. The celestial system
527 587 // is reinterpreted as the COMPTEL coordinate system, where the
528 588 // 90 degrees - zenith angle becomes the declination and the azimuth
... ... @@ -535,19 +595,22 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
535 595 geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
536 596  
537 597 // Compute geometrical factors for all (Chi, Psi). We assume here
538   - // that their size and definition is identical, which checked on
539   - // entry
  598 + // that their size and definition is identical, which was checked
  599 + // before
540 600 for (int index = 0; index < dri->m_dri.npix(); ++index) {
541 601  
542   - // Get sky direction for (Chi, Psi)
543   - GSkyDir sky = dri->m_dri.inx2dir(index);
  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);
544 605  
545 606 // Convert sky direction to COMPTEL coordinates
546 607 double theta = oad.theta(sky);
547 608 double phi = oad.phi(sky);
548 609  
549   - // Compute geometric factor summed over all D1, D2
550   - double geometry = dri->compute_geometry(oad.tjd(), theta, phi, select, status);
  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;
551 614  
552 615 // Compute Earth horizon angle as the distance between the Earth
553 616 // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
... ... @@ -562,32 +625,18 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
562 625 // Compute minimum Earth horizon angle for Phibar layer
563 626 double ehamin = double(iphibar) * dri->m_phibin + zetamin;
564 627  
565   - // Add up geometry if the Earth horizon angle is equal or
566   - // larger than the minimum
  628 + // Add up weighted geometry if the Earth horizon angle is
  629 + // equal or larger than the minimum
567 630 if (eha >= ehamin) {
568 631  
569 632 // Get data space index
570 633 int inx = index + iphibar * dri->m_dri.npix();
571 634  
572   - // Loop over DRWs and multiply with appopriate rates
  635 + // Loop over DRWs and multiply geometry with appopriate
  636 + // rates
573 637 for (int i = 0; i < size(); ++i) {
574   -
575   - // Set rate by interpolating the precomputed vectors
576   - double rate = 0.0;
577   - if (m_dris[i].m_ebounds.emin().MeV() > 4.3) {
578   - rate = times.interpolate(time, rates2);
579   - }
580   - else if (m_dris[i].m_ebounds.emin().MeV() > 1.8) {
581   - rate = times.interpolate(time, rates1);
582   - }
583   - else {
584   - rate = times.interpolate(time, rates0);
585   - }
586   -
587   - // Set DRW
588   - m_dris[i][inx] += geometry * rate;
589   -
590   - } // endfor: looped over DRWs
  638 + m_dris[i][inx] += geometry * rates[i];
  639 + }
591 640  
592 641 } // endif: Earth horizon cut passed
593 642  
... ... @@ -597,7 +646,7 @@ void GCOMDris::compute_drws(const GCOMObservation&amp; obs,
597 646  
598 647 } // endfor: looped over Orbit Aspect Data
599 648  
600   - // Normalisation Phibar layers of all DRWs to unity
  649 + // Normalise Phibar layers of all DRWs to unity
601 650 for (int i = 0; i < size(); ++i) {
602 651  
603 652 // Loop over Phibar layers
... ...
inst/com/src/GCOMObservation.cpp
  View file @b0b38ec
1 1 /***************************************************************************
2 2 * GCOMObservation.cpp - COMPTEL Observation class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2022 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,6 +55,7 @@ const GObservationRegistry g_obs_com_registry(&amp;g_obs_com_seed);
55 55 #define G_READ "GCOMObservation::read(GXmlElement&)"
56 56 #define G_WRITE "GCOMObservation::write(GXmlElement&)"
57 57 #define G_LOAD_DRB "GCOMObservation::load_drb(GFilename&)"
  58 +#define G_LOAD_DRW "GCOMObservation::load_drw(GFilename&)"
58 59 #define G_LOAD_DRG "GCOMObservation::load_drg(GFilename&)"
59 60 #define G_DRM "GCOMObservation::drm(GModels&)"
60 61 #define G_COMPUTE_DRB "GCOMObservation::compute_drb(std::string&, GCOMDri&, "\
... ... @@ -119,6 +120,7 @@ GCOMObservation::GCOMObservation(const GXmlElement&amp; xml) : GObservation()
119 120 *
120 121 * @param[in] dre Event cube.
121 122 * @param[in] drb Background cube.
  123 + * @param[in] drw Weighting cube.
122 124 * @param[in] drg Geometry cube.
123 125 * @param[in] drx Exposure map.
124 126 *
... ... @@ -128,6 +130,7 @@ GCOMObservation::GCOMObservation(const GXmlElement&amp; xml) : GObservation()
128 130 ***************************************************************************/
129 131 GCOMObservation::GCOMObservation(const GCOMDri& dre,
130 132 const GCOMDri& drb,
  133 + const GCOMDri& drw,
131 134 const GCOMDri& drg,
132 135 const GCOMDri& drx) : GObservation()
133 136 {
... ... @@ -137,6 +140,7 @@ GCOMObservation::GCOMObservation(const GCOMDri&amp; dre,
137 140 // Store DRIs
138 141 m_events = new GCOMEventCube(dre);
139 142 m_drb = drb;
  143 + m_drw = drw;
140 144 m_drg = drg;
141 145 m_drx = drx;
142 146  
... ... @@ -149,6 +153,7 @@ GCOMObservation::GCOMObservation(const GCOMDri&amp; dre,
149 153 m_name = "unknown";
150 154 m_drename = "";
151 155 m_drbname = "";
  156 + m_drwname = "";
152 157 m_drgname = "";
153 158 m_drxname = "";
154 159  
... ... @@ -162,6 +167,7 @@ GCOMObservation::GCOMObservation(const GCOMDri&amp; dre,
162 167 *
163 168 * @param[in] drename Event cube name.
164 169 * @param[in] drbname Background cube name.
  170 + * @param[in] drwname Weighting cube name.
165 171 * @param[in] drgname Geometry cube name.
166 172 * @param[in] drxname Exposure map name.
167 173 *
... ... @@ -169,6 +175,7 @@ GCOMObservation::GCOMObservation(const GCOMDri&amp; dre,
169 175 *
170 176 * DRE - Events cube
171 177 * DRB - Background model cube
  178 + * DRW - Weighting cube
172 179 * DRG - Geometry factors cube
173 180 * DRX - Exposure map
174 181 *
... ... @@ -176,6 +183,7 @@ GCOMObservation::GCOMObservation(const GCOMDri&amp; dre,
176 183 ***************************************************************************/
177 184 GCOMObservation::GCOMObservation(const GFilename& drename,
178 185 const GFilename& drbname,
  186 + const GFilename& drwname,
179 187 const GFilename& drgname,
180 188 const GFilename& drxname) : GObservation()
181 189 {
... ... @@ -183,7 +191,7 @@ GCOMObservation::GCOMObservation(const GFilename&amp; drename,
183 191 init_members();
184 192  
185 193 // Load observation
186   - load(drename, drbname, drgname, drxname);
  194 + load(drename, drbname, drwname, drgname, drxname);
187 195  
188 196 // Return
189 197 return;
... ... @@ -459,6 +467,7 @@ double GCOMObservation::npred(const GModel&amp; model) const
459 467 * <observation name="Crab" id="000001" instrument="COM">
460 468 * <parameter name="DRE" file="m50438_dre.fits"/>
461 469 * <parameter name="DRB" file="m34997_drg.fits"/>
  470 + * <parameter name="DRW" file="m34997_drw.fits"/>
462 471 * <parameter name="DRG" file="m34997_drg.fits"/>
463 472 * <parameter name="DRX" file="m32171_drx.fits"/>
464 473 * <parameter name="IAQ" value="UNH(1.0-3.0)MeV"/>
... ... @@ -513,6 +522,7 @@ void GCOMObservation::read(const GXmlElement&amp; xml)
513 522 // Get parameters
514 523 std::string drename = gammalib::xml_get_attr(G_READ, xml, "DRE", "file");
515 524 std::string drbname = gammalib::xml_get_attr(G_READ, xml, "DRB", "file");
  525 + std::string drwname = gammalib::xml_get_attr(G_READ, xml, "DRW", "file");
516 526 std::string drgname = gammalib::xml_get_attr(G_READ, xml, "DRG", "file");
517 527 std::string drxname = gammalib::xml_get_attr(G_READ, xml, "DRX", "file");
518 528 std::string iaqname = gammalib::xml_get_attr(G_READ, xml, "IAQ", "value");
... ... @@ -520,12 +530,13 @@ void GCOMObservation::read(const GXmlElement&amp; xml)
520 530 // Expand file names
521 531 drename = gammalib::xml_file_expand(xml, drename);
522 532 drbname = gammalib::xml_file_expand(xml, drbname);
  533 + drwname = gammalib::xml_file_expand(xml, drwname);
523 534 drgname = gammalib::xml_file_expand(xml, drgname);
524 535 drxname = gammalib::xml_file_expand(xml, drxname);
525 536 std::string iaqname_expanded = gammalib::xml_file_expand(xml, iaqname);
526 537  
527 538 // Load observation
528   - load(drename, drbname, drgname, drxname);
  539 + load(drename, drbname, drwname, drgname, drxname);
529 540  
530 541 // Load IAQ by trying first the expanded name as a FITS file and
531 542 // otherwise the unexpanded name
... ... @@ -597,6 +608,7 @@ void GCOMObservation::read(const GXmlElement&amp; xml)
597 608 * <observation name="Crab" id="000001" instrument="COM">
598 609 * <parameter name="DRE" file="m50438_dre.fits"/>
599 610 * <parameter name="DRB" file="m34997_drg.fits"/>
  611 + * <parameter name="DRW" file="m34997_drw.fits"/>
600 612 * <parameter name="DRG" file="m34997_drg.fits"/>
601 613 * <parameter name="DRX" file="m32171_drx.fits"/>
602 614 * <parameter name="IAQ" value="UNH(1.0-3.0)MeV"/>
... ... @@ -651,6 +663,10 @@ void GCOMObservation::write(GXmlElement&amp; xml) const
651 663 par = gammalib::xml_need_par(G_WRITE, xml, "DRB");
652 664 par->attribute("file", gammalib::xml_file_reduce(xml, m_drbname));
653 665  
  666 + // Set DRW parameter
  667 + par = gammalib::xml_need_par(G_WRITE, xml, "DRW");
  668 + par->attribute("file", gammalib::xml_file_reduce(xml, m_drwname));
  669 +
654 670 // Set DRG parameter
655 671 par = gammalib::xml_need_par(G_WRITE, xml, "DRG");
656 672 par->attribute("file", gammalib::xml_file_reduce(xml, m_drgname));
... ... @@ -702,15 +718,17 @@ void GCOMObservation::write(GXmlElement&amp; xml) const
702 718 *
703 719 * @param[in] drename Event cube name.
704 720 * @param[in] drbname Background cube name.
  721 + * @param[in] drwname Weighting cube name.
705 722 * @param[in] drgname Geometry cube name.
706 723 * @param[in] drxname Exposure map name.
707 724 *
708   - * Load event cube from DRE file, background model from DRB file, geometry
709   - * factors from DRG file and the exposure map from the DRX file. All files
710   - * are mandatory.
  725 + * Load event cube from DRE file, background model from DRB file, weigthing
  726 + * cube from DRW file, geometry factors from DRG file and the exposure map
  727 + * from the DRX file. All files are mandatory.
711 728 ***************************************************************************/
712 729 void GCOMObservation::load(const GFilename& drename,
713 730 const GFilename& drbname,
  731 + const GFilename& drwname,
714 732 const GFilename& drgname,
715 733 const GFilename& drxname)
716 734 {
... ... @@ -720,6 +738,9 @@ void GCOMObservation::load(const GFilename&amp; drename,
720 738 // Load DRB
721 739 load_drb(drbname);
722 740  
  741 + // Load DRW
  742 + load_drw(drwname);
  743 +
723 744 // Load DRG
724 745 load_drg(drgname);
725 746  
... ... @@ -1000,11 +1021,6 @@ std::string GCOMObservation::print(const GChatter&amp; chatter) const
1000 1021 result.append("\n"+gammalib::parformat("BVCs")+"undefined");
1001 1022 }
1002 1023  
1003   - // Append DRB, DRG and DRX
1004   - //result.append("\n"+m_drb.print(chatter));
1005   - //result.append("\n"+m_drg.print(chatter));
1006   - //result.append("\n"+m_drx.print(chatter));
1007   -
1008 1024 } // endif: chatter was not silent
1009 1025  
1010 1026 // Return result
... ... @@ -1034,10 +1050,12 @@ void GCOMObservation::init_members(void)
1034 1050 // Initialise members for binned observation
1035 1051 m_drename.clear();
1036 1052 m_drbname.clear();
  1053 + m_drwname.clear();
1037 1054 m_drgname.clear();
1038 1055 m_drxname.clear();
1039 1056 m_rspname.clear();
1040 1057 m_drb.clear();
  1058 + m_drw.clear();
1041 1059 m_drg.clear();
1042 1060 m_drx.clear();
1043 1061 m_ewidth = 0.0;
... ... @@ -1076,10 +1094,12 @@ void GCOMObservation::copy_members(const GCOMObservation&amp; obs)
1076 1094 // Copy members for binned observation
1077 1095 m_drename = obs.m_drename;
1078 1096 m_drbname = obs.m_drbname;
  1097 + m_drwname = obs.m_drwname;
1079 1098 m_drgname = obs.m_drgname;
1080 1099 m_drxname = obs.m_drxname;
1081 1100 m_rspname = obs.m_rspname;
1082 1101 m_drb = obs.m_drb;
  1102 + m_drw = obs.m_drw;
1083 1103 m_drg = obs.m_drg;
1084 1104 m_drx = obs.m_drx;
1085 1105 m_ewidth = obs.m_ewidth;
... ... @@ -1214,6 +1234,50 @@ void GCOMObservation::load_drb(const GFilename&amp; drbname)
1214 1234  
1215 1235  
1216 1236 /***********************************************************************//**
  1237 + * @brief Load weighting cube from DRW file
  1238 + *
  1239 + * @param[in] drwname DRW filename.
  1240 + *
  1241 + * @exception GException::invalid_value
  1242 + * DRW data space incompatible with DRE data space.
  1243 + *
  1244 + * Load the weighting cube from the primary image of the specified FITS file.
  1245 + ***************************************************************************/
  1246 +void GCOMObservation::load_drw(const GFilename& drwname)
  1247 +{
  1248 + // Continue only if filename is not empty
  1249 + if (!drwname.is_empty()) {
  1250 +
  1251 + // Open FITS file
  1252 + GFits fits(drwname);
  1253 +
  1254 + // Get image
  1255 + const GFitsImage& image = *fits.image("Primary");
  1256 +
  1257 + // Read weighting cube
  1258 + m_drw.read(image);
  1259 +
  1260 + // Close FITS file
  1261 + fits.close();
  1262 +
  1263 + // Check map dimensions
  1264 + if (!check_dri(m_drw)) {
  1265 + std::string msg = "DRW data cube \""+drwname+"\" incompatible with "
  1266 + "DRE data cube \""+m_drename+"\".";
  1267 + throw GException::invalid_value(G_LOAD_DRW, msg);
  1268 + }
  1269 +
  1270 + // Store DRW filename
  1271 + m_drwname = drwname;
  1272 +
  1273 + } // endif: DRW filename was empty
  1274 +
  1275 + // Return
  1276 + return;
  1277 +}
  1278 +
  1279 +
  1280 +/***********************************************************************//**
1217 1281 * @brief Load geometry factors from DRG file
1218 1282 *
1219 1283 * @param[in] drgname DRG filename.
... ... @@ -1245,7 +1309,7 @@ void GCOMObservation::load_drg(const GFilename&amp; drgname)
1245 1309 if (!check_dri(m_drg)) {
1246 1310 std::string msg = "DRG data cube \""+drgname+"\" incompatible with "
1247 1311 "DRE data cube \""+m_drename+"\".";
1248   - throw GException::invalid_value(G_LOAD_DRB, msg);
  1312 + throw GException::invalid_value(G_LOAD_DRG, msg);
1249 1313 }
1250 1314  
1251 1315 // Store DRG filename
... ...