Domenico Tiziani / gammalib

Commit b94ea3386ad20db69cbae8a98e54ba2605209b74

Authored by Jürgen Knödlseder
1 parent 640e93f4

GCTAOnOffObservation now supports any kind of background model (#2710)

Inline Side-by-side
Showing 8 changed files with 119 additions and 85 deletions
ChangeLog
  View file @b94ea33
1   -2018-11-12
  1 +2018-11-16
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
... ... @@ -10,6 +10,7 @@
10 10 to GCTAResponseIrf::lo_safe_thres() and GCTAResponseIrf::hi_safe_thres()
11 11 (#2716)
12 12 Add generic CTA response cache (#2715)
  13 + GCTAOnOffObservation now supports any kind of background model (#2710)
13 14 Fix precision issue with GCTABackground3D Monte-Carlo simulations (#2708)
14 15 Revise GCTAEdisp class interface (#2704)
15 16 GCTAEdisp returns now energy dispersion values per MeV (#2704)
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -12 November 2018
  3 +16 November 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -148,7 +148,7 @@ The following methods have been renamed:
148 148 - GCTAResponseIrf::hi_safe_thres() from GCTAResponseIrf::hi_save_thres()
149 149  
150 150 The arguments for the following methods have been changed:
151   -- GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
  151 +- GCTAOnOffObservation(GCTAObservation&, GModels&, std::string&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
152 152 - GCTAModelRadial::mc(GRan&)
153 153 - GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&)
154 154 - GCTAEdisp::mc(GRan&, GEnergy&, double&, double&, double&, double&)
... ... @@ -339,6 +339,14 @@ parameters. This means that now all other spatial models with frozen parameters
339 339 will be cached. The use of the cache can be toggled by internal data members.
340 340 The IRF cache was removed from GCTAEventList since it is no longer used (#2715).
341 341  
  342 +The GCTAOnOffObservation observation constructor now takes a model container
  343 +as argument and the name of a sky model component. The ARF computation is
  344 +performed by computing the IRF for the spatial model of the named sky model
  345 +component over the On region. In that way also extended sources can be properly
  346 +handled (#2292). The background model is now also extracted from the model
  347 +container and used for the computation of the alpha factors and the background
  348 +response, if the use_model_bkg argument is true (#2710).
  349 +
342 350 GCTABackground3D::init_mc_cache now pre-computes the Monte Carlo cache at energy
343 351 nodes of the background cube, and in addition, inserts 8 logarithmically spaced
344 352 nodes between each pair of nodes. This assures that the energy variations of the
... ... @@ -425,11 +433,6 @@ for the case that the model has no parameters. This is particularily needed for
425 433 the TS computation with the "wstat" statistic which is meaningful without any
426 434 model (#2312).
427 435  
428   -The GCTAOnOffObservation observation constructor now takes a spatial model
429   -component of type GModelSpatial instead of a sky direction as argument, and the
430   -ARF computation is performed by computing the IRF for the spatial model over the
431   -On region. In that way also extended source can be properly handled (#2292).
432   -
433 436 The energy bounds of a CTA observations are now taken into account when setting
434 437 up a GCTAOnOffObservation so that observation-specific energy thresholds are
435 438 correctly handled (#2429).
... ...
doc/source/admin/release_history/1.6.rst
  View file @b94ea33
... ... @@ -47,11 +47,14 @@ Bug fixes
47 47 Improvements
48 48 ------------
49 49  
  50 +
50 51 * [`2721 <https://cta-redmine.irap.omp.eu/issues/2721>`_] -
51 52 Revise ``GCTAEdispCube::operator()`` interface and disable stacked response
52 53 cutting in true energy
53 54 * [`2715 <https://cta-redmine.irap.omp.eu/issues/2715>`_] -
54 55 Add generic CTA response cache
  56 +* [`2710 <https://cta-redmine.irap.omp.eu/issues/2710>`_] -
  57 + ``GCTAOnOffObservation`` now supports any kind of background model
55 58 * [`2704 <https://cta-redmine.irap.omp.eu/issues/2704>`_] -
56 59 ``GCTAEdisp`` returns now energy dispersion values per MeV.
57 60 The ``GCTAEdisp`` interface was revised.
... ...
inst/cta/include/GCTAOnOffObservation.hpp
  View file @b94ea33
... ... @@ -75,12 +75,13 @@ public:
75 75 const GArf& arf,
76 76 const GRmf& rmf);
77 77 GCTAOnOffObservation(const GCTAObservation& obs,
78   - const GModelSpatial& spatial,
  78 + const GModels& models,
  79 + const std::string& srcname,
79 80 const GEbounds& etrue,
80 81 const GEbounds& ereco,
81 82 const GSkyRegions& on,
82 83 const GSkyRegions& off,
83   - const bool& use_irf_bkg = true);
  84 + const bool& use_model_bkg = true);
84 85 GCTAOnOffObservation(const GCTAOnOffObservation& obs);
85 86 virtual ~GCTAOnOffObservation(void);
86 87  
... ... @@ -128,10 +129,11 @@ protected:
128 129 GArf arf_stacked(const GArf& arf, const GEnergy& emin, const GEnergy& emax) const;
129 130 GRmf rmf_stacked(const GRmf& rmf, const GEnergy& emin, const GEnergy& emax) const;
130 131 void set(const GCTAObservation& obs,
131   - const GModelSpatial& spatial,
  132 + const GModels& models,
  133 + const std::string& srcname,
132 134 const GSkyRegions& on,
133 135 const GSkyRegions& off,
134   - const bool& use_irf_bkg);
  136 + const bool& use_model_bkg);
135 137 void compute_arf(const GCTAObservation& obs,
136 138 const GModelSpatial& spatial,
137 139 const GSkyRegions& on);
... ... @@ -140,11 +142,13 @@ protected:
140 142 const GSkyRegions& on);
141 143 void compute_bgd(const GCTAObservation& obs,
142 144 const GSkyRegions& off,
143   - const bool& use_irf_bkg);
  145 + const GModels& models,
  146 + const bool& use_model_bkg);
144 147 void compute_alpha(const GCTAObservation& obs,
145 148 const GSkyRegions& on,
146 149 const GSkyRegions& off,
147   - const bool& use_irf_bkg);
  150 + const GModels& models,
  151 + const bool& use_model_bkg);
148 152 void compute_rmf(const GCTAObservation& obs,
149 153 const GSkyRegions& on);
150 154 void apply_ebounds(const GCTAObservation& obs);
... ...
inst/cta/pyext/GCTAOnOffObservation.i
  View file @b94ea33
... ... @@ -82,12 +82,13 @@ public:
82 82 const GArf& arf,
83 83 const GRmf& rmf);
84 84 GCTAOnOffObservation(const GCTAObservation& obs,
85   - const GModelSpatial& spatial,
  85 + const GModels& models,
  86 + const std::string& srcname,
86 87 const GEbounds& etrue,
87 88 const GEbounds& ereco,
88 89 const GSkyRegions& on,
89 90 const GSkyRegions& off,
90   - const bool& use_irf_bkg = true);
  91 + const bool& use_bkg_model = true);
91 92 GCTAOnOffObservation(const GCTAOnOffObservation& obs);
92 93 virtual ~GCTAOnOffObservation(void);
93 94  
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @b94ea33
... ... @@ -360,9 +360,13 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
360 360 // Compute value
361 361 double value = spat * spec * temp;
362 362  
363   - // Apply deadtime correction
  363 + // Apply deadtime correction. This assumes that the IRF background rates
  364 + // are stored per livetime, and by multiplying these rates with the
  365 + // deadtime correction factor thay are given per ontime.
  366 + /*
364 367 double deadc = obs.deadc(event.time());
365 368 value *= deadc;
  369 + */
366 370  
367 371 // If gradients were requested then multiply factors to spectral and
368 372 // temporal gradients
... ... @@ -370,7 +374,7 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
370 374  
371 375 // Multiply factors to spectral gradients
372 376 if (spectral() != NULL) {
373   - double fact = spat * temp * deadc;
  377 + double fact = spat * temp; // * deadc;
374 378 if (fact != 1.0) {
375 379 for (int i = 0; i < spectral()->size(); ++i)
376 380 (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
... ... @@ -379,7 +383,7 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
379 383  
380 384 // Multiply factors to temporal gradients
381 385 if (temporal() != NULL) {
382   - double fact = spat * spec * deadc;
  386 + double fact = spat * spec; // * deadc;
383 387 if (fact != 1.0) {
384 388 for (int i = 0; i < temporal()->size(); ++i)
385 389 (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @b94ea33
... ... @@ -183,12 +183,13 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GPha&amp; pha_on,
183 183 * @brief CTA observation constructor
184 184 *
185 185 * @param[in] obs CTA observation.
186   - * @param[in] spatial Spatial source model.
  186 + * @param[in] models Models including source and background model.
  187 + * @param[in] srcname Name of source in models.
187 188 * @param[in] etrue True energy boundaries.
188 189 * @param[in] ereco Reconstructed energy boundaries.
189 190 * @param[in] on On regions.
190 191 * @param[in] off Off regions.
191   - * @param[in] use_irf_bkg Use IRF background template.
  192 + * @param[in] use_model_bkg Use model background.
192 193 *
193 194 * Constructs On/Off observation by filling the On and Off spectra and
194 195 * computing the Auxiliary Response File (ARF) and Redistribution Matrix
... ... @@ -196,19 +197,20 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GPha&amp; pha_on,
196 197 * reconstructed energy boundaries, as well as the definition of On and Off
197 198 * regions.
198 199 *
199   - * The @p use_irf_bkg flags controls whether the IRF background template
200   - * should be used for computations or not. This impacts the computations of
201   - * the @c BACKSCAL column in the On spectrum and the @c BACKRESP column in
  200 + * The @p use_model_bkg flag controls whether a background model should be
  201 + * used for the computations or not. This impacts the computations of the
  202 + * @c BACKSCAL column in the On spectrum and the @c BACKRESP column in
202 203 * the Off spectrum. See the compute_alpha() and compute_bgd() methods for
203 204 * more information on the applied formulae.
204 205 ***************************************************************************/
205 206 GCTAOnOffObservation::GCTAOnOffObservation(const GCTAObservation& obs,
206   - const GModelSpatial& spatial,
  207 + const GModels& models,
  208 + const std::string& srcname,
207 209 const GEbounds& etrue,
208 210 const GEbounds& ereco,
209 211 const GSkyRegions& on,
210 212 const GSkyRegions& off,
211   - const bool& use_irf_bkg)
  213 + const bool& use_model_bkg)
212 214 {
213 215 // Initialise private
214 216 init_members();
... ... @@ -222,7 +224,7 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GCTAObservation&amp; obs,
222 224 m_rmf = GRmf(etrue, ereco);
223 225  
224 226 // Set On/Off observation from CTA observation
225   - set(obs, spatial, on, off, use_irf_bkg);
  227 + set(obs, models, srcname, on, off, use_model_bkg);
226 228  
227 229 // Return
228 230 return;
... ... @@ -320,8 +322,10 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations&amp; obs) :
320 322 GEnergy emax = onoff->on_spec().emax_obs();
321 323  
322 324 // Get stacked ARF and RMF
323   - GArf arf_stacked = this->arf_stacked(onoff->arf(), emin, emax);
324   - GRmf rmf_stacked = this->rmf_stacked(onoff->rmf(), emin, emax);
  325 + //GArf arf_stacked = this->arf_stacked(onoff->arf(), emin, emax);
  326 + //GRmf rmf_stacked = this->rmf_stacked(onoff->rmf(), emin, emax);
  327 + GArf arf_stacked = onoff->arf();
  328 + GRmf rmf_stacked = onoff->rmf();
325 329  
326 330 // Get exposure for observation
327 331 double exposure = onoff->on_spec().exposure();
... ... @@ -1195,10 +1199,11 @@ GRmf GCTAOnOffObservation::rmf_stacked(const GRmf&amp; rmf,
1195 1199 * @brief Set On/Off observation from a CTA observation
1196 1200 *
1197 1201 * @param[in] obs CTA observation.
1198   - * @param[in] spatial Spatial source model.
  1202 + * @param[in] models Models including source and background model.
  1203 + * @param[in] srcname Name of soucre in models.
1199 1204 * @param[in] on On regions.
1200 1205 * @param[in] off Off regions.
1201   - * @param[in] use_irf_bkg Use IRF background template.
  1206 + * @param[in] use_model_bkg Use model background.
1202 1207 *
1203 1208 * @exception GException::invalid_value
1204 1209 * No CTA event list found in CTA observation.
... ... @@ -1207,17 +1212,18 @@ GRmf GCTAOnOffObservation::rmf_stacked(const GRmf&amp; rmf,
1207 1212 * that fall in the On and Off regions into the PHA spectra and by computing
1208 1213 * the corresponding ARF and RMF response functions.
1209 1214 *
1210   - * The @p use_irf_bkg flags controls whether the IRF background template
1211   - * should be used for computations or not. This impacts the computations of
1212   - * the `BACKSCAL` column in the On spectrum and the `BACKRESP` column in
1213   - * the Off spectrum. See the compute_alpha() and compute_bgd() methods for
1214   - * more information on the applied formulae.
  1215 + * The @p use_model_bkg flags controls whether the background model should
  1216 + * be used for computations or not. This impacts the computations of the
  1217 + * `BACKSCAL` column in the On spectrum and the `BACKRESP` column in the Off
  1218 + * spectrum. See the compute_alpha() and compute_bgd() methods for more
  1219 + * information on the applied formulae.
1215 1220 ***************************************************************************/
1216 1221 void GCTAOnOffObservation::set(const GCTAObservation& obs,
1217   - const GModelSpatial& spatial,
  1222 + const GModels& models,
  1223 + const std::string& srcname,
1218 1224 const GSkyRegions& on,
1219 1225 const GSkyRegions& off,
1220   - const bool& use_irf_bkg)
  1226 + const bool& use_model_bkg)
1221 1227 {
1222 1228 // Get CTA event list pointer
1223 1229 const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
... ... @@ -1228,6 +1234,27 @@ void GCTAOnOffObservation::set(const GCTAObservation&amp; obs,
1228 1234 throw GException::invalid_value(G_SET, msg);
1229 1235 }
1230 1236  
  1237 + // Get spatial component of source model
  1238 + const GModelSky* model = dynamic_cast<const GModelSky*>(models[srcname]);
  1239 + if (model == NULL) {
  1240 + std::string msg = "Model \""+srcname+"\" is not a sky model. Please "
  1241 + "specify the name of a sky model.";
  1242 + throw GException::invalid_value(G_SET, msg);
  1243 + }
  1244 + const GModelSpatial& spatial = *(model->spatial());
  1245 +
  1246 + // If background model is needed then extract background model components
  1247 + // from model container
  1248 + GModels bkg_models;
  1249 + if (use_model_bkg) {
  1250 + for (int i = 0; i < models.size(); ++i) {
  1251 + if ((dynamic_cast<const GModelData*>(models[i]) != NULL) &&
  1252 + (models[i]->classname().substr(0,4) == "GCTA")) {
  1253 + bkg_models.append(*models[i]);
  1254 + }
  1255 + }
  1256 + }
  1257 +
1231 1258 // Loop over all events
1232 1259 for (int i = 0; i < events->size(); ++i) {
1233 1260  
... ... @@ -1274,8 +1301,8 @@ void GCTAOnOffObservation::set(const GCTAObservation&amp; obs,
1274 1301 else {
1275 1302 compute_arf(obs, spatial, reg_on);
1276 1303 }
1277   - compute_bgd(obs, reg_off, use_irf_bkg);
1278   - compute_alpha(obs, reg_on, reg_off, use_irf_bkg);
  1304 + compute_bgd(obs, reg_off, bkg_models, use_model_bkg);
  1305 + compute_alpha(obs, reg_on, reg_off, bkg_models, use_model_bkg);
1279 1306 compute_rmf(obs, reg_on);
1280 1307  
1281 1308 // Apply reconstructed energy boundaries
... ... @@ -1538,11 +1565,8 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1538 1565 *
1539 1566 * @param[in] obs CTA observation.
1540 1567 * @param[in] off Off regions.
1541   - * @param[in] use_irf_bkg Use IRF background template.
1542   - *
1543   - * @exception GException::invalid_argument
1544   - * Observation does not contain relevant response or background
1545   - * information
  1568 + * @param[in] models Models including background model.
  1569 + * @param[in] use_model_bkg Use model background.
1546 1570 *
1547 1571 * Computes the background rate in units of
1548 1572 * \f${\rm events} \, {\rm MeV}^{-1} \, {\rm s}^{-1}\f$
... ... @@ -1552,7 +1576,7 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1552 1576 * All Off regions contained in @p off are expected to be sky region maps,
1553 1577 * i.e. of type GSkyRegionMap.
1554 1578 *
1555   - * If @p use_irf_bkg is @c true, the IRF background template will be used
  1579 + * If @p use_model_bkg is @c true, the IRF background template will be used
1556 1580 * for the computation, and `BACKRESP` will be computed using
1557 1581 *
1558 1582 * \f[
... ... @@ -1572,7 +1596,7 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1572 1596 * \f$\Omega_{{\rm off},p}\f$ is the solid angle in units of \f${\rm sr}\f$
1573 1597 * of the pixel index \f$p\f$ in the Off region \f${\rm off}\f$.
1574 1598 *
1575   - * If @p use_irf_bkg is @c false, `BACKRESP` will be computed using
  1599 + * If @p use_model_bkg is @c false, `BACKRESP` will be computed using
1576 1600 *
1577 1601 * \f[
1578 1602 * {\tt BACKRESP}(E_{\rm reco}) = \sum_{\rm off} \sum_p \Omega_{{\rm off},p}
... ... @@ -1586,7 +1610,8 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1586 1610 ***************************************************************************/
1587 1611 void GCTAOnOffObservation::compute_bgd(const GCTAObservation& obs,
1588 1612 const GSkyRegions& off,
1589   - const bool& use_irf_bkg)
  1613 + const GModels& models,
  1614 + const bool& use_model_bkg)
1590 1615 {
1591 1616 // Get reconstructed energies for the background rates
1592 1617 const GEbounds& ereco = m_off_spec.ebounds();
... ... @@ -1602,10 +1627,7 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1602 1627 GCTAPointing obspnt = obs.pointing();
1603 1628  
1604 1629 // Continue only if the IRF background template should be used
1605   - if (use_irf_bkg) {
1606   -
1607   - // Get CTA background
1608   - const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_COMPUTE_BGD, obs);
  1630 + if (use_model_bkg) {
1609 1631  
1610 1632 // Loop over regions
1611 1633 for (int k = 0; k < off.size(); ++k) {
... ... @@ -1633,13 +1655,11 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1633 1655 // Loop over energy bins
1634 1656 for (int i = 0; i < nreco; ++i) {
1635 1657  
1636   - // Get log10(E/TeV) of mean reconstructed bin energy
1637   - double logEreco = ereco.elogmean(i).log10TeV();
  1658 + // Set event
  1659 + GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
1638 1660  
1639 1661 // Get background rate in events/s/MeV
1640   - background[i] += bgd(logEreco,
1641   - pixinstdir.detx(),
1642   - pixinstdir.dety()) * pixsolid;
  1662 + background[i] += models.eval(event, obs) * pixsolid;
1643 1663  
1644 1664 } // endfor: looped over energy bins
1645 1665  
... ... @@ -1705,7 +1725,8 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1705 1725 * @param[in] obs CTA observation.
1706 1726 * @param[in] on On regions.
1707 1727 * @param[in] off Off regions.
1708   - * @param[in] use_irf_bkg Use IRF background template.
  1728 + * @param[in] models Models including background model.
  1729 + * @param[in] use_model_bkg Use model background.
1709 1730 *
1710 1731 * @exception GException::invalid_argument
1711 1732 * Observation does not contain relevant response or background
... ... @@ -1715,7 +1736,7 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1715 1736 * The \f$\alpha\f$ parameter gives the ratio between the On and Off region
1716 1737 * background acceptance multiplied by the On and Off region solid angles.
1717 1738 *
1718   - * If @p use_irf_bkg is @c true, the IRF background template will be used
  1739 + * If @p use_model_bkg is @c true, the IRF background template will be used
1719 1740 * for the computation, and \f$\alpha(E_{\rm reco})\f$ is given by
1720 1741 *
1721 1742 * \f[
... ... @@ -1740,7 +1761,7 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1740 1761 * \f${\rm sr}\f$ of the pixel index \f$p\f$ in the On or Off region
1741 1762 * \f${\rm on/off}\f$.
1742 1763 *
1743   - * If @p use_irf_bkg is @c false, the background acceptance is assumed
  1764 + * If @p use_model_bkg is @c false, the background acceptance is assumed
1744 1765 * constant and hence cancels out, which makes the \f$\alpha\f$ parameter
1745 1766 * independent of reconstructed energy \f$E_{\rm reco}\f$.
1746 1767 * The \f$\alpha\f$ parameter is then given by
... ... @@ -1756,7 +1777,8 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1756 1777 void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs,
1757 1778 const GSkyRegions& on,
1758 1779 const GSkyRegions& off,
1759   - const bool& use_irf_bkg)
  1780 + const GModels& models,
  1781 + const bool& use_model_bkg)
1760 1782 {
1761 1783 // Get reconstructed energy boundaries from RMF
1762 1784 const GEbounds& ereco = m_rmf.emeasured();
... ... @@ -1770,16 +1792,13 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs,
1770 1792  
1771 1793 // If IRF background templates shall be used then compute the
1772 1794 // energy dependent alpha factors
1773   - if (use_irf_bkg) {
1774   -
1775   - // Get CTA background
1776   - const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_COMPUTE_ALPHA, obs);
  1795 + if (use_model_bkg) {
1777 1796  
1778 1797 // Loop over reconstructed energies
1779 1798 for (int i = 0; i < nreco; ++i) {
1780 1799  
1781 1800 // Get mean log10 energy in TeV of bin
1782   - double logEreco = ereco.elogmean(i).log10TeV();
  1801 + //double logEreco = ereco.elogmean(i).log10TeV();
1783 1802  
1784 1803 // Initialise background rate totals
1785 1804 double aon = 0.0;
... ... @@ -1807,11 +1826,11 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs,
1807 1826 // Get solid angle subtended by this pixel
1808 1827 double pixsolid = on_map->map().solidangle(pixidx);
1809 1828  
  1829 + // Set event
  1830 + GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
  1831 +
1810 1832 // Add up acceptance
1811   - aon += bgd(logEreco,
1812   - pixinstdir.detx(),
1813   - pixinstdir.dety()) *
1814   - pixsolid;
  1833 + aon += models.eval(event, obs) * pixsolid;
1815 1834  
1816 1835 } // endfor: looped over all pixels in map
1817 1836  
... ... @@ -1839,11 +1858,11 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs,
1839 1858 // Get solid angle subtended by this pixel
1840 1859 double pixsolid = off_map->map().solidangle(pixidx);
1841 1860  
  1861 + // Set event
  1862 + GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
  1863 +
1842 1864 // Add up acceptance
1843   - aoff += bgd(logEreco,
1844   - pixinstdir.detx(),
1845   - pixinstdir.dety()) *
1846   - pixsolid;
  1865 + aoff += models.eval(event, obs) * pixsolid;
1847 1866  
1848 1867 } // endfor: looped over all pixels in map
1849 1868  
... ...
inst/cta/test/test_CTA.py
  View file @b94ea33
... ... @@ -305,16 +305,14 @@ class Test(gammalib.GPythonTestSuite):
305 305 """
306 306 Test On/Off analysis
307 307 """
308   - # Set Crab position
309   - crab = gammalib.GSkyDir()
310   - crab.radec_deg(83.6331, 22.0145)
311   -
312   - # Set point-source spatial model based on Crab position
313   - spatial = gammalib.GModelSpatialPointSource(crab)
  308 + # Load model container
  309 + models = gammalib.GModels(self._data + '/crab_irf.xml')
314 310  
315 311 # Create On region
  312 + ondir = gammalib.GSkyDir()
  313 + ondir.radec_deg(83.6331, 22.0145)
316 314 on = gammalib.GSkyRegions()
317   - on.append(gammalib.GSkyRegionCircle(crab, 0.2))
  315 + on.append(gammalib.GSkyRegionCircle(ondir, 0.2))
318 316  
319 317 # Create Off region
320 318 offdir = gammalib.GSkyDir()
... ... @@ -333,10 +331,11 @@ class Test(gammalib.GPythonTestSuite):
333 331 inobs = gammalib.GObservations(filename)
334 332 outobs = gammalib.GObservations()
335 333 for run in inobs:
336   - onoff = gammalib.GCTAOnOffObservation(run, spatial, etrue, ereco, on, off)
  334 + onoff = gammalib.GCTAOnOffObservation(run, models, 'Crab',
  335 + etrue, ereco, on, off)
337 336 outobs.append(onoff)
338 337  
339   - # Load model container and attach it to the observations
  338 + # Load On/Off models and attach them to the observations
340 339 models = gammalib.GModels(self._data + '/onoff_model.xml')
341 340 outobs.models(models)
342 341  
... ... @@ -454,7 +453,7 @@ class Test(gammalib.GPythonTestSuite):
454 453 region = gammalib.GSkyRegionCircle(dir, 0.2)
455 454 regs = gammalib.GSkyRegions()
456 455 regs.append(region)
457   - ptsrc = gammalib.GModelSpatialPointSource()
  456 + models = gammalib.GModels(self._data + '/crab_irf.xml')
458 457 expcube = gammalib.GCTACubeExposure('CAR','CEL',0.,0.,0.1,0.1,10,10,engs)
459 458 psfcube = gammalib.GCTACubePsf('CAR','CEL',0.,0.,0.1,0.1,10,10,engs,1.0,10)
460 459 bgdcube = gammalib.GCTACubeBackground('CAR','CEL',0.,0.,0.1,0.1,10,10,ebds)
... ... @@ -476,7 +475,7 @@ class Test(gammalib.GPythonTestSuite):
476 475 obs4 = gammalib.GCTAObservation()
477 476 obs4.events(cube)
478 477 obs4.response(rspcube2)
479   - obs5 = gammalib.GCTAOnOffObservation(obs1, ptsrc, ebds, ebds, regs, regs)
  478 + obs5 = gammalib.GCTAOnOffObservation(obs1, models, 'Crab', ebds, ebds, regs, regs)
480 479 radgauss = gammalib.GCTAModelRadialGauss(1.0)
481 480 radacc = gammalib.GCTAModelRadialAcceptance(radgauss, plaw)
482 481 multi = gammalib.GCTAModelSpatialMultiplicative()
... ...