GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAOnOffObservation.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAOnOffObservation.cpp - CTA On/Off observation class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2022 by Chia-Chun Lu & Christoph Deil *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAOnOffObservation.cpp
23  * @brief CTA On/Off observation class implementation
24  * @author Chia-Chun Lu & Christoph Deil & Pierrick Martin
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <typeinfo>
32 #include "GObservationRegistry.hpp"
33 #include "GTools.hpp"
34 #include "GIntegral.hpp"
35 #include "GMatrixSparse.hpp"
36 #include "GModels.hpp"
37 #include "GModelSky.hpp"
38 #include "GModelData.hpp"
39 #include "GModelSpatial.hpp"
40 #include "GModelSpectral.hpp"
41 #include "GModelTemporal.hpp"
42 #include "GSource.hpp"
43 #include "GSkyRegions.hpp"
44 #include "GSkyRegionMap.hpp"
45 #include "GSkyRegionCircle.hpp"
46 #include "GOptimizerPars.hpp"
47 #include "GObservations.hpp"
48 #include "GCTAObservation.hpp"
49 #include "GCTAEventAtom.hpp"
50 #include "GCTAEventCube.hpp"
51 #include "GCTAResponseIrf.hpp"
52 #include "GCTAAeff2D.hpp"
53 #include "GCTACubeBackground.hpp"
54 #include "GCTAOnOffObservation.hpp"
55 #include "GCTASupport.hpp"
56 
57 /* __ OpenMP section _____________________________________________________ */
58 #ifdef _OPENMP
59 #include <omp.h>
60 #endif
61 
62 /* __ Globals ____________________________________________________________ */
63 const GCTAOnOffObservation g_onoff_obs_cta_seed(true, "CTAOnOff");
64 const GCTAOnOffObservation g_onoff_obs_hess_seed(true, "HESSOnOff");
65 const GCTAOnOffObservation g_onoff_obs_magic_seed(true, "MAGICOnOff");
66 const GCTAOnOffObservation g_onoff_obs_veritas_seed(true, "VERITASOnOff");
67 const GCTAOnOffObservation g_onoff_obs_astri_seed(true, "ASTRIOnOff");
68 const GCTAOnOffObservation g_onoff_obs_fact_seed(true, "FACTOnOff");
69 const GObservationRegistry g_onoff_obs_cta_registry(&g_onoff_obs_cta_seed);
70 const GObservationRegistry g_onoff_obs_hess_registry(&g_onoff_obs_hess_seed);
71 const GObservationRegistry g_onoff_obs_magic_registry(&g_onoff_obs_magic_seed);
72 const GObservationRegistry g_onoff_obs_veritas_registry(&g_onoff_obs_veritas_seed);
73 const GObservationRegistry g_onoff_obs_astri_registry(&g_onoff_obs_astri_seed);
74 const GObservationRegistry g_onoff_obs_fact_registry(&g_onoff_obs_fact_seed);
75 
76 /* __ Method name definitions ____________________________________________ */
77 #define G_CONSTRUCTOR1 "GCTAOnOffObservation::GCTAOnOffObservation(GPha&, "\
78  "GPha&, GArf&, GRmf&)"
79 #define G_CONSTRUCTOR2 "GCTAOnOffObservation(GObservations& obs)"
80 #define G_RESPONSE_SET "GCTAOnOffObservation::response(GResponse&)"
81 #define G_RESPONSE_GET "GCTAOnOffObservation::response()"
82 #define G_WRITE "GCTAOnOffObservation::write(GXmlElement&)"
83 #define G_READ "GCTAOnOffObservation::read(GXmlElement&)"
84 #define G_LIKELIHOOD "GCTAOnOffObservation::likelihood(GModels&, "\
85  "GOptimizerPars&, GMatrixSparse&, "\
86  "GVector&, double&, double&)"
87 #define G_SET "GCTAOnOffObservation::set(GCTAObservation&, GModelSpatial&, "\
88  "GSkyRegionMap&, GSkyRegionMap&)"
89 #define G_COMPUTE_ARF "GCTAOnOffObservation::compute_arf(GCTAObservation&, "\
90  "GModelSpatial&, GSkyRegionMap&)"
91 #define G_COMPUTE_ARF_CUT "GCTAOnOffObservation::compute_arf_cut("\
92  "GCTAObservation&, GModelSpatial&, GSkyRegionMap&)"
93 #define G_COMPUTE_BGD "GCTAOnOffObservation::compute_bgd(GCTAObservation&, "\
94  "GSkyRegionMap&)"
95 #define G_COMPUTE_ALPHA "GCTAOnOffObservation::compute_alpha("\
96  "GCTAObservation&, GSkyRegionMap&, GSkyRegionMap&)"
97 #define G_COMPUTE_RMF "GCTAOnOffObservation::compute_rmf(GCTAObservation&, "\
98  "GSkyRegionMap&)"
99 #define G_MODEL_BACKGROUND "GCTAOnOffObservation::model_background(GModels&)"
100 #define G_ARF_RADIUS_CUT "GCTAOnOffObservation::arf_radius_cut("\
101  "GCTAObservation&, GSkyRegions&)"
102 
103 /* __ Constants __________________________________________________________ */
104 const double minmod = 1.0e-100; //!< Minimum model value
105 const double minerr = 1.0e-100; //!< Minimum statistical error
106 
107 /* __ Macros _____________________________________________________________ */
108 
109 /* __ Coding definitions _________________________________________________ */
110 
111 /* __ Debug definitions __________________________________________________ */
112 //#define G_LIKELIHOOD_DEBUG //!< Debug likelihood computation
113 //#define G_N_GAMMA_DEBUG //!< Debug N_gamma computation
114 
115 
116 /*==========================================================================
117  = =
118  = Constructors/destructors =
119  = =
120  ==========================================================================*/
121 
122 /***********************************************************************//**
123  * @brief Void constructor
124  *
125  * Constructs empty On/Off observation.
126  ***************************************************************************/
128 {
129  // Initialise private members
130  init_members();
131 
132  // Return
133  return;
134 }
135 
136 
137 /***********************************************************************//**
138  * @brief Instrument constructor
139  *
140  * @param[in] dummy Dummy flag.
141  * @param[in] instrument Instrument name.
142  *
143  * Constructs an empty CTA On/Off observation for a given instrument.
144  *
145  * This method is specifically used allocation of global class instances for
146  * observation registry. By specifying explicit instrument names it is
147  * possible to use the "CTA" module are for other Imaging Air Cherenkov
148  * Telescopes. So far, the following instrument codes are supported:
149  * "CTAOnOff", "HESSOnOff", "VERITASOnOff", "MAGICOnOff", "ASTRIOnOff".
150  ***************************************************************************/
152  const std::string& instrument) :
153  GObservation()
154 {
155  // Initialise members
156  init_members();
157 
158  // Set instrument
160 
161  // Return
162  return;
163 }
164 
165 
166 /***********************************************************************//**
167  * @brief Copy constructor
168  *
169  * @param[in] obs On/Off observation.
170  ***************************************************************************/
172  GObservation(obs)
173 {
174  // Initialise private
175  init_members();
176 
177  // Copy members
178  copy_members(obs);
179 
180  // Return
181  return;
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief CTA observation constructor
187  *
188  * @param[in] pha_on On spectrum.
189  * @param[in] pha_off Off spectrum.
190  * @param[in] arf Auxiliary Response File.
191  * @param[in] rmf Redistribution Matrix File.
192  *
193  * Constructs On/Off observation from On and Off spectra, an Auxiliary
194  * Response File and a Redistribution Matrix File.
195  ***************************************************************************/
197  const GPha& pha_off,
198  const GArf& arf,
199  const GRmf& rmf)
200 {
201  // Initialise private
202  init_members();
203 
204  // Set data members
205  m_on_spec = pha_on;
206  m_off_spec = pha_off;
207  m_arf = arf;
208  m_rmf = rmf;
209 
210  // Set the ontime, livetime and deadtime correction
211  set_exposure();
212 
213  // Check consistency of On/Off observation
215 
216  // Return
217  return;
218 }
219 
220 
221 /***********************************************************************//**
222  * @brief CTA observation constructor (same On and Off observation)
223  *
224  * @param[in] obs CTA observation.
225  * @param[in] models Models including source and background model.
226  * @param[in] srcname Name of source in models.
227  * @param[in] etrue True energy boundaries.
228  * @param[in] ereco Reconstructed energy boundaries.
229  * @param[in] on On regions.
230  * @param[in] off Off regions.
231  * @param[in] use_model_bkg Use model background.
232  *
233  * Constructs On/Off observation from one osbervation by filling the On and
234  * Off spectra and computing the Auxiliary Response File (ARF) and
235  * Redistribution Matrix File (RMF).
236  *
237  * The method requires the specification of the true and reconstructed energy
238  * boundaries, as well as the definition of On and Off regions.
239  *
240  * The @p use_model_bkg flag controls whether a background model should be
241  * used for the computations or not. This impacts the computations of the
242  * @c BACKSCAL column in the On spectrum and the @c BACKRESP column in
243  * the Off spectrum. See the compute_alpha() and compute_bgd() methods for
244  * more information on the applied formulae.
245  ***************************************************************************/
247  const GModels& models,
248  const std::string& srcname,
249  const GEbounds& etrue,
250  const GEbounds& ereco,
251  const GSkyRegions& on,
252  const GSkyRegions& off,
253  const bool& use_model_bkg)
254 {
255  // Initialise private
256  init_members();
257 
258  // Initialise spectra
259  m_on_spec = GPha(ereco);
260  m_off_spec = GPha(ereco);
261 
262  // Initialise response information
263  m_arf = GArf(etrue);
264  m_rmf = GRmf(etrue, ereco);
265 
266  // Set On/Off observation from CTA observation
267  // In this case On and Off are the same observation
268  set(obs, obs, models, srcname, on, off, use_model_bkg);
269 
270  // Return
271  return;
272 }
273 
274 
275 /***********************************************************************//**
276  * @brief CTA observation constructor (different On and Off observations)
277  *
278  * @param[in] obs_on On CTA observation.
279  * @param[in] obs_off Off CTA observation.
280  * @param[in] models Models including source and background models.
281  * @param[in] srcname Name of source in models.
282  * @param[in] etrue True energy boundaries.
283  * @param[in] ereco Reconstructed energy boundaries.
284  * @param[in] on On regions.
285  * @param[in] off Off regions.
286  * @param[in] use_model_bkg Use model background.
287  *
288  * Constructs On/Off observation from two osbervations by filling the On and
289  * Off spectra and computing the Auxiliary Response File (ARF) and
290  * Redistribution Matrix File (RMF).
291  *
292  * The method requires the specification of the true and reconstructed energy
293  * boundaries, as well as the definition of On and Off regions.
294  *
295  * The @p use_model_bkg flag controls whether a background model should be
296  * used for the computations or not. This impacts the computations of the
297  * @c BACKSCAL column in the On spectrum and the @c BACKRESP column in
298  * the Off spectrum. See the compute_alpha() and compute_bgd() methods for
299  * more information on the applied formulae.
300  ***************************************************************************/
302  const GCTAObservation& obs_off,
303  const GModels& models,
304  const std::string& srcname,
305  const GEbounds& etrue,
306  const GEbounds& ereco,
307  const GSkyRegions& on,
308  const GSkyRegions& off,
309  const bool& use_model_bkg)
310 {
311  // Initialise private
312  init_members();
313 
314  // Initialise spectra
315  m_on_spec = GPha(ereco);
316  m_off_spec = GPha(ereco);
317 
318  // Initialise response information
319  m_arf = GArf(etrue);
320  m_rmf = GRmf(etrue, ereco);
321 
322  // Set On/Off observation from CTA observation
323  // In this case On and Off are the same observation
324  set(obs_on, obs_off, models, srcname, on, off, use_model_bkg);
325 
326  // Return
327  return;
328 }
329 
330 
331 /***********************************************************************//**
332  * @brief CTA On/Off observation stacking constructor
333  *
334  * @param[in] obs Observation container.
335  *
336  * GException::invalid_value
337  * Incompatible On/Off observation definition.
338  *
339  * Constructs On/Off observation by stacking all On/Off observation in the
340  * observation container into a single On/Off observation.
341  *
342  * The constructor uses the following formulae:
343  *
344  * \f[
345  * N^{\rm on}(E_{\rm reco}) = \sum_i N^{\rm on}_i(E_{\rm reco})
346  * \f]
347  *
348  * \f[
349  * N^{\rm off}(E_{\rm reco}) = \sum_i N^{\rm off}_i(E_{\rm reco})
350  * \f]
351  *
352  * \f[
353  * \alpha(E_{\rm reco}) = \frac{\sum_i \alpha_i(E_{\rm reco})
354  * B_i(E_{\rm reco}) \tau_off_i}
355  * {\sum_i B_i(E_{\rm reco}) \tau_off_i}
356  * \f]
357  *
358  * \f[
359  * B(E_{\rm reco}) = \frac{\sum_i B_i(E_{\rm reco}) \tau_off_i}{\sum_i \tau_off_i}
360  * \f]
361  *
362  * \f[
363  * ARF(E_{\rm true}) = \frac{\sum_i ARF_i(E_{\rm true}) \tau_on_i}
364  * {\sum_i \tau_i}
365  * \f]
366  *
367  * \f[
368  * RMF(E_{\rm true}, E_{\rm reco}) =
369  * \frac{\sum_i RMF_i(E_{\rm true}, E_{\rm reco}) ARF_i(E_{\rm true}) \tau_on_i}
370  * {\sum_i ARF_i(E_{\rm true}) \tau_on_i}
371  * \f]
372  *
373  * where
374  * \f$N^{\rm on}_i(E_{\rm reco})\f$ is the On spectrum of observation \f$i\f$,
375  * \f$N^{\rm off}_i(E_{\rm reco})\f$ is the Off spectrum of observation
376  * \f$i\f$,
377  * \f$\alpha_i(E_{\rm reco})\f$ is the background scaling of observation
378  * \f$i\f$,
379  * \f$B_i(E_{\rm reco})\f$ is the background rate for observation \f$i\f$,
380  * \f$ARF_i(E_{\rm true})\f$ is the Auxiliary Response File of observation
381  * \f$i\f$,
382  * \f$RMF_i(E_{\rm true}, E_{\rm reco})\f$ is the Redistribution Matrix File
383  * of observation \f$i\f$, and
384  * \f$\tau_on_i\f$ is the livetime of On observation \f$i\f$.
385  * \f$\tau_off_i\f$ is the livetime of Off observation \f$i\f$.
386  *
387  * The method extracts the instrument name from the first On/Off observation
388  * in the observation container and only stacks subsequent On/Off
389  * observations with the same instrument name.
390  ***************************************************************************/
392  GObservation()
393 {
394  // Initialise private
395  init_members();
396 
397  // Signal first On/Off observation
398  bool first = true;
399 
400  // Initialise total exposure
401  double total_exposure_on = 0.0;
402  double total_exposure_off = 0.0;
403 
404  // Allocate observation energy range
405  GEnergy emin_obs;
406  GEnergy emax_obs;
407 
408  // Allocate instrument
409  std::string instrument;
410 
411  // Loop over all observation in container
412  for (int i = 0; i < obs.size(); ++i) {
413 
414  // Get pointer to On/Off observation
415  const GCTAOnOffObservation* onoff =
416  dynamic_cast<const GCTAOnOffObservation*>(obs[i]);
417 
418  // Skip observation if it is not a On/Off observation
419  if (onoff == NULL) {
420  continue;
421  }
422 
423  // If the instrument name is empty then set the instrument now.
424  // Otherwise, skip the observation if it does not correspond
425  // to the same instrument.
426  if (instrument.empty()) {
427  instrument = onoff->instrument();
428  }
429  else if (instrument != onoff->instrument()) {
430  continue;
431  }
432 
433  // Check consistency of On/Off observation
435 
436  // Get energy boundaries of observation
437  GEnergy emin = onoff->on_spec().emin_obs();
438  GEnergy emax = onoff->on_spec().emax_obs();
439 
440  // Get stacked ARF and RMF
441  GArf arf_stacked = onoff->arf();
442  GRmf rmf_stacked = onoff->rmf();
443 
444  // Get exposure for observation
445  double exposure_on = onoff->on_spec().exposure();
446  double exposure_off = onoff->off_spec().exposure();
447 
448  // If this is the first On/Off observation then store the data to
449  // initialise the data definition
450  if (first) {
451 
452  // Store PHA, ARF and RMF, set exposure, ontime and livetime
453  m_on_spec = onoff->on_spec();
454  m_off_spec = onoff->off_spec();
455  m_arf = arf_stacked * exposure_on;
456  m_rmf = rmf_stacked;
457  total_exposure_on = exposure_on;
458  total_exposure_off = exposure_off;
459  m_ontime = onoff->ontime();
460  m_livetime = onoff->livetime();
461 
462  // Compute number of background events/MeV from BACKRESP column
463  // and store result intermediately into BACKRESP column of off
464  // spectum
465  std::vector<double>& backresp = m_off_spec["BACKRESP"];
466  for (int k = 0; k < backresp.size(); ++k) {
467  backresp[k] *= exposure_off;
468  }
469 
470  // Compute background scaling factor contribution
471  for (int k = 0; k < m_on_spec.size(); ++k) {
472  double background = onoff->off_spec()["BACKRESP"][k];
473  double alpha = onoff->on_spec().backscal(k);
474  double scale = alpha * background * exposure_off;
475  m_on_spec.backscal(k,scale);
476  }
477 
478  // Compute RMF contribution
479  for (int itrue = 0; itrue < m_rmf.ntrue(); ++itrue) {
480  double arf = m_arf[itrue];
481  for (int ireco = 0; ireco < m_rmf.nmeasured(); ++ireco) {
482  m_rmf(itrue,ireco) *= arf;
483  }
484  }
485 
486  // Set observation energy range
487  emin_obs = emin;
488  emax_obs = emax;
489 
490  // Signal that the On/Off definition has been set
491  first = false;
492 
493  }
494 
495  // ... otherwise stack data
496  else {
497 
498  // Check consistency of On spectrum
499  if (m_on_spec.ebounds() != onoff->on_spec().ebounds()) {
500  std::string msg = "Incompatible energy binning of On spectrum.";
502  }
503 
504  // Check consistency of Off spectrum
505  if (m_off_spec.ebounds() != onoff->off_spec().ebounds()) {
506  std::string msg = "Incompatible energy binning of Off spectrum.";
508  }
509 
510  // Check consistency of ARF
511  if (m_arf.ebounds() != onoff->arf().ebounds()) {
512  std::string msg = "Incompatible energy binning of ARF.";
514  }
515 
516  // Check consistency of RMF
517  if (m_rmf.etrue() != onoff->rmf().etrue()) {
518  std::string msg = "Incompatible true energy binning of RMF.";
520  }
521  if (m_rmf.emeasured() != onoff->rmf().emeasured()) {
522  std::string msg = "Incompatible measured energy binning of RMF.";
524  }
525 
526  // Add On and Off spectrum
527  m_on_spec += onoff->on_spec();
528  m_off_spec += onoff->off_spec();
529 
530  // Compute background scaling factor contribution
531  for (int k = 0; k < m_on_spec.size(); ++k) {
532  double background = onoff->off_spec()["BACKRESP"][k];
533  double alpha = onoff->on_spec().backscal(k);
534  double scale = alpha * background * exposure_off;
535  m_on_spec.backscal(k, m_on_spec.backscal(k) + scale);
536  }
537 
538  // Add ARF
539  m_arf += arf_stacked * exposure_on;
540 
541  // Add number of background events/MeV from BACKRESP column
542  const std::vector<double>& src = onoff->off_spec()["BACKRESP"];
543  std::vector<double>& dst = m_off_spec["BACKRESP"];
544  for (int k = 0; k < dst.size(); ++k) {
545  dst[k] += src[k] * exposure_off;
546  }
547 
548  // Add RMF
549  for (int itrue = 0; itrue < m_rmf.ntrue(); ++itrue) {
550  double arf = arf_stacked[itrue] * exposure_on;
551  for (int ireco = 0; ireco < m_rmf.nmeasured(); ++ireco) {
552  m_rmf(itrue,ireco) += rmf_stacked(itrue,ireco) * arf;
553  }
554  }
555 
556  // Add exposure, ontime and livetime
557  total_exposure_on += exposure_on;
558  total_exposure_off += exposure_off;
559  m_ontime += onoff->ontime();
560  m_livetime += onoff->livetime();
561 
562  // Update observation energy range
563  if (emin < emin_obs) {
564  emin_obs = emin;
565  }
566  if (emax > emax_obs) {
567  emax_obs = emax;
568  }
569 
570  } // endelse: stacked data
571 
572  } // endfor: looped over observations
573 
574  // Compute stacked background scaling factor
575  for (int k = 0; k < m_on_spec.size(); ++k) {
576  // at this point norm is already multiplied by exposure
577  // the ratio is performed later
578  double norm = m_off_spec["BACKRESP"][k];
579  if (norm > 0.0) {
580  double scale = m_on_spec.backscal(k) / norm;
581  m_on_spec.backscal(k,scale);
582  }
583  else {
584  m_on_spec.backscal(k,0.0);
585  }
586  }
587 
588  // Compute RMF
589  for (int itrue = 0; itrue < m_rmf.ntrue(); ++itrue) {
590  double arf = m_arf[itrue];
591  if (arf > 0.0) {
592  for (int ireco = 0; ireco < m_rmf.nmeasured(); ++ireco) {
593  m_rmf(itrue,ireco) /= arf;
594  }
595  }
596  }
597 
598  // Compute ARF
599  if (total_exposure_on > 0.0) {
600  m_arf /= total_exposure_on;
601  }
602 
603  // Compute background events/MeV/sec and store them in BACKRESP column
604  if (total_exposure_off > 0.0) {
605  std::vector<double>& backresp = m_off_spec["BACKRESP"];
606  for (int k = 0; k < backresp.size(); ++k) {
607  backresp[k] /= total_exposure_off;
608  }
609  }
610 
611  // Set energy boundaries of observation
612  m_on_spec.emin_obs(emin_obs);
613  m_on_spec.emax_obs(emax_obs);
614  m_off_spec.emin_obs(emin_obs);
615  m_off_spec.emax_obs(emax_obs);
616 
617  // Set instrument name
619 
620  // Return
621  return;
622 }
623 
624 
625 /***********************************************************************//**
626  * @brief Destructor
627  ***************************************************************************/
629 {
630  // Free members
631  free_members();
632 
633  // Return
634  return;
635 }
636 
637 
638 /*==========================================================================
639  = =
640  = Operators =
641  = =
642  ==========================================================================*/
643 
644 /***********************************************************************//**
645  * @brief Assignment operator
646  *
647  * @param[in] obs On/Off observation.
648  * @return On/Off observation.
649  *
650  * Assigns one On/Off observation to another On/Off observation object.
651  ***************************************************************************/
653 {
654  // Execute only if object is not identical
655  if (this != &obs) {
656 
657  // Copy base class members
658  this->GObservation::operator=(obs);
659 
660  // Free members
661  free_members();
662 
663  // Initialise private members for clean destruction
664  init_members();
665 
666  // Copy members
667  copy_members(obs);
668 
669  } // endif: object was not identical
670 
671  // Return
672  return *this;
673 }
674 
675 
676 /*==========================================================================
677  = =
678  = Public methods =
679  = =
680  ==========================================================================*/
681 
682 /***********************************************************************//**
683  * @brief Clear instance
684  *
685  * Clears the On/Off observation. All class members will be set to the
686  * initial state. Any information that was present in the object before will
687  * be lost.
688  ***************************************************************************/
690 {
691  // Free class members
692  free_members();
693 
694  // Initialise members
695  init_members();
696 
697  // Return
698  return;
699 }
700 
701 
702 /***********************************************************************//**
703  * @brief Clone instance
704  *
705  * @return Pointer to deep copy of On/Off observation.
706  *
707  * Returns a pointer to a deep copy of an On/Off observation.
708  **************************************************************************/
710 {
711  return new GCTAOnOffObservation(*this);
712 }
713 
714 
715 /***********************************************************************//**
716  * @brief Set response function
717  *
718  * @param[in] rsp Response function.
719  *
720  * @exception GException::invalid_argument
721  * Invalid response class specified.
722  *
723  * Sets the response function for the On/Off observation.
724  ***************************************************************************/
726 {
727  // Retrieve CTA response pointer
728  const GCTAResponse* ptr = gammalib::cta_rsp(G_RESPONSE_SET, rsp);
729 
730  // Free existing response only if it differs from current response. This
731  // prevents unintential deallocation of the argument
732  if ((m_response != NULL) && (m_response != &rsp)) {
733  delete m_response;
734  }
735 
736  // Clone response function
737  m_response = ptr->clone();
738 
739  // Return
740  return;
741 }
742 
743 
744 /***********************************************************************//**
745  * @brief Return pointer to CTA response function
746  *
747  * @return Pointer to CTA response function.
748  *
749  * @exception GException::invalid_value
750  * No valid response found in CTA observation.
751  *
752  * Returns a pointer to the CTA response function. An exception is thrown if
753  * the pointer is not valid, hence the user does not need to verify the
754  * validity of the pointer.
755  ***************************************************************************/
757 {
758  // Throw an exception if the response pointer is not valid
759  if (m_response == NULL) {
760  std::string msg = "No valid response function found in CTA On/Off "
761  "observation.\n";
763  }
764 
765  // Return pointer
766  return m_response;
767 }
768 
769 
770 /***********************************************************************//**
771  * @brief Read On/Off observation from an XML element
772  *
773  * @param[in] xml XML element.
774  *
775  * @exception GException::invalid_value
776  * Invalid statistic attribute encountered
777  *
778  * Reads information for an On/Off observation from an XML element. The
779  * expected format of the XML element is
780  *
781  * <observation name="..." id="..." instrument="...">
782  * <parameter name="Pha_on" file="..."/>
783  * <parameter name="Pha_off" file="..."/>
784  * <parameter name="Arf" file="..."/>
785  * <parameter name="Rmf" file="..."/>
786  * </observation>
787  *
788  * Optionally, the statistic used for maximum likelihood fitting can be
789  * specified:
790  *
791  * <observation name="..." id="..." instrument="..." statistic="...">
792  *
793  ***************************************************************************/
795 {
796  // clean object
797  clear();
798 
799  // Extract instrument name
800  m_instrument = xml.attribute("instrument");
801 
802  // Read in user defined statistic for this observation
803  if (xml.attribute("statistic") != "") {
804 
805  // Extract statistic value
806  std::string statistic = gammalib::toupper(xml.attribute("statistic"));
807 
808  // If statistic is not POISSON, CSTAT or WSTAT than throw an exception
809  if ((statistic != "POISSON") &&
810  (statistic != "CSTAT") &&
811  (statistic != "WSTAT")) {
812  std::string msg = "Invalid statistic \""+statistic+"\" encountered "
813  "in observation definition XML file for "
814  "\""+m_instrument+"\" observation with identifier "
815  "\""+xml.attribute("id")+"\". Only \"POISSON\" "
816  ", \"CSTAT\" or \"WSTAT\" are supported.";
817  throw GException::invalid_value(G_READ, msg);
818  }
819 
820  // Save statistic value
821  this->statistic(xml.attribute("statistic"));
822 
823  }
824 
825  // Get file names
826  std::string pha_on = gammalib::xml_get_attr(G_READ, xml, "Pha_on", "file");
827  std::string pha_off = gammalib::xml_get_attr(G_READ, xml, "Pha_off", "file");
828  std::string arf = gammalib::xml_get_attr(G_READ, xml, "Arf", "file");
829  std::string rmf = gammalib::xml_get_attr(G_READ, xml, "Rmf", "file");
830 
831  // Expand file names
832  pha_on = gammalib::xml_file_expand(xml, pha_on);
833  pha_off = gammalib::xml_file_expand(xml, pha_off);
834  arf = gammalib::xml_file_expand(xml, arf);
835  rmf = gammalib::xml_file_expand(xml, rmf);
836 
837  // Load files
838  m_on_spec.load(pha_on);
839  m_off_spec.load(pha_off);
840  m_arf.load(arf);
841  m_rmf.load(rmf);
842 
843  // Set the ontime, livetime and deadtime correction
844  set_exposure();
845 
846  // Check consistency of On/Off observation
848 
849  // Return
850  return;
851 }
852 
853 
854 /***********************************************************************//**
855  * @brief write observation to an xml element
856  *
857  * @param[in] xml XML element.
858  *
859  * Writes information for an On/Off observation into an XML element. The
860  * expected format of the XML element is
861  *
862  * <observation name="..." id="..." instrument="..." statistic="...">
863  * <parameter name="Pha_on" file="..."/>
864  * <parameter name="Pha_off" file="..."/>
865  * <parameter name="Arf" file="..."/>
866  * <parameter name="Rmf" file="..."/>
867  * </observation>
868  *
869  * The actual files described in the XML elements are not written.
870  ***************************************************************************/
872 {
873  // Allocate XML element pointer
874  GXmlElement* par;
875 
876  // Set Pha_on parameter
877  par = gammalib::xml_need_par(G_WRITE, xml, "Pha_on");
879 
880  // Set Pha_off parameter
881  par = gammalib::xml_need_par(G_WRITE, xml, "Pha_off");
883 
884  // Set Arf parameter
885  par = gammalib::xml_need_par(G_WRITE, xml, "Arf");
886  par->attribute("file", gammalib::xml_file_reduce(xml, m_arf.filename()));
887 
888  // Set Rmf parameter
889  par = gammalib::xml_need_par(G_WRITE, xml, "Rmf");
890  par->attribute("file", gammalib::xml_file_reduce(xml, m_rmf.filename()));
891 
892  // Add user defined statistic attributes
893  if (statistic() != "") {
894  xml.attribute("statistic", statistic());
895  }
896 
897  // Return
898  return;
899 }
900 
901 
902 /***********************************************************************//**
903  * @brief Evaluate log-likelihood function for On/Off analysis
904  *
905  * @param[in] models Models.
906  * @param[in,out] gradient Pointer to gradients.
907  * @param[in,out] curvature Pointer to curvature matrix.
908  * @param[in,out] npred Pointer to Npred value.
909  * @return Log-likelihood value.
910  *
911  * @exception GException::invalid_value
912  * Invalid statistic encountered.
913  ***************************************************************************/
915  GVector* gradient,
916  GMatrixSparse* curvature,
917  double* npred) const
918 {
919  // Initialise likelihood value
920  double value = 0.0;
921 
922  // Extract statistic for this observation
923  std::string statistic = gammalib::toupper(this->statistic());
924 
925  // Poisson statistic with modeled background
926  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
927  value = likelihood_cstat(models, gradient, curvature, npred);
928  }
929 
930  // ... or Poisson statistic with measured background
931  else if (statistic == "WSTAT") {
932  value = likelihood_wstat(models, gradient, curvature, npred);
933  }
934 
935  // ... or unsupported
936  else {
937  std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
938  "Either specify \"POISSON\", \"CSTAT\" or "
939  "\"WSTAT\".";
941  }
942 
943  // Return likelihood
944  return value;
945 }
946 
947 
948 /***********************************************************************
949  * @brief Compute predicted gamma-ray counts for given model
950  *
951  * @param[in] models Model container.
952  * @return Model PHA with predicted gamma-ray counts.
953  *
954  * Returns spectrum of predicted number of source events in the On
955  * regions.
956  ***********************************************************************/
958 {
959  // Get number of model parameters in model container
960  int npars = models.npars();
961 
962  // Create dummy gradient vectors to provide required argument to N_gamma
963  GVector dummy_grad(npars);
964 
965  // Initialise output spectrum
966  GEbounds ereco = m_on_spec.ebounds();
967  GPha gammas = GPha(ereco);
968 
969  // Loop over all energy bins
970  for (int i = 0; i < m_on_spec.size(); ++i) {
971  gammas.fill(ereco.emean(i), N_gamma(models, i, &dummy_grad));
972  }
973 
974  // Return predicted gamma-ray counts spectrum
975  return gammas;
976 }
977 
978 
979 /***********************************************************************
980  * @brief Compute predicted background counts PHA for given model
981  *
982  * @param[in] models Model container.
983  * @return Model background PHA.
984  *
985  * Returns spectrum of predicted number of background events in the Off
986  * regions. The computation method changed depending on the statistic used
987  * to match the spectrum used for likelihood fitting.
988  * To be noted: for the WSTAT statistic the model background depends on
989  * on the model adopted for gamma rays, that therefore need to be passed
990  * to this function.
991  ***********************************************************************/
993 {
994  // Get number of model parameters in model container
995  int npars = models.npars();
996 
997  // Create dummy gradient vectors to provide required argument to N_gamma
998  GVector dummy_grad(npars);
999 
1000  // Initialise output spectrum
1001  GEbounds ereco = m_on_spec.ebounds();
1002  GPha bgds = GPha(ereco);
1003 
1004  // Extract statistic for this observation
1005  std::string statistic = gammalib::toupper(this->statistic());
1006 
1007  // Loop over all energy bins
1008  for (int i = 0; i < m_on_spec.size(); ++i) {
1009 
1010  // Initialise variable to store number of background counts
1011  double nbgd = 0.0;
1012 
1013  // Choose background evaluation method based on fit statistics
1014 
1015  // CSTAT
1016  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
1017  nbgd = N_bgd(models, i, &dummy_grad);
1018  }
1019 
1020  // ... or Poisson statistic with measured background
1021  else if (statistic == "WSTAT") {
1022 
1023  // Get number of On and Off counts
1024  double non = m_on_spec[i];
1025  double noff = m_off_spec[i];
1026 
1027  // Get background scaling
1028  double alpha = m_on_spec.backscal(i);
1029 
1030  // Get number of gamma and background events (and corresponding
1031  // spectral model gradients)
1032  double ngam = N_gamma(models, i, &dummy_grad);
1033 
1034  // Initialise variables for likelihood computation
1035  double nonpred = 0.0;
1036  double dlogLdsky = 0.0;
1037  double d2logLdsky2 = 0.0;
1038 
1039  // Perform likelihood profiling and derive likelihood value
1040  // The number of background counts is updated to the profile value
1041  wstat_value(non, noff, alpha, ngam, nonpred, nbgd,
1042  dlogLdsky,d2logLdsky2);
1043 
1044  } // endelse: WSTAT
1045 
1046  // ... or unsupported
1047  else {
1048  std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
1049  "Either specify \"POISSON\", \"CSTAT\" or "
1050  "\"WSTAT\".";
1052  }
1053 
1054  // Fill background spectrum
1055  bgds.fill(ereco.emean(i), nbgd);
1056 
1057  } // endfor: looped over energy bins
1058 
1059  // Return model count spectrum
1060  return bgds;
1061 }
1062 
1063 
1064 /***********************************************************************//**
1065  * @brief Print On/Off observation information
1066  *
1067  * @param[in] chatter Chattiness.
1068  * @return String containing On/Off observation information.
1069  ***************************************************************************/
1070 std::string GCTAOnOffObservation::print(const GChatter& chatter) const
1071 {
1072  // Initialise result string
1073  std::string result;
1074 
1075  // Continue only if chatter is not silent
1076  if (chatter != SILENT) {
1077 
1078  // Append header
1079  result.append("=== GCTAOnOffObservation ===");
1080 
1081  // Append parameters
1082  result.append("\n"+gammalib::parformat("Name")+m_name);
1083  result.append("\n"+gammalib::parformat("Identifier")+m_id);
1084  result.append("\n"+gammalib::parformat("Instrument")+instrument());
1085  result.append("\n"+gammalib::parformat("Statistic")+statistic());
1086  result.append("\n"+gammalib::parformat("Ontime"));
1087  result.append(gammalib::str(ontime())+" s");
1088  result.append("\n"+gammalib::parformat("Livetime"));
1089  result.append(gammalib::str(livetime())+" s");
1090  result.append("\n"+gammalib::parformat("Deadtime correction"));
1091  result.append(gammalib::str(m_deadc));
1092 
1093  // Append spectra, ARF and RMF
1094  result.append("\n"+m_on_spec.print(gammalib::reduce(chatter)));
1095  result.append("\n"+m_off_spec.print(gammalib::reduce(chatter)));
1096  result.append("\n"+m_arf.print(gammalib::reduce(chatter)));
1097  result.append("\n"+m_rmf.print(gammalib::reduce(chatter)));
1098  }
1099 
1100  // Return result
1101  return result;
1102 }
1103 
1104 
1105 /*==========================================================================
1106  = =
1107  = Private methods =
1108  = =
1109  ==========================================================================*/
1110 
1111 /***********************************************************************//**
1112  * @brief Initialise class members
1113  ***************************************************************************/
1115 {
1116  // Initialise members
1117  m_instrument = "CTAOnOff";
1118  m_response = NULL;
1119  m_ontime = 0.0;
1120  m_livetime = 0.0;
1121  m_deadc = 1.0;
1122  m_on_spec.clear();
1123  m_off_spec.clear();
1124  m_arf.clear();
1125  m_rmf.clear();
1126 
1127  // Return
1128  return;
1129 }
1130 
1131 
1132 /***********************************************************************//**
1133  * @brief Copy class members
1134  *
1135  * @param[in] obs CTA On/Off observation.
1136  ***************************************************************************/
1138 {
1139  // Copy attributes
1140  m_instrument = obs.m_instrument;
1141  m_ontime = obs.m_ontime;
1142  m_livetime = obs.m_livetime;
1143  m_deadc = obs.m_deadc;
1144  m_on_spec = obs.m_on_spec;
1145  m_off_spec = obs.m_off_spec;
1146  m_arf = obs.m_arf;
1147  m_rmf = obs.m_rmf;
1148 
1149  // Clone members
1150  m_response = (obs.m_response != NULL) ? obs.m_response->clone() : NULL;
1151 
1152  // Return
1153  return;
1154 }
1155 
1156 
1157 
1158 /***********************************************************************//**
1159  * @brief Delete class members
1160  ***************************************************************************/
1162 {
1163  // Return
1164  return;
1165 }
1166 
1167 
1168 /***********************************************************************//**
1169  * @brief Set ontime, livetime and deadtime correction factor
1170  ***************************************************************************/
1172 {
1173  // Set the ontime, livetime and deadtime correction
1176  m_deadc = 1.0;
1177 
1178  // Return
1179  return;
1180 }
1181 
1182 
1183 /***********************************************************************//**
1184  * @brief Check consistency of data members
1185  *
1186  * @param[in] method Calling method.
1187  *
1188  * @exception GException::invalid_value
1189  * Inconsistency found.
1190  *
1191  * Checks the consistency of data members and throw exceptions in case that
1192  * inconsistencies are found.
1193  ***************************************************************************/
1194 void GCTAOnOffObservation::check_consistency(const std::string& method) const
1195 {
1196  // Check that On spectrum is consistent with Off spectrum
1197  if (m_on_spec.ebounds() != m_off_spec.ebounds()) {
1198  std::string msg = "On and Off spectra are incompatible.";
1199  throw GException::invalid_value(method, msg);
1200  }
1201 
1202  // Check that On spectrum is consistent with RMF
1203  if (m_on_spec.ebounds() != m_rmf.emeasured()) {
1204  std::string msg = "Redistribution Matrix File is incompatible with "
1205  "spectra.";
1206  throw GException::invalid_value(method, msg);
1207  }
1208 
1209  // Check that ARF is consistent with RMF
1210  if (m_arf.ebounds() != m_rmf.etrue()) {
1211  std::string msg = "Redistribution Matrix File is incompatible with "
1212  "Auxiliary Response File.";
1213  throw GException::invalid_value(method, msg);
1214  }
1215 
1216  // Return
1217  return;
1218 }
1219 
1220 
1221 /***********************************************************************//**
1222  * @brief Return ARF for stacking
1223  *
1224  * @param[in] arf Auxiliary Response File.
1225  * @param[in] emin Minimum observation energy.
1226  * @param[in] emax Maximum observation energy.
1227  * @return Auxiliary Response File for stacking.
1228  *
1229  * Returns an Auxiliary Response File for stacking that has all elements with
1230  * true energies outside the interval [@p emin, @p emax] clipped to zero.
1231  * This prevents spill over from one observation into another.
1232  ***************************************************************************/
1234  const GEnergy& emin,
1235  const GEnergy& emax) const
1236 {
1237  // Copy ARF
1238  GArf arf_stacked = arf;
1239 
1240  // Get true energy boundaries of ARF
1241  const GEbounds& etrue = arf_stacked.ebounds();
1242 
1243  // Set all ARF bins outside energy interval to zero
1244  for (int itrue = 0; itrue < etrue.size(); ++itrue) {
1245  if ((etrue.emax(itrue) < emin) || (etrue.emin(itrue) > emax)) {
1246  arf_stacked[itrue] = 0.0;
1247  }
1248  }
1249 
1250  // Return ARF
1251  return arf_stacked;
1252 }
1253 
1254 
1255 /***********************************************************************//**
1256  * @brief Return RMF for stacking
1257  *
1258  * @param[in] rmf Redistribution Matrix File.
1259  * @param[in] emin Minimum observation energy.
1260  * @param[in] emax Maximum observation energy.
1261  * @return Redistribution Matrix File for stacking.
1262  *
1263  * Returns a Redistribution Matrix File for stacking that has all matrix
1264  * elements with true energies outside the interval [@p emin, @p emax]
1265  * clipped to zero. This prevents spill over from one observation into
1266  * another.
1267  *
1268  * To correct for the missing events at the edges, the Redistribution Matrix
1269  * File is renormalized so that for each reconstructed energy the sum over
1270  * the true energies is the same as before the clipping.
1271  ***************************************************************************/
1273  const GEnergy& emin,
1274  const GEnergy& emax) const
1275 {
1276  // Copy RMF
1277  GRmf rmf_stacked = rmf;
1278 
1279  // Get true energy boundaries of RMF
1280  const GEbounds& etrue = rmf_stacked.etrue();
1281  int nreco = rmf_stacked.emeasured().size();
1282 
1283  // Get RMF totals for all reconstructed energy bins
1284  std::vector<double> sums(nreco, 0.0);
1285  for (int ireco = 0; ireco < nreco; ++ireco) {
1286  for (int itrue = 0; itrue < etrue.size(); ++itrue) {
1287  sums[ireco] += rmf_stacked(itrue, ireco);
1288  }
1289  }
1290 
1291  // Set all RMF bins outside energy interval to zero (clipping)
1292  for (int itrue = 0; itrue < etrue.size(); ++itrue) {
1293  if ((etrue.emax(itrue) < emin) || (etrue.emin(itrue) > emax)) {
1294  for (int ireco = 0; ireco < nreco; ++ireco) {
1295  rmf_stacked(itrue, ireco) = 0.0;
1296  }
1297  }
1298  }
1299 
1300  // Renormalize RMF
1301  for (int ireco = 0; ireco < nreco; ++ireco) {
1302  double sum = 0.0;
1303  for (int itrue = 0; itrue < etrue.size(); ++itrue) {
1304  sum += rmf_stacked(itrue, ireco);
1305  }
1306  if (sum > 0.0) {
1307  double norm = sums[ireco] / sum;
1308  for (int itrue = 0; itrue < etrue.size(); ++itrue) {
1309  rmf_stacked(itrue, ireco) *= norm;
1310  }
1311  }
1312  }
1313 
1314  // Return RMF
1315  return rmf_stacked;
1316 }
1317 
1318 
1319 /***********************************************************************//**
1320  * @brief Set On/Off observation from a CTA observation
1321  *
1322  * @param[in] obs_on CTA On observation.
1323  * @param[in] obs_off CTA Off observation.
1324  * @param[in] models Models including source and background model.
1325  * @param[in] srcname Name of soucre in models.
1326  * @param[in] on On regions.
1327  * @param[in] off Off regions.
1328  * @param[in] use_model_bkg Use model background.
1329  *
1330  * @exception GException::invalid_value
1331  * No CTA event list found in CTA observation.
1332  *
1333  * Sets an On/Off observation from a CTA observation by filling the events
1334  * that fall in the On and Off regions into the PHA spectra and by computing
1335  * the corresponding ARF and RMF response functions.
1336  *
1337  * The @p use_model_bkg flags controls whether the background model should
1338  * be used for computations or not. This impacts the computations of the
1339  * `BACKSCAL` column in the On spectrum and the `BACKRESP` column in the Off
1340  * spectrum. See the compute_alpha() and compute_bgd() methods for more
1341  * information on the applied formulae.
1342  ***************************************************************************/
1344  const GCTAObservation& obs_off,
1345  const GModels& models,
1346  const std::string& srcname,
1347  const GSkyRegions& on,
1348  const GSkyRegions& off,
1349  const bool& use_model_bkg)
1350 {
1351  // Get CTA event list pointers
1352  const GCTAEventList* events_on = dynamic_cast<const GCTAEventList*>(obs_on.events());
1353  if (events_on == NULL) {
1354  std::string msg = "No event list found in CTA On observation \""+
1355  obs_on.name()+"\" (ID="+obs_on.id()+"). ON/OFF observation "
1356  "can only be filled from event list.";
1357  throw GException::invalid_value(G_SET, msg);
1358  }
1359  const GCTAEventList* events_off = dynamic_cast<const GCTAEventList*>(obs_off.events());
1360  if (events_off == NULL) {
1361  std::string msg = "No event list found in CTA Off observation \""+
1362  obs_off.name()+"\" (ID="+obs_off.id()+"). ON/OFF observation "
1363  "can only be filled from event list.";
1364  throw GException::invalid_value(G_SET, msg);
1365  }
1366 
1367  // Get spatial component of source model
1368  const GModelSky* model = dynamic_cast<const GModelSky*>(models[srcname]);
1369  if (model == NULL) {
1370  std::string msg = "Model \""+srcname+"\" is not a sky model. Please "
1371  "specify the name of a sky model.";
1372  throw GException::invalid_value(G_SET, msg);
1373  }
1374  const GModelSpatial& spatial = *(model->spatial());
1375 
1376  // If background model is needed then extract background model components
1377  // from model container
1378  GModels bkg_models;
1379  if (use_model_bkg) {
1380  for (int i = 0; i < models.size(); ++i) {
1381  if ((dynamic_cast<const GModelData*>(models[i]) != NULL) &&
1382  (models[i]->classname().substr(0,4) == "GCTA")) {
1383  bkg_models.append(*models[i]);
1384  }
1385  }
1386  }
1387 
1388  // Loop over all events and fill On spectrum
1389  for (int i = 0; i < events_on->size(); ++i) {
1390 
1391  // Get measured event direction
1392  const GCTAEventAtom* atom = (*events_on)[i];
1393  GSkyDir dir = atom->dir().dir();
1394 
1395  // Fill in spectrum according to region containment
1396  if (on.contains(dir)) {
1397  m_on_spec.fill(atom->energy());
1398  }
1399 
1400  } // endfor: looped over all On events
1401 
1402  // Loop over all events and fill Off spectrum
1403  for (int i = 0; i < events_off->size(); ++i) {
1404 
1405  // Get measured event direction
1406  const GCTAEventAtom* atom = (*events_off)[i];
1407  GSkyDir dir = atom->dir().dir();
1408 
1409  // Fill in spectrum according to region containment
1410  if (off.contains(dir)) {
1411  m_off_spec.fill(atom->energy());
1412  }
1413 
1414  } // endfor: looped over all Off events
1415 
1416  // Store the livetime as exposures of the spectra
1417  m_on_spec.exposure(obs_on.livetime());
1418  m_off_spec.exposure(obs_off.livetime());
1419 
1420  // Store the ontime, livetime and deadtime correction in the On observation
1421  // These variables are used so far only for display
1422  m_ontime = obs_on.ontime();
1423  m_livetime = obs_on.livetime();
1424  m_deadc = obs_on.deadc();
1425 
1426  // Convert all regions into region maps
1427  GSkyRegions reg_on;
1428  GSkyRegions reg_off;
1429  for (int i = 0; i < on.size(); ++i) {
1430  reg_on.append(GSkyRegionMap(on[i]));
1431  }
1432  for (int i = 0; i < off.size(); ++i) {
1433  reg_off.append(GSkyRegionMap(off[i]));
1434  }
1435 
1436  // Get effective area radius cut
1437  double rad_max = arf_rad_max(obs_on, on);
1438 
1439  // Compute response components
1440  if (rad_max > 0.0) {
1441  compute_arf_cut(obs_on, spatial, reg_on);
1442  }
1443  else {
1444  compute_arf(obs_on, spatial, reg_on);
1445  }
1446  compute_bgd(obs_off, reg_off, bkg_models, use_model_bkg);
1447  compute_alpha(obs_on, obs_off, reg_on, reg_off, bkg_models, use_model_bkg);
1448  compute_rmf(obs_on, reg_on);
1449 
1450  // Apply reconstructed energy boundaries
1451  apply_ebounds(obs_on);
1452  apply_ebounds(obs_off);
1453 
1454  // Set observation energy band
1455  m_on_spec.emin_obs(obs_on.ebounds().emin());
1456  m_on_spec.emax_obs(obs_on.ebounds().emax());
1457  m_off_spec.emin_obs(obs_off.ebounds().emin());
1458  m_off_spec.emax_obs(obs_off.ebounds().emax());
1459 
1460  // Get mission, instrument and response name
1461  std::string mission;
1462  std::string instrument;
1463  std::string rspname;
1464  const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs_on.response());
1465  if (rsp != NULL) {
1466  mission = gammalib::toupper(rsp->caldb().mission());
1467  instrument = gammalib::toupper(rsp->caldb().instrument());
1468  rspname = gammalib::toupper(rsp->rspname());
1469  }
1470 
1471  // Set header keywords
1472  GFitsHeader arf_header;
1473  GFitsHeader pha_header;
1474  GFitsHeader rmf_header;
1475  if (!mission.empty()) {
1476  arf_header.append(GFitsHeaderCard("TELESCOP", mission, "Telescope"));
1477  pha_header.append(GFitsHeaderCard("TELESCOP", mission, "Telescope"));
1478  rmf_header.append(GFitsHeaderCard("TELESCOP", mission, "Telescope"));
1479  }
1480  if (!instrument.empty()) {
1481  arf_header.append(GFitsHeaderCard("INSTRUME", instrument, "Instrument"));
1482  pha_header.append(GFitsHeaderCard("INSTRUME", instrument, "Instrument"));
1483  rmf_header.append(GFitsHeaderCard("INSTRUME", instrument, "Instrument"));
1484  }
1485  if (!rspname.empty()) {
1486  arf_header.append(GFitsHeaderCard("RESPNAME", rspname, "Response name"));
1487  pha_header.append(GFitsHeaderCard("RESPNAME", rspname, "Response name"));
1488  rmf_header.append(GFitsHeaderCard("RESPNAME", rspname, "Response name"));
1489  }
1490  if (rad_max > 0.0) {
1491  arf_header.append(GFitsHeaderCard("RAD_MAX", rad_max, "[deg] Applied radius cut"));
1492  }
1493 
1494  // Store header keywords
1495  m_on_spec.header(pha_header);
1496  m_off_spec.header(pha_header);
1497  m_arf.header(arf_header);
1498  m_rmf.header(rmf_header);
1499 
1500  // Set instrument name
1501  m_instrument = obs_on.instrument() + "OnOff";
1502 
1503  // Return
1504  return;
1505 }
1506 
1507 
1508 /***********************************************************************//**
1509  * @brief Compute ARF of On/Off observation
1510  *
1511  * @param[in] obs CTA observation.
1512  * @param[in] spatial Spatial source model.
1513  * @param[in] on On regions.
1514  *
1515  * @exception GException::invalid_argument
1516  * No CTA response found in CTA observation.
1517  *
1518  * Computes the ARF for an On/Off observation by integration over the IRF
1519  * for the specified @p spatial source model over the @p on regions.
1520  *
1521  * All On regions contained in @p on are expected to be sky region maps,
1522  * i.e. of type GSkyRegionMap.
1523  ***************************************************************************/
1525  const GModelSpatial& spatial,
1526  const GSkyRegions& on)
1527 {
1528  // Get reconstructed energy boundaries from on ARF
1529  const GEbounds& etrue = m_arf.ebounds();
1530  int ntrue = etrue.size();
1531 
1532  // Continue only if there are ARF bins
1533  if (ntrue > 0) {
1534 
1535  // Get CTA IRF response
1537 
1538  // Set dummy time
1539  const GTime time;
1540 
1541  // Save original energy dispersion application status
1542  bool save_edisp = rsp.apply_edisp();
1543 
1544  // Switch-off application of energy dispersion
1545  const_cast<GCTAResponseIrf&>(rsp).apply_edisp(false);
1546 
1547  // Loop over true energies
1548  for (int i = 0; i < ntrue; ++i) {
1549 
1550  // Get mean energy of bin
1551  GEnergy energy = etrue.elogmean(i);
1552 
1553  // Set source
1554  GSource source("", const_cast<GModelSpatial*>(&spatial), energy, time);
1555 
1556  // Initialize effective area for this bin
1557  m_arf[i] = 0.0;
1558 
1559  // Loop over regions
1560  for (int k = 0; k < on.size(); ++k) {
1561 
1562  // Get pointer on sky region map
1563  const GSkyRegionMap* on_map = static_cast<const GSkyRegionMap*>(on[k]);
1564 
1565  // Loop over pixels in On region map and integrate effective
1566  // area
1567  for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {
1568 
1569  // Get pixel index
1570  int pixidx = on_map->nonzero_indices()[j];
1571 
1572  // Get direction to pixel center
1573  GCTAInstDir pixdir = GCTAInstDir(on_map->map().inx2dir(pixidx));
1574 
1575  // Get solid angle subtended by this pixel
1576  double pixsolid = on_map->map().solidangle(pixidx);
1577 
1578  // Set event
1579  GCTAEventBin event;
1580  event.dir(pixdir);
1581  event.energy(energy);
1582  event.time(time);
1583 
1584  // Get ARF value. We need to devide by the deadtime
1585  // correction since the IRF method multiplies with it.
1586  double irf = rsp.irf_spatial(event, source, obs) / m_deadc;
1587 
1588  // Add up effective area
1589  m_arf[i] += irf * pixsolid;
1590 
1591  } // endfor: looped over all pixels in region map
1592 
1593  } // endfor: looped over all regions
1594 
1595  } // endfor: looped over true energies
1596 
1597  // Put back original energy dispersion application status
1598  const_cast<GCTAResponseIrf&>(rsp).apply_edisp(save_edisp);
1599 
1600  } // endif: there were energy bins
1601 
1602  // Return
1603  return;
1604 }
1605 
1606 
1607 /***********************************************************************//**
1608  * @brief Compute ARF of On/Off observation for a IRF with radius cut
1609  *
1610  * @param[in] obs CTA observation.
1611  * @param[in] spatial Spatial source model.
1612  * @param[in] on On regions.
1613  *
1614  * @exception GException::invalid_argument
1615  * No CTA response found in CTA observation.
1616  *
1617  * Computes the ARF for an On/Off observation by computing the average
1618  * effective area over the @p on regions for the specified @p spatial source
1619  * model.
1620  *
1621  * All On regions contained in @p on are expected to be sky region maps,
1622  * i.e. of type GSkyRegionMap.
1623  ***************************************************************************/
1625  const GModelSpatial& spatial,
1626  const GSkyRegions& on)
1627 {
1628  // Get reconstructed energy boundaries from on ARF
1629  const GEbounds& etrue = m_arf.ebounds();
1630  int ntrue = etrue.size();
1631 
1632  // Continue only if there are ARF bins
1633  if (ntrue > 0) {
1634 
1635  // Get CTA IRF response
1637 
1638  // Get CTA observation pointing direction, zenith, and azimuth
1639  GCTAPointing obspnt = obs.pointing();
1640  GSkyDir obsdir = obspnt.dir();
1641  double zenith = obspnt.zenith();
1642  double azimuth = obspnt.azimuth();
1643 
1644  // Loop over true energies
1645  for (int i = 0; i < ntrue; ++i) {
1646 
1647  // Get mean energy of bin
1648  double logEtrue = etrue.elogmean(i).log10TeV();
1649 
1650  // Initialize effective area and solid angle sum for this bin
1651  m_arf[i] = 0.0;
1652  double sum = 0.0;
1653 
1654  // Loop over regions
1655  for (int k = 0; k < on.size(); ++k) {
1656 
1657  // Get pointer on sky region map
1658  const GSkyRegionMap* on_map = static_cast<const GSkyRegionMap*>(on[k]);
1659 
1660  // Loop over pixels in On region map and integrate effective
1661  // area
1662  for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {
1663 
1664  // Get pixel index
1665  int pixidx = on_map->nonzero_indices()[j];
1666 
1667  // Get direction to pixel center
1668  GSkyDir pixdir = on_map->map().inx2dir(pixidx);
1669 
1670  // Get solid angle subtended by this pixel
1671  double pixsolid = on_map->map().solidangle(pixidx);
1672 
1673  // Compute position of pixel centre in instrument coordinates
1674  double theta = obsdir.dist(pixdir);
1675  double phi = obsdir.posang(pixdir); // Celestial system
1676 
1677  // Add up effective area
1678  m_arf[i] += rsp.aeff(theta, phi,
1679  zenith, azimuth,
1680  logEtrue) * pixsolid;
1681 
1682  // Sum up solid angles
1683  sum += pixsolid;
1684 
1685  } // endfor: looped over all pixels in region map
1686 
1687  } // endfor: looped over all regions
1688 
1689  // Divide effective area by solid angle so that ARF contains the
1690  // average effective area over the On region
1691  if (sum > 0.0) {
1692  m_arf[i] /= sum;
1693  }
1694  else {
1695  m_arf[i] = 0.0;
1696  }
1697 
1698  } // endfor: looped over true energies
1699 
1700  } // endif: there were energy bins
1701 
1702  // Return
1703  return;
1704 }
1705 
1706 
1707 /***********************************************************************//**
1708  * @brief Compute background rate in Off regions
1709  *
1710  * @param[in] obs CTA observation.
1711  * @param[in] off Off regions.
1712  * @param[in] models Models including background model.
1713  * @param[in] use_model_bkg Use model background.
1714  *
1715  * Computes the background rate in units of
1716  * \f${\rm events} \, {\rm MeV}^{-1} \, {\rm s}^{-1}\f$
1717  * in the Off regions and stores the result as additional column with name
1718  * `BACKRESP` in the Off spectrum.
1719  *
1720  * All Off regions contained in @p off are expected to be sky region maps,
1721  * i.e. of type GSkyRegionMap.
1722  *
1723  * If @p use_model_bkg is @c true, the IRF background template will be used
1724  * for the computation, and `BACKRESP` will be computed using
1725  *
1726  * \f[
1727  * {\tt BACKRESP}(E_{\rm reco}) = \sum_{\rm off} \sum_p
1728  * {\tt BKG}_{{\rm off},p}(E_{\rm reco})
1729  * \times \Omega_{{\rm off},p}
1730  * \f]
1731  *
1732  * where \f${\rm off}\f$ is the index of the Off region and \f$p\f$ is the
1733  * pixel in the Off region (note that each Off region is transformed into a
1734  * region map and \f$p\f$ indicates the pixels of this map).
1735  * \f${\tt BKG}_{{\rm off},p}(E_{\rm reco})\f$ is the background rate as
1736  * provided by the IRF background template in units of
1737  * \f${\rm events} \, {\rm MeV}^{-1} \, {\rm s}^{-1} \, {\rm sr}^{-1}\f$
1738  * for a reconstructed energy \f$E_{\rm reco}\f$ and a pixel index \f$p\f$
1739  * in the Off region \f${\rm off}\f$.
1740  * \f$\Omega_{{\rm off},p}\f$ is the solid angle in units of \f${\rm sr}\f$
1741  * of the pixel index \f$p\f$ in the Off region \f${\rm off}\f$.
1742  *
1743  * If @p use_model_bkg is @c false, `BACKRESP` will be computed using
1744  *
1745  * \f[
1746  * {\tt BACKRESP}(E_{\rm reco}) = \sum_{\rm off} \sum_p \Omega_{{\rm off},p}
1747  * \f]
1748  *
1749  * which actually assumes that the background rate is constant over the
1750  * field of view.
1751  *
1752  * @todo Integrate background rate over energy bin instead of computing the
1753  * rate at the energy bin centre.
1754  ***************************************************************************/
1756  const GSkyRegions& off,
1757  const GModels& models,
1758  const bool& use_model_bkg)
1759 {
1760  // Get reconstructed energies for the background rates
1761  const GEbounds& ereco = m_off_spec.ebounds();
1762  int nreco = ereco.size();
1763 
1764  // Continue only if there are energy bins
1765  if (nreco > 0) {
1766 
1767  // Initialise background rates to zero
1768  std::vector<double> background(nreco, 0.0);
1769 
1770  // Get CTA observation pointing direction
1771  GCTAPointing obspnt = obs.pointing();
1772 
1773  // Continue only if the IRF background template should be used
1774  if (use_model_bkg) {
1775 
1776  // Loop over regions
1777  for (int k = 0; k < off.size(); ++k) {
1778 
1779  // Get pointer on sky region map
1780  const GSkyRegionMap* off_map =
1781  static_cast<const GSkyRegionMap*>(off[k]);
1782 
1783  // Loop over pixels in Off region map and integrate
1784  // background rate
1785  for (int j = 0; j < off_map->nonzero_indices().size(); ++j) {
1786 
1787  // Get pixel index
1788  int pixidx = off_map->nonzero_indices()[j];
1789 
1790  // Get direction to pixel center
1791  GSkyDir pixdir = off_map->map().inx2dir(pixidx);
1792 
1793  // Translate sky direction into instrument direction
1794  GCTAInstDir pixinstdir = obspnt.instdir(pixdir);
1795 
1796  // Get solid angle subtended by this pixel
1797  double pixsolid = off_map->map().solidangle(pixidx);
1798 
1799  // Loop over energy bins
1800  for (int i = 0; i < nreco; ++i) {
1801 
1802  // Set event
1803  GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
1804 
1805  // Get background rate in events/s/MeV
1806  background[i] += models.eval(event, obs) * pixsolid;
1807 
1808  } // endfor: looped over energy bins
1809 
1810  } // endfor: looped over all pixels in map
1811 
1812  } // endfor: looped over all regions
1813 
1814  } // endif: IRF background template was used
1815 
1816  // ... otherwise
1817  else {
1818 
1819  // Loop over regions
1820  for (int k = 0; k < off.size(); ++k) {
1821 
1822  // Get pointer on sky region map
1823  const GSkyRegionMap* off_map =
1824  static_cast<const GSkyRegionMap*>(off[k]);
1825 
1826  // Loop over pixels in Off region map and integrate
1827  // background rate
1828  for (int j = 0; j < off_map->nonzero_indices().size(); ++j) {
1829 
1830  // Get pixel index
1831  int pixidx = off_map->nonzero_indices()[j];
1832 
1833  // Get direction to pixel center
1834  GSkyDir pixdir = off_map->map().inx2dir(pixidx);
1835 
1836  // Translate sky direction into instrument direction
1837  GCTAInstDir pixinstdir = obspnt.instdir(pixdir);
1838 
1839  // Get solid angle subtended by this pixel
1840  double pixsolid = off_map->map().solidangle(pixidx);
1841 
1842  // Loop over energy bins
1843  for (int i = 0; i < nreco; ++i) {
1844 
1845  // Get background rate in events/s/MeV
1846  background[i] += pixsolid;
1847 
1848  } // endfor: looped over energy bins
1849 
1850  } // endfor: looped over all pixels in map
1851 
1852  } // endfor: looped over all regions
1853 
1854  } // endelse: no IRF background template was used
1855 
1856  // Append background vector to Off spectrum
1857  m_off_spec.append("BACKRESP", background);
1858 
1859  } // endif: there were spectral bins
1860 
1861  // Return
1862  return;
1863 }
1864 
1865 
1866 /***********************************************************************//**
1867  * @brief Compute vector of alpha parameters
1868  *
1869  * @param[in] obs_on On CTA observation.
1870  * @param[in] obs_off Off CTA observation.
1871  * @param[in] on On regions.
1872  * @param[in] off Off regions.
1873  * @param[in] models Models including background model.
1874  * @param[in] use_model_bkg Use model background.
1875  *
1876  * @exception GException::invalid_argument
1877  * Observation does not contain relevant response or background
1878  * information
1879  *
1880  * Compute the \f$\alpha\f$ parameters for all reconstructed energy bins.
1881  * The \f$\alpha\f$ parameter gives the ratio between the On and Off region
1882  * background acceptance multiplied by the On and Off region solid angles.
1883  *
1884  * If @p use_model_bkg is @c true, the IRF background template will be used
1885  * for the computation, and \f$\alpha(E_{\rm reco})\f$ is given by
1886  *
1887  * \f[
1888  * \alpha(E_{\rm reco}) =
1889  * \frac{\sum_{\rm on} \sum_p {\tt BKG}_{{\rm on},p}(E_{\rm reco})
1890  * \times \Omega_{{\rm on},p}}
1891  * {\sum_{\rm off} \sum_p {\tt BKG}_{{\rm off},p}(E_{\rm reco})
1892  * \times \Omega_{{\rm off},p}}
1893  * \f]
1894  *
1895  * where the nominator sums over the On regions, indicated by the index
1896  * \f${\rm on}\f$, and the denominator sums over the Off regions, indicated
1897  * by the index \f${\rm off}\f$. Each On or Off region is defined by a
1898  * region sky map of type GSkyRegionMap, and the pixels of these maps
1899  * are index by \f$p\f$.
1900  * \f${\tt BKG}_{{\rm on/off},p}(E_{\rm reco})\f$ is the background rate as
1901  * provided by the IRF background template in units of
1902  * \f${\rm events} \, {\rm MeV}^{-1} \, {\rm s}^{-1} \, {\rm sr}^{-1}\f$
1903  * for a reconstructed energy \f$E_{\rm reco}\f$ and a pixel index \f$p\f$
1904  * in the On or Off region \f${\rm on/off}\f$.
1905  * \f$\Omega_{{\rm on/off},p}\f$ is the solid angle in units of
1906  * \f${\rm sr}\f$ of the pixel index \f$p\f$ in the On or Off region
1907  * \f${\rm on/off}\f$.
1908  *
1909  * If @p use_model_bkg is @c false, the background acceptance is assumed
1910  * constant and hence cancels out, which makes the \f$\alpha\f$ parameter
1911  * independent of reconstructed energy \f$E_{\rm reco}\f$.
1912  * The \f$\alpha\f$ parameter is then given by
1913  *
1914  * \f[
1915  * \alpha = \frac{\sum_{\rm on} \sum_p \Omega_{{\rm on},p}}
1916  * {\sum_{\rm off} \sum_p \Omega_{{\rm off},p}}
1917  * \f]
1918  *
1919  * @todo Compute alpha by integrating the background rate over the energy
1920  * bins and not at the bin centre.
1921  ***************************************************************************/
1923  const GCTAObservation& obs_off,
1924  const GSkyRegions& on,
1925  const GSkyRegions& off,
1926  const GModels& models,
1927  const bool& use_model_bkg)
1928 {
1929  // Get reconstructed energy boundaries from RMF
1930  const GEbounds& ereco = m_rmf.emeasured();
1931  int nreco = ereco.size();
1932 
1933  // Continue only if there are reconstructed energy bins
1934  if (nreco > 0) {
1935 
1936  // Get pointing directions of CTA On and Off observations
1937  GCTAPointing obs_on_pnt = obs_on.pointing();
1938  GCTAPointing obs_off_pnt = obs_off.pointing();
1939 
1940  // If IRF background templates shall be used then compute the
1941  // energy dependent alpha factors
1942  if (use_model_bkg) {
1943 
1944  // Loop over reconstructed energies
1945  for (int i = 0; i < nreco; ++i) {
1946 
1947  // Initialise background rate totals
1948  double aon = 0.0;
1949  double aoff = 0.0;
1950 
1951  // Loop over On regions
1952  for (int k = 0; k < on.size(); ++k) {
1953 
1954  // Get pointer on sky region map
1955  const GSkyRegionMap* on_map =
1956  static_cast<const GSkyRegionMap*>(on[k]);
1957 
1958  // Loop over pixels in On region map and integrate acceptance
1959  for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {
1960 
1961  // Get pixel index
1962  int pixidx = on_map->nonzero_indices()[j];
1963 
1964  // Get direction to pixel center
1965  GSkyDir pixdir = on_map->map().inx2dir(pixidx);
1966 
1967  // Translate sky direction into instrument direction
1968  GCTAInstDir pixinstdir = obs_on_pnt.instdir(pixdir);
1969 
1970  // Get solid angle subtended by this pixel
1971  double pixsolid = on_map->map().solidangle(pixidx);
1972 
1973  // Set event
1974  GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
1975 
1976  // Add up acceptance
1977  aon += models.eval(event, obs_on) * pixsolid;
1978 
1979  } // endfor: looped over all pixels in map
1980 
1981  } // endfor: looped over regions
1982 
1983  // Loop over Off regions
1984  for (int k = 0; k < off.size(); ++k) {
1985 
1986  // Get pointer on sky region map
1987  const GSkyRegionMap* off_map =
1988  static_cast<const GSkyRegionMap*>(off[k]);
1989 
1990  // Loop over pixels in Off region map and integrate acceptance
1991  for (int j = 0; j < off_map->nonzero_indices().size(); ++j) {
1992 
1993  // Get pixel index
1994  int pixidx = off_map->nonzero_indices()[j];
1995 
1996  // Get direction to pixel center
1997  GSkyDir pixdir = off_map->map().inx2dir(pixidx);
1998 
1999  // Translate sky direction into instrument direction
2000  GCTAInstDir pixinstdir = obs_off_pnt.instdir(pixdir);
2001 
2002  // Get solid angle subtended by this pixel
2003  double pixsolid = off_map->map().solidangle(pixidx);
2004 
2005  // Set event
2006  GCTAEventAtom event(pixinstdir, ereco.elogmean(i), GTime());
2007 
2008  // Add up acceptance
2009  aoff += models.eval(event, obs_off) * pixsolid;
2010 
2011  } // endfor: looped over all pixels in map
2012 
2013  } // endfor: looped over all regions
2014 
2015  // Compute alpha for this energy bin
2016  double alpha = (aoff > 0.0) ? aon/aoff : 1.0;
2017 
2018  // Correct for different livetime in On and Off
2019  if (obs_off.livetime() > 0.0) {
2020  alpha *= obs_on.livetime() / obs_off.livetime();
2021  }
2022 
2023  // Set background scaling in On spectra
2024  m_on_spec.backscal(i, alpha);
2025 
2026  } // endfor: looped over reconstructed energies
2027 
2028  } // endif: IRF background templates were used
2029 
2030  // ... otherwise compute energy independent alpha factor
2031  else {
2032 
2033  // Initialise background rate totals
2034  double aon = 0.0;
2035  double aoff = 0.0;
2036 
2037  // Loop over On regions
2038  for (int k = 0; k < on.size(); ++k) {
2039 
2040  // Get pointer on sky region map
2041  const GSkyRegionMap* on_map =
2042  static_cast<const GSkyRegionMap*>(on[k]);
2043 
2044  // Loop over pixels in On region map and integrate acceptance
2045  for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {
2046 
2047  // Get pixel index
2048  int pixidx = on_map->nonzero_indices()[j];
2049 
2050  // Get direction to pixel center
2051  GSkyDir pixdir = on_map->map().inx2dir(pixidx);
2052 
2053  // Translate sky direction into instrument direction
2054  GCTAInstDir pixinstdir = obs_on_pnt.instdir(pixdir);
2055 
2056  // Add up solid angle subtended by this pixel
2057  aon += on_map->map().solidangle(pixidx);
2058 
2059  } // endfor: looped over all pixels in map
2060 
2061  } // endfor: looped over regions
2062 
2063  // Loop over Off regions
2064  for (int k = 0; k < off.size(); ++k) {
2065 
2066  // Get pointer on sky region map
2067  const GSkyRegionMap* off_map =
2068  static_cast<const GSkyRegionMap*>(off[k]);
2069 
2070  // Loop over pixels in Off region map and integrate acceptance
2071  for (int j = 0; j < off_map->nonzero_indices().size(); ++j) {
2072 
2073  // Get pixel index
2074  int pixidx = off_map->nonzero_indices()[j];
2075 
2076  // Get direction to pixel center
2077  GSkyDir pixdir = off_map->map().inx2dir(pixidx);
2078 
2079  // Translate sky direction into instrument direction
2080  GCTAInstDir pixinstdir = obs_off_pnt.instdir(pixdir);
2081 
2082  // Add up solid angle subtended by this pixel
2083  aoff += off_map->map().solidangle(pixidx);
2084 
2085  } // endfor: looped over all pixels in map
2086 
2087  } // endfor: looped over all regions
2088 
2089  // Compute alpha from acceptance for this energy bin
2090  double alpha = (aoff > 0.0) ? aon/aoff : 1.0;
2091 
2092  // Correct for different livetime in On and Off
2093  if (obs_off.livetime() > 0.0) {
2094  alpha *= obs_on.livetime() / obs_off.livetime();
2095  }
2096 
2097  // Set background scaling in On spectra
2098  for (int i = 0; i < nreco; ++i) {
2099  m_on_spec.backscal(i, alpha);
2100  }
2101 
2102  } // endelse: computed energy independent alpha factor
2103 
2104  } // endif: there were energy bins
2105 
2106  // Return
2107  return;
2108 }
2109 
2110 
2111 /***********************************************************************//**
2112  * @brief Compute RMF of On/Off observation
2113  *
2114  * @param[in] obs CTA observation.
2115  * @param[in] on On regions.
2116  *
2117  * @exception GException::invalid_argument
2118  * Observation does not contain IRF response
2119  *
2120  * Compute the energy redistribution matrix for an On/Off observation. The
2121  * method requires that the RMF energy axes have been defined before.
2122  ***************************************************************************/
2124  const GSkyRegions& on)
2125 {
2126  // Get true and reconstructed energy boundaries from RMF
2127  const GEbounds& etrue = m_rmf.etrue();
2128  const GEbounds& ereco = m_rmf.emeasured();
2129  int ntrue = etrue.size();
2130  int nreco = ereco.size();
2131 
2132  // Continue only if there are RMF bins
2133  if (ntrue > 0 && nreco > 0) {
2134 
2135  // Get CTA IRF response
2137 
2138  // Get CTA observation pointing direction, zenith, and azimuth
2139  GCTAPointing obspnt = obs.pointing();
2140  GSkyDir obsdir = obspnt.dir();
2141  double zenith = obspnt.zenith();
2142  double azimuth = obspnt.azimuth();
2143 
2144  // Initialise RMF matrix
2145  for (int itrue = 0; itrue < ntrue; ++itrue) {
2146  for (int ireco = 0; ireco < nreco; ++ireco) {
2147  m_rmf(itrue, ireco) = 0.0;
2148  }
2149  }
2150 
2151  // Initialise weight matrix
2152  GMatrixSparse weight(ntrue, nreco);
2153 
2154  // Loop over On regions
2155  for (int k = 0; k < on.size(); ++k) {
2156 
2157  // Get pointer on sky region map
2158  const GSkyRegionMap* on_map = static_cast<const GSkyRegionMap*>(on[k]);
2159 
2160  // Loop over pixels in On region map and integrate acceptance
2161  for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {
2162 
2163  // Get pixel index
2164  int pixidx = on_map->nonzero_indices()[j];
2165 
2166  // Get direction to pixel center
2167  GSkyDir pixdir = on_map->map().inx2dir(pixidx);
2168 
2169  // Compute position of pixel centre in instrument coordinates
2170  double theta = obsdir.dist(pixdir);
2171  double phi = obsdir.posang(pixdir); // Celestial system
2172 
2173  // Loop over true energy
2174  for (int itrue = 0; itrue < ntrue; ++itrue) {
2175 
2176  // Compute log10 of true energy in TeV
2177  double logEtrue = etrue.elogmean(itrue).log10TeV();
2178 
2179  // Get effective area for weighting
2180  double aeff = rsp.aeff(theta, phi,
2181  zenith, azimuth,
2182  logEtrue);
2183 
2184  // Loop over reconstructed energy
2185  for (int ireco = 0; ireco < nreco; ++ireco) {
2186 
2187  // Get RMF value
2188  double value = rsp.edisp()->prob_erecobin(ereco.emin(ireco),
2189  ereco.emax(ireco),
2190  etrue.elogmean(itrue),
2191  theta);
2192 
2193  // Update RMF value and weight
2194  m_rmf(itrue, ireco) += value * aeff;
2195  weight(itrue, ireco) += aeff;
2196 
2197  } // endfor: looped over reconstructed energy
2198 
2199  } // endfor: looped over true energy
2200 
2201  } // endfor: looped over all pixels in map
2202 
2203  } // endfor: looped over all regions
2204 
2205  // Normalise RMF matrix
2206  for (int itrue = 0; itrue < ntrue; ++itrue) {
2207  for (int ireco = 0; ireco < nreco; ++ireco) {
2208  if (weight(itrue, ireco) > 0.0) {
2209  m_rmf(itrue, ireco) /= weight(itrue, ireco);
2210  }
2211  }
2212  }
2213 
2214  } // endif: there were energy bins
2215 
2216  // Return
2217  return;
2218 }
2219 
2220 
2221 /***********************************************************************//**
2222  * @brief Apply CTA observation energy boundaries to On/Off observation
2223  *
2224  * @param[in] obs CTA observation.
2225  *
2226  * Applies CTA energy boundaries to all histograms used for the On/Off
2227  * analysis in.
2228  *
2229  * For the PHA On and Off spectra, all bins are set to zero that do not fully
2230  * overlap with the CTA observation energy boundaries. Specifically,
2231  * partially overlapping bins are excluded. Some margin is applied that
2232  * effectively reduces the width of the PHA energy bin, which should cope
2233  * with any rounding errors.
2234  *
2235  * For the background response, stored in the Off spectrum, all bins are set
2236  * to zero that do not fully overlap with the CTA observation energy
2237  * boundaries.
2238  *
2239  * For the RMF, all reconstructued energy bins are set to zero that do not
2240  * fully overlap with the CTA observation energy boundaries. True energy
2241  * bins remain unchanged to properly account for energy migration.
2242  ***************************************************************************/
2244 {
2245  // Set energy margin
2246  const GEnergy energy_margin(0.01, "GeV");
2247 
2248  // Get true and reconstructed energy boundaries from RMF
2249  const GEbounds& etrue = m_rmf.etrue();
2250  const GEbounds& ereco = m_rmf.emeasured();
2251  int ntrue = etrue.size();
2252  int nreco = ereco.size();
2253 
2254  // Get energy boundaries of observations
2255  GEbounds eobs = obs.ebounds();
2256 
2257  // Get reference to background response
2258  std::vector<double>& background = m_off_spec["BACKRESP"];
2259 
2260  // Apply energy boundaries in reconstructed energy
2261  for (int ireco = 0; ireco < nreco; ++ireco) {
2262  if (!eobs.contains(ereco.emin(ireco) + energy_margin,
2263  ereco.emax(ireco) - energy_margin)) {
2264  m_on_spec[ireco] = 0.0;
2265  m_off_spec[ireco] = 0.0;
2266  background[ireco] = 0.0;
2267  for (int itrue = 0; itrue < ntrue; ++itrue) {
2268  m_rmf(itrue, ireco) = 0.0;
2269  }
2270  }
2271  }
2272 
2273  // Return
2274  return;
2275 }
2276 
2277 
2278 /***********************************************************************//**
2279  * @brief Check if effective area IRF has a radius cut
2280  *
2281  * @param[in] obs CTA observation.
2282  * @param[in] on On regions.
2283  * @return Radius cut value in degrees (zero for no radius cut).
2284  *
2285  * @exception GException::invalid_argument
2286  * IRF has a radius cut that is incompatible with the On regions
2287  *
2288  * Checks if the effective area IRF has a radius cut. If a radius cut is
2289  * found the cut value is checked against the radii of the On regions. If
2290  * the On regions are not circular, or if they have a radius that differs
2291  * from the IRF cut value, an exception is thrown.
2292  ***************************************************************************/
2294  const GSkyRegions& on) const
2295 {
2296  // Initialise radius cut value
2297  double rad_max = 0.0;
2298 
2299  // Get pointer on CTA IRF response. Continue only if a valid response
2300  // was found
2301  const GCTAResponseIrf* rsp =
2302  dynamic_cast<const GCTAResponseIrf*>(obs.response());
2303  if (rsp != NULL) {
2304 
2305  // Get pointer to CTA 2D effective area. Continue only if a 2D
2306  // effective area was found
2307  const GCTAAeff2D* aeff = dynamic_cast<const GCTAAeff2D*>(rsp->aeff());
2308  if (aeff != NULL) {
2309 
2310  // Get cut radius. Continue only if cut radius is positive
2311  rad_max = aeff->rad_max();
2312  if (rad_max > 0.0) {
2313 
2314  // Verify that all On regions are circular regions and make
2315  // sure that they have, within some margin, the same radius
2316  // than the cut radius.
2317  for (int i = 0; i < on.size(); ++i) {
2318 
2319  // Check that region is a circular region
2320  const GSkyRegionCircle* region =
2321  dynamic_cast<const GSkyRegionCircle*>(on[i]);
2322  if (region == NULL) {
2323  std::string msg = "An effective area IRF with a theta "
2324  "cut was specified, but the On region "
2325  "is not a circle, hence the "
2326  "consistency of the event selection "
2327  "could not be checked. Please specify "
2328  "a circular On region if an effective "
2329  "area with a theta cut is used.";
2331  }
2332 
2333  // Check that the cut radius is consistent
2334  if (std::abs(region->radius()-rad_max) > 1.0e-3) {
2335  std::string msg = "An effective area IRF with a theta "
2336  "cut of "+gammalib::str(rad_max)+ " deg "
2337  "was specified but an On region with "
2338  "a radius of "+
2339  gammalib::str(region->radius())+" deg "
2340  "was encountered. Please specify On "
2341  "regions with a radius of "+
2342  gammalib::str(rad_max)+" deg for this "
2343  "IRF.";
2345  }
2346 
2347  } // endfor: looped over On regions
2348 
2349  } // endif: cut radius was positive
2350 
2351  } // endif: 2D effective area was found
2352 
2353  } // endif: CTA response was found
2354 
2355  // Return radius cut flag
2356  return rad_max;
2357 }
2358 
2359 
2360 /***********************************************************************
2361  * @brief Compute \f$N_{\gamma}\f$ value and model parameter gradients
2362  *
2363  * @param[in] models Model container.
2364  * @param[in] ibin Energy bin number.
2365  * @param[in,out] grad Model gradient vector.
2366  * @returns Predicted number of source events.
2367  *
2368  * @exception GException::invalid_value
2369  * Source model is not a point source model
2370  *
2371  * Returns the predicted number of source events \f$N_{\gamma}\f$
2372  * in the On regions for a given energy bin. The method computes also
2373  *
2374  * \f[
2375  * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
2376  * \f]
2377  *
2378  * which are the gradients in the predicted number of source events with
2379  * respect to all model parameters.
2380  ***********************************************************************/
2382  const int& ibin,
2383  GVector* grad) const
2384 {
2385  // Get total number of model parameters
2386  int npars = models.npars();
2387 
2388  // Initialize results
2389  double value = 0.0;
2390  for (int i = 0; i < npars; ++i) {
2391  (*grad)[i] = 0.0;
2392  }
2393 
2394  // Continue only if bin number is in range and there are model parameters
2395  if ((ibin >= 0) && (ibin < m_on_spec.size()) && (npars > 0)) {
2396 
2397  // Initialise parameter index
2398  int ipar = 0;
2399 
2400  // Loop over models
2401  for (int j = 0; j < models.size(); ++j) {
2402 
2403  // Get model pointer. Fall through if pointer is not valid
2404  const GModel* mptr = models[j];
2405  if (mptr == NULL) {
2406  continue;
2407  }
2408 
2409  // Fall through if model does not apply to specific instrument
2410  // and observation identifier
2411  if (!mptr->is_valid(instrument(), id())) {
2412  ipar += mptr->size();
2413  continue;
2414  }
2415 
2416  // Fall through if this model component is a not sky component
2417  const GModelSky* sky = dynamic_cast<const GModelSky*>(mptr);
2418  if (sky == NULL) {
2419  ipar += mptr->size();
2420  continue;
2421  }
2422 
2423  // Spectral component (the useful one)
2424  GModelSpectral* spectral = sky->spectral();
2425  if (spectral != NULL) {
2426 
2427  // Debug code
2428  #if defined(G_N_GAMMA_DEBUG)
2429  double rmf_sum = 0.0;
2430  #endif
2431 
2432  // Get indices of spectral gradients. We need this to update
2433  // only the spectral parameter gradients later.
2434  std::vector<int> inx_spec;
2435  for (int k = 0; k < mptr->size(); ++k) {
2436  const GModelPar& par = (*mptr)[k];
2437  if (spectral->has_par(par.name())) {
2438  inx_spec.push_back(k);
2439  }
2440  }
2441 
2442  // Set instrument scale factor
2443  double scale = (sky->has_scales()) ? sky->scale(instrument()).value() : 1.0;
2444 
2445  // Loop over true energy bins
2446  for (int itrue = 0; itrue < m_arf.size(); ++itrue) {
2447 
2448  // Get ARF value. Continue only if it is positive
2449  double arf = m_arf[itrue];
2450  if (arf <= 0.0) {
2451  continue;
2452  }
2453 
2454  // Get RMF value. Continue only if it is positive
2455  double rmf = m_rmf(itrue, ibin);
2456  if (rmf <= 0.0) {
2457  continue;
2458  }
2459 
2460  // Debug code
2461  #if defined(G_N_GAMMA_DEBUG)
2462  rmf_sum += rmf;
2463  #endif
2464 
2465  // Get true energy bin properties
2466  GEnergy etruemin = m_arf.ebounds().emin(itrue);
2467  GEnergy etruemax = m_arf.ebounds().emax(itrue);
2468  GEnergy etruemean = m_arf.ebounds().elogmean(itrue);
2469  double etruewidth = m_arf.ebounds().ewidth(itrue).MeV();
2470 
2471  // Compute normalisation factors. The gradient normalisation
2472  // factor uses the energy bin width, since the eval() method
2473  // returns a differential flux
2474  double exposure = m_on_spec.exposure();
2475  double norm_scale = arf * exposure * rmf;
2476  double norm_flux = norm_scale * scale;
2477  double norm_grad = norm_flux * etruewidth;
2478 
2479  // Determine number of gamma-ray events in model by
2480  // computing the flux over the true energy bin in
2481  // ph/cm2/s and multiplying this by effective area (cm2),
2482  // livetime (s), redistribution probability and scale
2483  double flux = spectral->flux(etruemin, etruemax);
2484  value += flux * norm_flux;
2485 
2486  // Determine the model gradients at the current true
2487  // energy. The eval() method needs a time in case that the
2488  // spectral model has a time dependence. We simply use a
2489  // dummy time here.
2490  spectral->eval(etruemean, GTime(), true);
2491 
2492  // Compute the parameter gradients for spectral parameters
2493  for (int k = 0; k < inx_spec.size(); ++k) {
2494  const GModelPar& par = (*mptr)[inx_spec[k]];
2495  if (par.is_free() && ipar+k < npars) {
2496  (*grad)[ipar+inx_spec[k]] += par.factor_gradient() * norm_grad;
2497  }
2498  }
2499 
2500  // Optionally compute scaling parameter gradient
2501  if (sky->has_scales()) {
2502  for (int k = 0; k < mptr->size(); ++k) {
2503  const GModelPar& par = (*mptr)[k];
2504  if ((par.name() == instrument()) && par.is_free() &&
2505  (ipar+k < npars)) {
2506  (*grad)[ipar+k] += par.scale() * flux * norm_scale;
2507  }
2508  }
2509  }
2510 
2511  } // endfor: looped over true energy bins
2512 
2513  // Debug code
2514  #if defined(G_N_GAMMA_DEBUG)
2515  std::cout << "sum(RMF) = " << rmf_sum << std::endl;
2516  #endif
2517 
2518  } // endif: pointer to spectral component was not NULL
2519 
2520  // Increment the parameter index
2521  ipar += mptr->size();
2522 
2523  } // endfor: looped over model components
2524 
2525  } // endif: bin number is in the range and model container is not empty
2526 
2527  // Return number of gamma-ray events
2528  return value;
2529 }
2530 
2531 
2532 /***********************************************************************
2533  * @brief Compute \f$N_{\rm bgd}\f$ value and model parameter gradients
2534  *
2535  * @param[in] models Model container.
2536  * @param[in] ibin Energy bin index.
2537  * @param[in,out] grad Model gradient vector.
2538  * @return Predicted number of background events in Off regions.
2539  *
2540  * Returns the predicted number of background events \f$N_{\rm bgd}\f$
2541  * in the Off regions for a given energy bin. The method computes also
2542  *
2543  * \f[
2544  * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
2545  * \f]
2546  *
2547  * which are the gradients in the predicted number of background events
2548  * with respect to all model parameters.
2549  ***********************************************************************/
2551  const int& ibin,
2552  GVector* grad) const
2553 {
2554  // Get total number of model parameters
2555  int npars = models.npars();
2556 
2557  // Initialize results
2558  double value = 0.0;
2559  for (int i = 0; i < npars; ++i) {
2560  (*grad)[i] = 0.0;
2561  }
2562 
2563  // Continue only if bin number is valid and if there are model parameters
2564  if ((ibin >= 0) && (ibin < m_on_spec.size()) && (npars > 0)) {
2565 
2566  // Get reconstructed energy bin mean and width
2567  GEnergy emean = m_on_spec.ebounds().elogmean(ibin);
2568  double ewidth = m_on_spec.ebounds().ewidth(ibin).MeV();
2569 
2570  // Perform log-log interpolation of background rate (events/MeV/s)
2571  // at reconstructed energy
2572  double background = m_off_spec["BACKRESP"][ibin];
2573 
2574  // Continue only if background rate is positive
2575  if (background > 0.0) {
2576 
2577  // Compute normalisation factor (events)
2578  double exposure = m_off_spec.exposure();
2579  double norm = background * exposure * ewidth;
2580 
2581  // Loop over models
2582  for (int j = 0, ipar = 0; j < models.size(); ++j) {
2583 
2584  // Get model pointer. Fall through if pointer is not valid
2585  const GModel* mptr = models[j];
2586  if (mptr == NULL) {
2587  continue;
2588  }
2589 
2590  // Fall through if model does not apply to specific instrument
2591  // and observation identifier.
2592  if (!mptr->is_valid(this->instrument(), this->id())) {
2593  ipar += mptr->size();
2594  continue;
2595  }
2596 
2597  // Fall through if model is not a background component
2598  const GModelData* bgd = dynamic_cast<const GModelData*>(mptr);
2599  if (bgd == NULL) {
2600  ipar += mptr->size();
2601  continue;
2602  }
2603 
2604  // Extract model dependent spectral component
2605  const GModelSpectral* spectral = gammalib::cta_model_spectral(*bgd);
2606 
2607  // Get value from spectral component
2608  if (spectral != NULL) {
2609 
2610  // Determine the number of background events in model by
2611  // computing the model normalization at the mean value of
2612  // the energy bin and multiplying the normalisation with
2613  // the number of background events. The eval() method needs
2614  // a time in case that the spectral model has a time
2615  // dependence. We simply use a dummy time here.
2616  value += spectral->eval(emean, GTime(), true) * norm;
2617 
2618  // Compute the parameter gradients for all model parameters
2619  for (int k = 0; k < mptr->size(); ++k) {
2620  const GModelPar& par = (*mptr)[k];
2621  if (par.is_free() && ipar+k < npars) {
2622  (*grad)[ipar+k] += par.factor_gradient() * norm;
2623  }
2624  }
2625 
2626  } // endif: pointer to spectral component was not NULL
2627 
2628  // Increment the parameter index
2629  ipar += mptr->size();
2630 
2631  } // endfor: looped over model components
2632 
2633  } // endif: background rate was positive
2634 
2635  } // endif: bin number is in the range and model container is not empty
2636 
2637  // Return
2638  return value;
2639 }
2640 
2641 
2642 /***********************************************************************//**
2643  * @brief Evaluate log-likelihood function for On/Off analysis in the
2644  * case of Poisson signal with modeled Poisson background
2645  *
2646  * @param[in] models Models.
2647  * @param[in,out] gradient Pointer to gradients.
2648  * @param[in,out] curvature Pointer to curvature matrix.
2649  * @param[in,out] npred Pointer to Npred value.
2650  * @return Log-likelihood value.
2651  *
2652  * @exception GException::invalid_value
2653  * There are no model parameters.
2654  *
2655  * Computes the log-likelihood value for the On/Off observation. The
2656  * method loops over the energy bins to update the function value, its
2657  * derivatives and the curvature matrix. The number of On counts
2658  * \f$N_{\rm on}\f$ and Off counts \f$N_{\rm off}\f$ are taken from the
2659  * On and Off spectra, the expected number of gamma-ray events
2660  * \f$N_{\gamma}\f$ and background events \f$N_{\rm bgd}\f$ are
2661  * computed from the spectral models of the relevant components in the
2662  * model container (spatial and temporal components are ignored so far).
2663  * See the N_gamma() and N_bgd() methods for details about the model
2664  * computations.
2665  *
2666  * The log-likelihood is given by
2667  *
2668  * \f[
2669  * \ln L = \sum_i \ln L_i
2670  * \f]
2671  *
2672  * where the sum is taken over all energy bins \f$i\f$ and
2673  *
2674  * \f[
2675  * \ln L_i = - N_{\rm on} \ln N_{\rm pred} + N_{\rm pred}
2676  * - N_{\rm off} \ln N_{\rm bgd} + N_{\rm bgd}
2677  * \f]
2678  *
2679  * with
2680  *
2681  * \f[
2682  * N_{\rm pred} = N_{\gamma} + \alpha N_{\rm bgd}
2683  * \f]
2684  *
2685  * being the total number of predicted events for an energy bin in the On
2686  * region,
2687  * \f$N_{\rm on}\f$ is the total number of observed events for an energy
2688  * bin in the On region,
2689  * \f$N_{\rm off}\f$ is the total number of observed events for an energy
2690  * bin in the Off region, and
2691  * \f$N_{\rm bgd}\f$ is the predicted number of background events for an
2692  * energy bin in the Off region.
2693  *
2694  * The log-likelihood gradient with respect to sky model parameters
2695  * \f$p_{\rm sky}\f$ is given by
2696  *
2697  * \f[
2698  * \left( \frac{\partial \ln L_i}{\partial p_{\rm sky}} \right) =
2699  * \left( 1 - \frac{N_{\rm on}}{N_{\rm pred}} \right)
2700  * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
2701  * \f]
2702  *
2703  * and with respect to background model parameters \f$p_{\rm bgd}\f$ is
2704  * given by
2705  *
2706  * \f[
2707  * \left( \frac{\partial \ln L_i}{\partial p_{\rm bgd}} \right) =
2708  * \left( 1 + \alpha - \frac{N_{\rm off}}{N_{\rm bgd}} -
2709  * \frac{\alpha N_{\rm on}}{N_{\rm pred}} \right)
2710  * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
2711  * \f]
2712  *
2713  * The curvature matrix elements are given by
2714  *
2715  * \f[
2716  * \left( \frac{\partial^2 \ln L_i}{\partial^2 p_{\rm sky}} \right) =
2717  * \left( \frac{N_{\rm on}}{N_{\rm pred}^2} \right)
2718  * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)^2
2719  * \f]
2720  *
2721  * \f[
2722  * \left( \frac{\partial^2 \ln L_i}{\partial p_{\rm sky}
2723  * \partial p_{\rm bgd}} \right) =
2724  * \left( \frac{\alpha N_{\rm on}}{N_{\rm pred}^2} \right)
2725  * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
2726  * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
2727  * \f]
2728  *
2729  * \f[
2730  * \left( \frac{\partial^2 \ln L_i}{\partial p_{\rm bgd}
2731  * \partial p_{\rm sky}} \right) =
2732  * \left( \frac{\alpha N_{\rm on}}{N_{\rm pred}^2} \right)
2733  * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
2734  * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
2735  * \f]
2736  *
2737  * \f[
2738  * \left( \frac{\partial^2 \ln L_i}{\partial^2 p_{\rm bgd}} \right) =
2739  * \left( \frac{N_{\rm off}}{N_{\rm bgd}^2} +
2740  * \frac{\alpha^2 N_{\rm on}}{N_{\rm pred}^2} \right)
2741  * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)^2
2742  * \f]
2743  ***********************************************************************/
2745  GVector* gradient,
2746  GMatrixSparse* curvature,
2747  double* npred) const
2748 {
2749  // Timing measurement
2750  #if defined(G_LIKELIHOOD_DEBUG)
2751  #ifdef _OPENMP
2752  double t_start = omp_get_wtime();
2753  #else
2754  clock_t t_start = clock();
2755  #endif
2756  #endif
2757 
2758  // Initialise statistics
2759  #if defined(G_LIKELIHOOD_DEBUG)
2760  int n_bins = m_on_spec.size();
2761  int n_used = 0;
2762  int n_small_model = 0;
2763  int n_zero_data = 0;
2764  double sum_data = 0.0;
2765  double sum_model = 0.0;
2766  #endif
2767 
2768  // Initialise log-likelihood value
2769  double value = 0.0;
2770 
2771  // Get number of model parameters in model container
2772  int npars = models.npars();
2773 
2774  // Create model gradient vectors for sky and background parameters
2775  GVector sky_grad(npars);
2776  GVector bgd_grad(npars);
2777 
2778  // Allocate working array
2779  GVector colvar(npars);
2780 
2781  // Loop over all energy bins
2782  for (int i = 0; i < m_on_spec.size(); ++i) {
2783 
2784  // Reinitialize working arrays
2785  for (int j = 0; j < npars; ++j) {
2786  sky_grad[j] = 0.0;
2787  bgd_grad[j] = 0.0;
2788  }
2789 
2790  // Get number of On and Off counts
2791  double non = m_on_spec[i];
2792  double noff = m_off_spec[i];
2793 
2794  // Get background scaling
2795  double alpha = m_on_spec.backscal(i);
2796 
2797  // Get number of gamma and background events (and corresponding
2798  // spectral model gradients)
2799  double ngam = N_gamma(models, i, &sky_grad);
2800  double nbgd = N_bgd(models, i, &bgd_grad);
2801 
2802  // Compute and updated predicted number of events
2803  double nonpred = ngam + alpha * nbgd;
2804  *npred += nonpred;
2805 
2806  // Skip bin if model is too small (avoids -Inf or NaN gradients)
2807  if ((nbgd <= minmod) || (nonpred <= minmod)) {
2808  #if defined(G_LIKELIHOOD_DEBUG)
2809  n_small_model++;
2810  #endif
2811  continue;
2812  }
2813 
2814  // Now we have all predicted gamma and background events for
2815  // current energy bin. Update the log(likelihood)
2816  value += -non * std::log(nonpred) + nonpred - noff * std::log(nbgd) + nbgd;
2817 
2818  // Update statistics
2819  #if defined(G_LIKELIHOOD_DEBUG)
2820  n_used++;
2821  sum_data += non;
2822  sum_model += nonpred;
2823  #endif
2824 
2825  // Fill derivatives
2826  double fa = non/nonpred;
2827  double fb = fa/nonpred;
2828  double fc = alpha * fb;
2829  double fd = fc * alpha + noff/(nbgd*nbgd);
2830  double sky_factor = 1.0 - fa;
2831  double bgd_factor = 1.0 + alpha - alpha * fa - noff/nbgd;
2832 
2833  // Loop over all parameters
2834  for (int j = 0; j < npars; ++j) {
2835 
2836  // If spectral model for sky component is non-zero and
2837  // non-infinite then handle sky component gradients and
2838  // second derivatives including at least a sky component ...
2839  if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
2840 
2841  // Gradient
2842  (*gradient)[j] += sky_factor * sky_grad[j];
2843 
2844  // Hessian (from first-order derivatives only)
2845  for (int k = 0; k < npars; ++k) {
2846 
2847  // If spectral model for sky component is non-zero and
2848  // non-infinite then we have the curvature element
2849  // of a sky component
2850  if (sky_grad[k] != 0.0 &&
2851  !gammalib::is_infinite(sky_grad[k])) {
2852  colvar[k] = sky_grad[j] * sky_grad[k] * fb;
2853  }
2854 
2855  // ... else if spectral model for background component
2856  // is non-zero and non-infinite then we have the mixed
2857  // curvature element between a sky and a background
2858  // component
2859  else if (bgd_grad[k] != 0.0 &&
2860  !gammalib::is_infinite(bgd_grad[k])) {
2861  colvar[k] = sky_grad[j] * bgd_grad[k] * fc;
2862  }
2863 
2864  // ...else neither sky nor background
2865  else {
2866  colvar[k] = 0.0;
2867  }
2868 
2869  } // endfor: Hessian computation
2870 
2871  // Update matrix
2872  curvature->add_to_column(j, colvar);
2873 
2874  } // endif: spectral model is non-zero and non-infinite
2875 
2876  // ... otherwise if spectral model for background component is
2877  // non-zero and non-infinite then handle background component
2878  // gradients and second derivatives including at least a
2879  // background component
2880  else if (bgd_grad[j] != 0.0 &&
2881  !gammalib::is_infinite(bgd_grad[j])) {
2882 
2883  // Gradient
2884  (*gradient)[j] += bgd_factor * bgd_grad[j];
2885 
2886  // Hessian (from first-order derivatives only)
2887  for (int k = 0; k < npars; ++k) {
2888 
2889  // If spectral model for sky component is non-zero and
2890  // non-infinite then we have the mixed curvature element
2891  // between a sky and a background component
2892  if (sky_grad[k] != 0.0 &&
2893  !gammalib::is_infinite(sky_grad[k])) {
2894  colvar[k] = bgd_grad[j] * sky_grad[k] * fc;
2895  }
2896 
2897  // ... else if spectral model for background component
2898  // is non-zero and non-infinite then we have the
2899  // curvature element of a background component
2900  else if (bgd_grad[k] != 0.0 &&
2901  !gammalib::is_infinite(bgd_grad[k])) {
2902  colvar[k] = bgd_grad[j] * bgd_grad[k] * fd;
2903  }
2904 
2905  // ... else neither sky nor background
2906  else {
2907  colvar[k] = 0.0;
2908  }
2909 
2910  } // endfor: Hessian computation
2911 
2912  // Update matrix
2913  curvature->add_to_column(j, colvar);
2914 
2915  } // endif: spectral model for background component valid
2916 
2917  } // endfor: looped over all parameters for derivatives computation
2918 
2919  } // endfor: looped over energy bins
2920 
2921  // Dump statistics
2922  #if defined(G_LIKELIHOOD_DEBUG)
2923  std::cout << "Number of parameters: " << npars << std::endl;
2924  std::cout << "Number of bins: " << n_bins << std::endl;
2925  std::cout << "Number of bins used for computation: " << n_used << std::endl;
2926  std::cout << "Number of bins excluded due to small model: ";
2927  std::cout << n_small_model << std::endl;
2928  std::cout << "Number of bins with zero data: " << n_zero_data << std::endl;
2929  std::cout << "Sum of data (On): " << sum_data << std::endl;
2930  std::cout << "Sum of model (On): " << sum_model << std::endl;
2931  std::cout << "Statistic: " << value << std::endl;
2932  #endif
2933 
2934  // Optionally dump gradient and curvature matrix
2935  #if defined(G_LIKELIHOOD_DEBUG)
2936  std::cout << *gradient << std::endl;
2937  std::cout << *curvature << std::endl;
2938  #endif
2939 
2940  // Timing measurement
2941  #if defined(G_LIKELIHOOD_DEBUG)
2942  #ifdef _OPENMP
2943  double t_elapse = omp_get_wtime()-t_start;
2944  #else
2945  double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
2946  #endif
2947  std::cout << "GCTAOnOffObservation::optimizer::likelihood_cstat: CPU usage = "
2948  << t_elapse << " sec" << std::endl;
2949  #endif
2950 
2951  // Return
2952  return value;
2953 }
2954 
2955 
2956 
2957 /***********************************************************************//**
2958  * @brief Evaluate log-likelihood function for On/Off analysis in the
2959  * case of Poisson signal with measured Poisson background
2960  *
2961  * @param[in] models Models.
2962  * @param[in,out] gradient Pointer to gradients.
2963  * @param[in,out] curvature Pointer to curvature matrix.
2964  * @param[in,out] npred Pointer to Npred value.
2965  * @return Log-likelihood value.
2966  *
2967  * Computes the log-likelihood value for the On/Off observation. The method
2968  * loops over the energy bins to update the function value, its derivatives
2969  * and the curvature matrix.
2970  *
2971  * The number of On counts \f$n_{\rm on}\f$ and Off counts \f$n_{\rm off}\f$
2972  * are taken from the On and Off spectra, the expected number of gamma-ray
2973  * events \f$\mu_s\f$ is computed from the spectral models of the relevant
2974  * components in the model container (spatial and temporal components
2975  * are ignored so far). See the N_gamma() method for details about the
2976  * model computations.
2977  *
2978  * The estimated number of background counts \f$\mu_b\f$ is derived based
2979  * on the measurement in the Off region by analytical minimization of
2980  * the Poisson likelihood, i.e., it is treated as a nuisance parameter.
2981  * See the wstat_value function for the derivation.
2982  ***********************************************************************/
2984  GVector* gradient,
2985  GMatrixSparse* curvature,
2986  double* npred) const
2987 {
2988  // Debug option: dump header
2989  #if defined(G_LIKELIHOOD_DEBUG)
2990  std::cout << "GCTAOnOffObservation::likelihood_wstat: enter" << std::endl;
2991  #endif
2992 
2993  // Timing measurement
2994  #if defined(G_LIKELIHOOD_DEBUG)
2995  #ifdef _OPENMP
2996  double t_start = omp_get_wtime();
2997  #else
2998  clock_t t_start = clock();
2999  #endif
3000  #endif
3001 
3002  // Initialise statistics
3003  #if defined(G_LIKELIHOOD_DEBUG)
3004  int n_bins = m_on_spec.size();
3005  int n_used = 0;
3006  double sum_data = 0.0;
3007  double sum_model = 0.0;
3008  #endif
3009 
3010  // Initialise log-likelihood value
3011  double value = 0.0;
3012 
3013  // Get number of model parameters in model container
3014  int npars = models.npars();
3015 
3016  // Create model gradient vectors for sky parameters
3017  GVector sky_grad(npars);
3018 
3019  // Allocate working array
3020  GVector colvar(npars);
3021 
3022  // Loop over all energy bins
3023  for (int i = 0; i < m_on_spec.size(); ++i) {
3024 
3025  // Reinitialize working arrays
3026  for (int j = 0; j < npars; ++j) {
3027  sky_grad[j] = 0.0;
3028  }
3029 
3030  // Get number of On and Off counts
3031  double non = m_on_spec[i];
3032  double noff = m_off_spec[i];
3033 
3034  // Get background scaling
3035  double alpha = m_on_spec.backscal(i);
3036 
3037  // Get number of gamma events (and corresponding spectral model
3038  // gradients)
3039  double ngam = N_gamma(models, i, &sky_grad);
3040 
3041  // Initialise variables for likelihood computation
3042  double nonpred = 0.0;
3043  double nbgd = 0.0;
3044  double dlogLdsky = 0.0;
3045  double d2logLdsky2 = 0.0;
3046 
3047  // Perform likelihood profiling and derive likelihood value
3048  double logL = wstat_value(non, noff, alpha, ngam,
3049  nonpred, nbgd, dlogLdsky, d2logLdsky2);
3050 
3051  // Update global log-likelihood and Npred
3052  value += logL;
3053  *npred += nonpred;
3054 
3055  // Update statistics
3056  #if defined(G_LIKELIHOOD_DEBUG)
3057  n_used++;
3058  sum_data += non;
3059  sum_model += nonpred;
3060  #endif
3061 
3062  // Loop over all parameters
3063  for (int j = 0; j < npars; ++j) {
3064 
3065  // If spectral model for sky component is non-zero and
3066  // non-infinite then handle sky component gradients and
3067  // second derivatives including at least a sky component ...
3068  if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
3069 
3070  // Gradient
3071  (*gradient)[j] += dlogLdsky * sky_grad[j];
3072 
3073  // Hessian (from first-order derivatives only)
3074  for (int k = 0; k < npars; ++k) {
3075 
3076  // If spectral model for sky component is non-zero and
3077  // non-infinite then we have the curvature element
3078  // of a sky component
3079  if (sky_grad[k] != 0.0 &&
3080  !gammalib::is_infinite(sky_grad[k])) {
3081  colvar[k] = sky_grad[j] * sky_grad[k] * d2logLdsky2;
3082  }
3083 
3084  // ...... else if spectral model for background component
3085  // or neither sky nor background
3086  else {
3087  colvar[k] = 0.0;
3088  }
3089 
3090  } // endfor: Hessian computation
3091 
3092  // Update matrix
3093  curvature->add_to_column(j, colvar);
3094 
3095  } // endif: spectral model is non-zero and non-infinite
3096 
3097  } // endfor: looped over all parameters for derivatives computation
3098 
3099  } // endfor: looped over energy bins
3100 
3101  // Dump statistics
3102  #if defined(G_LIKELIHOOD_DEBUG)
3103  std::cout << "Number of parameters: " << npars << std::endl;
3104  std::cout << "Number of bins: " << n_bins << std::endl;
3105  std::cout << "Number of bins used for computation: " << n_used << std::endl;
3106  std::cout << "Sum of data (On): " << sum_data << std::endl;
3107  std::cout << "Sum of model (On): " << sum_model << std::endl;
3108  std::cout << "Statistic: " << value << std::endl;
3109  #endif
3110 
3111  // Optionally dump gradient and curvature matrix
3112  #if defined(G_LIKELIHOOD_DEBUG)
3113  std::cout << *gradient << std::endl;
3114  std::cout << *curvature << std::endl;
3115  #endif
3116 
3117  // Timing measurement
3118  #if defined(G_LIKELIHOOD_DEBUG)
3119  #ifdef _OPENMP
3120  double t_elapse = omp_get_wtime()-t_start;
3121  #else
3122  double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
3123  #endif
3124  std::cout << "GCTAOnOffObservation::optimizer::likelihood_wstat: CPU usage = "
3125  << t_elapse << " sec" << std::endl;
3126  #endif
3127 
3128  // Debug option: dump trailer
3129  #if defined(G_LIKELIHOOD_DEBUG)
3130  std::cout << "GCTAOnOffObservation::likelihood_wstat: exit" << std::endl;
3131  #endif
3132 
3133  // Return
3134  return value;
3135 }
3136 
3137 
3138 /***********************************************************************//**
3139  * @brief Evaluate log-likelihood value in energy bin for On/Off analysis
3140  * by profiling over the number of background counts
3141  *
3142  * @param[in] non number of On counts
3143  * @param[in] noff number of Off counts
3144  * @param[in] alpha background scaling rate
3145  * @param[in] ngam number of predicted gamma-ray counts
3146  * @param[in,out] nonpred number of predicted On counts
3147  * @param[in,out] nbgd number of predicted background counts
3148  * @param[in,out] dlogLdsky first derivative of log-like w.r.t. sky pars
3149  * @param[in,out] d2logLdsky2 second derivative of log-like w.r.t. sky pars
3150  * @return Log-likelihood value.
3151  *
3152  * Computes the log-likelihood value for the On/Off observation in an energy bin
3153  * by treating the number of background counts as nuisance parameter. The method
3154  * performs an analytical minimisation of the Poisson likelihood and updates
3155  * the relevant values.
3156  * In the general case, the log-likelihood function is computed using
3157  *
3158  * \f[
3159  * L = \mu_s + (1+\alpha) \mu_b - n_{\rm on} - n_{\rm off} -
3160  * n_{\rm on} \left( \ln(\mu_s + \alpha \mu_b) - \ln(n_{\rm on})
3161  * \right) -
3162  * n_{\rm off} \left( \ln(\mu_b) - \ln(n_{\rm off}) \right)
3163  * \f]
3164  *
3165  * where
3166  *
3167  * \f[
3168  * \mu_b = \frac{C+D}{2\alpha(1+\alpha)}
3169  * \f]
3170  *
3171  * are the estimated number of background counts with
3172  *
3173  * \f[
3174  * C = \alpha (n_{\rm on} + n_{\rm off}) - (1 + \alpha) \mu_s
3175  * \f]
3176  *
3177  * and
3178  *
3179  * \f[
3180  * D = \sqrt{C^2 + 4 (1 + \alpha) \, \alpha \, n_{\rm off} \, \mu_s}
3181  * \f]
3182  *
3183  * \
3184  *
3185  * The first derivative of the log-likelihood function is given by
3186  *
3187  * \f[
3188  * \frac{\delta L}{\delta \mu_s} =
3189  * 1 + (1 + \alpha) \frac{\delta \mu_b}{\delta \mu_s} -
3190  * \frac{n_{\rm on}}{\mu_s + \alpha \mu_b} \left( 1 + \alpha
3191  * \frac{\delta \mu_b}{\delta \mu_s} \right) -
3192  * \frac{n_{\rm off}}{\mu_b} \frac{\delta \mu_b}{\delta \mu_s}
3193  * \f]
3194  *
3195  * with
3196  *
3197  * \f[
3198  * \frac{\delta \mu_b}{\delta \mu_s} =
3199  * \frac{n_{\rm off} - C}{D} - \frac{1}{2 \alpha}
3200  * \f]
3201  *
3202  * The second derivative of the log-likelihood function is given by
3203  *
3204  * \f[
3205  * \frac{\delta^2 L}{\delta \mu_s^2} =
3206  * \frac{n_{\rm on}}{(\mu_s + \alpha \mu_b)^2} \left( 1 + \alpha
3207  * \frac{\delta \mu_b}{\delta \mu_s} \right) +
3208  * \frac{\delta^2 \mu_b}{\delta \mu_s^2} \left( (1 + \alpha) -
3209  * \frac{\alpha n_{\rm on}}{\mu_s + \alpha \mu_b} -
3210  * \frac{n_{\rm off}}{\mu_b} \right) +
3211  * \frac{\delta \mu_b}{\delta \mu_s} \left(
3212  * \frac{\alpha n_{\rm on}}{(\mu_s + \alpha \mu_b)^2} \left(
3213  * 1 + \alpha \frac{\delta \mu_b}{\delta \mu_s} \right) +
3214  * \frac{n_{\rm off}}{\mu_b^2} \frac{\delta \mu_b}{\delta \mu_s}
3215  * \right)
3216  * \f]
3217  *
3218  * with
3219  *
3220  * \f[
3221  * \frac{\delta^2 \mu_b}{\delta \mu_s^2} =
3222  * \frac{-1}{2 \alpha} \left(
3223  * \frac{1}{D} \frac{\delta C}{\delta \mu_s} +
3224  * \frac{2 \alpha \, n_{\rm off} - C}{D^2} \frac{\delta D}{\delta \mu_s}
3225  * \right)
3226  * \f]
3227  *
3228  * where
3229  *
3230  * \f[
3231  * \frac{\delta C}{\delta \mu_s} = -(1 + \alpha)
3232  * \f]
3233  *
3234  * and
3235  *
3236  * \f[
3237  * \frac{\delta D}{\delta \mu_s} =
3238  * \frac{4 (1 + \alpha) \, \alpha \, n_{\rm off} - 2 \, C \, (1 + \alpha)}
3239  * {2D}
3240  * \f]
3241  *
3242  * In the special case \f$n_{\rm on}=n_{\rm off}=0\f$ the formal
3243  * background estimate becomes negative, hence we set \f$\mu_b=0\f$ and
3244  * the log-likelihood function becomes
3245  *
3246  * \f[
3247  * L = \mu_s
3248  * \f]
3249  *
3250  * the first derivative
3251  *
3252  * \f[
3253  * \frac{\delta L}{\delta \mu_s} = 1
3254  * \f]
3255  *
3256  * and the second derivative
3257  *
3258  * \f[
3259  * \frac{\delta^2 L}{\delta \mu_s^2} = 0
3260  * \f]
3261  *
3262  * In the special case \f$n_{\rm on}=0\f$ and \f$n_{\rm off}>0\f$
3263  * the log-likelihood function becomes
3264  *
3265  * \f[
3266  * L = \mu_s + n_{\rm off} \ln(1 + \alpha)
3267  * \f]
3268  *
3269  * the background estimate
3270  *
3271  * \f[
3272  * \mu_b = \frac{n_{\rm off}}{1+\alpha}
3273  * \f]
3274  *
3275  * the first derivative
3276  *
3277  * \f[
3278  * \frac{\delta L}{\delta \mu_s} = 1
3279  * \f]
3280  *
3281  * and the second derivative
3282  *
3283  * \f[
3284  * \frac{\delta^2 L}{\delta \mu_s^2} = 0
3285  * \f]
3286  *
3287  * In the special case \f$n_{\rm on}>0\f$ and \f$n_{\rm off}=0\f$
3288  * the background estimate becomes
3289  *
3290  * \f[
3291  * \mu_b = \frac{n_{\rm on}}{1+\alpha} - \frac{\mu_s}{\alpha}
3292  * \f]
3293  *
3294  * which is positive for
3295  *
3296  * \f[
3297  * \mu_s < n_{\rm on} \frac{\alpha}{1+\alpha}
3298  * \f]
3299  *
3300  * For positive \f$\mu_b\f$ the log-likelihood function is given by
3301  *
3302  * \f[
3303  * L = -\frac{\mu_s}{\alpha}
3304  * - n_{\rm on} \ln \left(\frac{\alpha}{1 + \alpha} \right)
3305  * \f]
3306  *
3307  * the first derivative
3308  *
3309  * \f[
3310  * \frac{\delta L}{\delta \mu_s} = -\frac{1}{\alpha}
3311  * \f]
3312  *
3313  * and the second derivative
3314  *
3315  * \f[
3316  * \frac{\delta^2 L}{\delta \mu_s^2} = 0
3317  * \f]
3318  *
3319  * For negative \f$\mu_b\f$ we set \f$\mu_b=0\f$ and the log-likelihood
3320  * function becomes
3321  *
3322  * \f[
3323  * L = \mu_s - n_{\rm on} -
3324  * n_{\rm on} \left( \ln(\mu_s) - \ln(n_{\rm on}) \right)
3325  * \f]
3326  *
3327  * the first derivative
3328  *
3329  * \f[
3330  * \frac{\delta L}{\delta \mu_s} = 1 - \frac{n_{\rm on}}{\mu_s}
3331  * \f]
3332  *
3333  * and the second derivative
3334  *
3335  * \f[
3336  * \frac{\delta^2 L}{\delta \mu_s^2} = \frac{1}{\mu_s^2}
3337  * \f]
3338  *
3339  * Note that the fit results may be biased and the statistical errors
3340  * overestimated if for some bins \f$n_{\rm on}=0\f$ and/or
3341  * \f$n_{\rm off}=0\f$ (i.e. if the special cases are encountered).
3342  ***********************************************************************/
3343 double GCTAOnOffObservation::wstat_value(const double& non,
3344  const double& noff,
3345  const double& alpha,
3346  const double& ngam,
3347  double& nonpred,
3348  double& nbgd,
3349  double& dlogLdsky,
3350  double& d2logLdsky2) const
3351 {
3352  // Initialise log-likelihood value
3353  double logL;
3354 
3355  // Precompute some values
3356  double alphap1 = alpha + 1.0;
3357  double alpharat = alpha / alphap1;
3358 
3359  // Calculate number of background events, profile likelihood value
3360  // and likelihood derivatives
3361  // Special case noff = 0
3362  if (noff == 0.0) {
3363 
3364  // Case A: non = 0. In this case nbgd < 0 hence we set nbgd = 0
3365  if (non == 0.0) {
3366  nbgd = 0;
3367  nonpred = ngam;
3368  logL = ngam;
3369  dlogLdsky = 1.0;
3370  d2logLdsky2 = 0.0;
3371  }
3372 
3373  // Case B: nbgd is positive
3374  else if (ngam < non * alpharat) {
3375  nbgd = non / alphap1 - ngam / alpha;
3376  nonpred = ngam + alpha * nbgd;
3377  logL = -ngam / alpha - non * std::log(alpharat);
3378  dlogLdsky = -1.0/alpha;
3379  d2logLdsky2 = 0.0;
3380  }
3381 
3382  // Case C: nbgd is zero or negative, hence set nbgd = 0
3383  else {
3384  nbgd = 0;
3385  nonpred = ngam;
3386  logL = ngam + non * (std::log(non) - std::log(ngam) - 1.0);
3387  dlogLdsky = 1.0 - non / ngam;
3388  d2logLdsky2 = non / (ngam * ngam);
3389  }
3390 
3391  } // endif: noff = 0
3392 
3393  // Special case non = 0
3394  else if (non == 0.0) {
3395  nbgd = noff / alphap1;
3396  nonpred = ngam + alpharat * noff;
3397  logL = ngam + noff * std::log(alphap1);
3398  dlogLdsky = 1.0;
3399  d2logLdsky2 = 0.0;
3400  } // endif: non = 0
3401 
3402  // General case
3403  else {
3404 
3405  // Compute log-likelihood value
3406  double alphat2 = 2.0 * alpha;
3407  double C = alpha * (non + noff) - alphap1 * ngam;
3408  double D = std::sqrt(C*C + 4.0 * alpha * alphap1 * noff * ngam);
3409  nbgd = (C + D) / (alphat2 * alphap1);
3410  nonpred = ngam + alpha * nbgd;
3411  logL = ngam + alphap1 * nbgd - non - noff -
3412  non * (std::log(nonpred) - std::log(non)) -
3413  noff * (std::log(nbgd) - std::log(noff));
3414 
3415  // Compute derivatives
3416  double f0 = alphat2 * noff - C;
3417  double dCds = -alphap1;
3418  double dDds = (C * dCds + 2.0 * alphap1 * alpha * noff) / D;
3419  double dbds = (f0 / D - 1.0) / alphat2;
3420  double d2bds2 = (-dCds / D - f0 / (D*D) * dDds) / alphat2;
3421  double f1 = alphap1 - alpha*non/nonpred - noff/nbgd;
3422  double f2 = nonpred * nonpred;
3423  double dpds = 1.0 + alpha * dbds;
3424  double f3 = non / f2 * dpds;
3425  dlogLdsky = 1.0 - non / nonpred + dbds * f1;
3426  d2logLdsky2 = f3 + d2bds2 * f1 +
3427  dbds * (alpha * f3 + noff / (nbgd*nbgd) * dbds);
3428 
3429  } // endelse: general case
3430 
3431  // Compile option: Check for NaN/Inf
3432  #if defined(G_NAN_CHECK)
3433  if (gammalib::is_notanumber(logL) || gammalib::is_infinite(logL) ||
3434  gammalib::is_notanumber(nonpred) || gammalib::is_infinite(nonpred) ||
3436  std::cout << "*** ERROR: GCTAOnOffObservation::wstat_value";
3437  std::cout << "(noff=" << noff;
3438  std::cout << ", alpha=" << alpha;
3439  std::cout << ", ngam=" << ngam << "):";
3440  std::cout << " NaN/Inf encountered";
3441  std::cout << " (logL=" << logL;
3442  std::cout << ", nonpred=" << nonpred;
3443  std::cout << ", nbgd=" << nbgd;
3444  std::cout << ", dlogLdsky=" << dlogLdsky;
3445  std::cout << ", d2logLdsky2=" << d2logLdsky2;
3446  std::cout << " )" << std::endl;
3447  }
3448  #endif
3449 
3450  // Return
3451  return logL;
3452 }
double m_livetime
Livetime of On observation (seconds)
virtual double likelihood_wstat(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for On/Off analysis in the case of Poisson signal with measured Pois...
std::string print(const GChatter &chatter=NORMAL) const
Print Redistribution Matrix File.
Definition: GRmf.cpp:866
const std::string & statistic(void) const
Return optimizer statistic.
double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Evaluate sum of all models.
Definition: GModels.cpp:911
Abstract model class.
Definition: GModel.hpp:100
#define G_ARF_RADIUS_CUT
std::string m_instrument
Instrument name.
void apply_ebounds(const GCTAObservation &obs)
Apply CTA observation energy boundaries to On/Off observation.
const double & factor_gradient(void) const
Return parameter factor gradient.
void init_members(void)
Initialise class members.
void append(const std::string &name, const std::vector< double > &column)
Append additional column to spectrum.
Definition: GPha.cpp:574
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
virtual bool apply_edisp(void) const
Signal if energy dispersion should be applied.
std::string m_name
Observation name.
virtual const GCTAResponse * response(void) const
Return pointer to CTA response function.
virtual ~GCTAOnOffObservation(void)
Destructor.
CTA event atom class definition.
const std::string & name(void) const
Return parameter name.
virtual GCTAOnOffObservation * clone(void) const
Clone instance.
void emin_obs(const GEnergy &emin_obs)
Set minimum energy of observations.
Definition: GPha.hpp:369
Abstract spectral model base class.
Sparse matrix class interface definition.
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:311
virtual const GCTAInstDir & dir(void) const
Return instrument direction of event bin.
Model container class definition.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
virtual void clear(void)
Clear instance.
GCTAEventBin class interface definition.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GPha.hpp:282
Sky regions container class definition.
GEbounds ebounds(void) const
Get energy boundaries.
virtual double livetime(void) const
Return livetime.
const GCTAOnOffObservation g_onoff_obs_hess_seed(true,"HESSOnOff")
bool contains(const GEnergy &eng) const
Checks whether energy boundaries contain energy.
Definition: GEbounds.cpp:1180
#define G_CONSTRUCTOR2
const GCTAOnOffObservation g_onoff_obs_astri_seed(true,"ASTRIOnOff")
CTA instrument response function class definition.
void compute_bgd(const GCTAObservation &obs, const GSkyRegions &off, const GModels &models, const bool &use_model_bkg)
Compute background rate in Off regions.
CTA event list class.
XML element node class.
Definition: GXmlElement.hpp:48
virtual void response(const GResponse &rsp)
Set response function.
Redistribution Matrix File class.
Definition: GRmf.hpp:55
GPha model_gamma(const GModels &models) const
const double & zenith(void) const
Return pointing zenith angle.
const std::string & rspname(void) const
Return response name.
const GCTAOnOffObservation g_onoff_obs_magic_seed(true,"MAGICOnOff")
std::string m_id
Observation identifier.
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
Abstract spatial model base class interface definition.
void backscal(const int &index, const double &backscal)
Set background scaling factor.
Definition: GPha.cpp:658
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
const GCTAOnOffObservation g_onoff_obs_fact_seed(true,"FACTOnOff")
Time class.
Definition: GTime.hpp:55
int nmeasured(void) const
Return number of measured energy bins in redistribution matrix.
Definition: GRmf.hpp:207
Gammalib tools definition.
const GPha & off_spec(void) const
Return Off spectrum.
virtual double wstat_value(const double &non, const double &noff, const double &alpha, const double &ngam, double &nonpred, double &nbgd, double &dlogLdsky, double &d2logLdsky2) const
Evaluate log-likelihood value in energy bin for On/Off analysis by profiling over the number of backg...
GEnergy emean(const int &index) const
Returns mean energy for a given energy interval.
Definition: GEbounds.cpp:1095
const GCTAOnOffObservation g_onoff_obs_veritas_seed(true,"VERITASOnOff")
const GArf & arf(void) const
Return Auxiliary Response File.
Interface for the circular sky region class.
virtual GCTAResponse * clone(void) const =0
Clones object.
void fill(const GEnergy &energy, const double &value=1.0)
Fill spectrum with a value.
Definition: GPha.cpp:771
GCTAResponse * m_response
Pointer to IRFs.
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
virtual double deadc(const GTime &time=GTime()) const
Return deadtime correction factor.
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1946
#define G_LIKELIHOOD
GEnergy ewidth(const int &index) const
Returns energy interval width.
Definition: GEbounds.cpp:1154
Source class definition.
void dir(const GSkyDir &dir)
Set sky direction.
bool is_free(void) const
Signal if parameter is free.
virtual double irf_spatial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response integrated over the spatial model.
Definition: GResponse.cpp:596
Auxiliary Response File class.
Definition: GArf.hpp:54
GArf m_arf
Auxiliary Response Function vector.
void id(const std::string &id)
Set observation identifier.
void load(const GFilename &filename)
Load Pulse Height Analyzer spectrum.
Definition: GPha.cpp:806
virtual void add_to_column(const int &column, const GVector &vector)
Add vector column into matrix.
#define G_WRITE
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & radius(void) const
Return circular region radius (in degrees)
CTA cube background class definition.
void map(const GSkyMap &map)
Set sky map.
void compute_alpha(const GCTAObservation &obs_on, const GCTAObservation &obs_off, const GSkyRegions &on, const GSkyRegions &off, const GModels &models, const bool &use_model_bkg)
Compute vector of alpha parameters.
const GEbounds & emeasured(void) const
Return measured energy boundaries.
Definition: GRmf.hpp:265
const GFilename & filename(void) const
Return file name.
Definition: GRmf.hpp:294
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
virtual double ontime(void) const
Return ontime.
const GFitsHeader & header(void) const
Return FITS header.
Definition: GPha.hpp:558
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
const double minerr
Minimum statistical error.
Implements FITS header card interface.
bool is_valid(const std::string &instruments, const std::string &ids) const
Verifies if model is valid for a given instrument and identifier.
Definition: GModel.cpp:570
Definition of support function used by CTA classes.
#define G_SET
Model parameter class.
Definition: GModelPar.hpp:87
double N_gamma(const GModels &models, const int &ibin, GVector *grad) const
const GRmf & rmf(void) const
Return Redistribution Matrix File.
void set(const GCTAObservation &obs_on, const GCTAObservation &obs_off, const GModels &models, const std::string &srcname, const GSkyRegions &on, const GSkyRegions &off, const bool &use_model_bkg)
Set On/Off observation from a CTA observation.
Interface for FITS header class.
Definition: GFitsHeader.hpp:49
const GFilename & filename(void) const
Return file name.
Definition: GPha.hpp:546
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
Abstract temporal model base class interface definition.
Model container class.
Definition: GModels.hpp:152
virtual double livetime(void) const
Return livetime.
const std::vector< int > & nonzero_indices(void) const
Get non-zero index vector.
CTA 2D effective area class.
Definition: GCTAAeff2D.hpp:55
virtual void read(const GXmlElement &xml)
Read On/Off observation from an XML element.
const std::string & id(void) const
Return observation identifier.
GCTAOnOffObservation & operator=(const GCTAOnOffObservation &obs)
Assignment operator.
void set_exposure(void)
Set ontime, livetime and deadtime correction factor.
CTA On/Off observation class definition.
Energy boundaries container class.
Definition: GEbounds.hpp:60
GFitsHeaderCard & append(const GFitsHeaderCard &card)
Append or update header card.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void copy_members(const GCTAOnOffObservation &obs)
Copy class members.
std::string classname(void) const
Return class name.
Definition: GModels.hpp:217
const GCTAResponse * cta_rsp(const std::string &origin, const GResponse &rsp)
Retrieve CTA response from generic observation.
bool has_par(const std::string &name) const
Checks if parameter name exists.
virtual void write(GXmlElement &xml) const
write observation to an xml element
const GCTAAeff * aeff(void) const
Return pointer to effective area response.
virtual double model(const GModels &models, const GEvent &event, GVector *gradients=NULL) const
Return model value and (optionally) gradients.
const GCTAOnOffObservation g_onoff_obs_cta_seed(true,"CTAOnOff")
#define G_RESPONSE_GET
std::string print(const GChatter &chatter=NORMAL) const
Print Auxiliary Response File.
Definition: GArf.cpp:834
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
Abstract data model base class interface definition.
CTA 2D effective area class definition.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
virtual std::string instrument(void) const
Return instrument name.
Interface definition for the observation registry class.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
virtual double likelihood_cstat(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for On/Off analysis in the case of Poisson signal with modeled Poiss...
Observation registry class definition.
CTA pointing class.
Abstract data model class.
Definition: GModelData.hpp:55
virtual std::string instrument(void) const
Return instrument.
int ntrue(void) const
Return number of true energy bins in redistribution matrix.
Definition: GRmf.hpp:193
const GModelSpectral * cta_model_spectral(const GModelData &model)
Retrieve spectral component from CTA background model.
GChatter
Definition: GTypemaps.hpp:33
Interface for a sky region map.
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1126
Optimizer parameters base class definition.
int size(void) const
Return number of regions in container.
const double & rad_max(void) const
Return radius cut value.
Definition: GCTAAeff2D.hpp:164
const GPha & on_spec(void) const
Return On spectrum.
Abstract observation base class.
int size(void) const
Return number of observations in container.
void clear(void)
Clear object.
Definition: GRmf.cpp:315
bool has_scales(void) const
Signals that model has scales.
Definition: GModel.hpp:355
GRmf m_rmf
Redistribution matrix.
CTA observation class interface definition.
#define G_COMPUTE_ARF
GPha model_background(const GModels &models) const
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:369
void clear(void)
Clear object.
Definition: GPha.cpp:447
Observation container class.
const double minmod
Minimum model value.
CTA instrument response function class.
double m_deadc
Deadtime correction (livetime/ontime)
CTA instrument response function class.
void exposure(const double &exposure)
Set exposure time.
Definition: GPha.hpp:340
Sky region container class.
Definition: GSkyRegions.hpp:56
Abstract spectral model base class interface definition.
void compute_arf_cut(const GCTAObservation &obs, const GModelSpatial &spatial, const GSkyRegions &on)
Compute ARF of On/Off observation for a IRF with radius cut.
const double & azimuth(void) const
Return pointing azimuth angle.
const GFitsHeader & header(void) const
Return FITS header.
Definition: GRmf.hpp:306
void name(const std::string &name)
Set observation name.
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1151
void compute_arf(const GCTAObservation &obs, const GModelSpatial &spatial, const GSkyRegions &on)
Compute ARF of On/Off observation.
#define G_RESPONSE_SET
void free_members(void)
Delete class members.
#define G_CONSTRUCTOR1
#define G_COMPUTE_RMF
Sky model class interface definition.
double N_bgd(const GModels &models, const int &ibin, GVector *grad) const
Sky model class.
Definition: GModelSky.hpp:122
virtual std::string print(const GChatter &chatter=NORMAL) const
Print On/Off observation information.
const GFitsHeader & header(void) const
Return FITS header.
Definition: GArf.hpp:266
std::string print(const GChatter &chatter=NORMAL) const
Print Pulse Height Analyzer spectrum.
Definition: GPha.cpp:1163
GCTAOnOffObservation(void)
Void constructor.
virtual double ontime(void) const
Return ontime.
Sky region map class interface definition.
int size(void) const
Return number of bins in spectrum.
Definition: GPha.hpp:254
double value(void) const
Return parameter value.
CTA event atom class.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
Sparse matrix class definition.
CTA event bin container class interface definition.
GPha m_off_spec
Off counts spectrum.
void check_consistency(const std::string &method) const
Check consistency of data members.
void clear(void)
Clear object.
Definition: GArf.cpp:395
virtual GObservation & operator=(const GObservation &obs)
Assignment operator.
virtual GEvents * events(void)
Return events.
virtual double likelihood(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for On/Off analysis.
virtual double prob_erecobin(const GEnergy &ereco_min, const GEnergy &ereco_max, const GEnergy &etrue, const double &theta) const =0
GRmf rmf_stacked(const GRmf &rmf, const GEnergy &emin, const GEnergy &emax) const
Return RMF for stacking.
void compute_rmf(const GCTAObservation &obs, const GSkyRegions &on)
Compute RMF of On/Off observation.
int npars(void) const
Return number of model parameters in container.
Definition: GModels.cpp:853
virtual int size(void) const
Return number of events in list.
GPha m_on_spec
On counts spectrum.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
const GFilename & filename(void) const
Return file name.
Definition: GArf.hpp:254
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const GSkyDir & dir(void) const
Return pointing sky direction.
std::string xml_get_attr(const std::string &origin, const GXmlElement &xml, const std::string &name, const std::string &attribute)
Return attribute value for a given parameter in XML element.
Definition: GTools.cpp:1738
Abstract spatial model base class.
#define G_MODEL_BACKGROUND
double m_ontime
Ontime of On observation (seconds)
Vector class.
Definition: GVector.hpp:46
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
virtual double flux(const GEnergy &emin, const GEnergy &emax) const =0
const GEbounds & etrue(void) const
Return true energy boundaries.
Definition: GRmf.hpp:251
CTA observation class.
Abstract instrument response base class.
Definition: GResponse.hpp:77
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
const GCTAEdisp * edisp(void) const
Return pointer to energy dispersion.
GArf arf_stacked(const GArf &arf, const GEnergy &emin, const GEnergy &emax) const
Return ARF for stacking.
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
Pulse Height Analyzer class.
Definition: GPha.hpp:64
Integration class interface definition.
#define G_READ
Observations container class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
int size(void) const
Return number of spectral bins.
Definition: GArf.hpp:209
void load(const GFilename &filename)
Load Auxiliary Response File.
Definition: GArf.cpp:552
bool contains(const GSkyDir &dir) const
Check if direction is contained in one of the regions.
Circular sky region class interface definition.
int size(void) const
Return number of models in container.
Definition: GModels.hpp:259
void caldb(const GCaldb &caldb)
Set calibration database.
GModel * append(const GModel &model)
Append model to container.
Definition: GModels.cpp:420
double arf_rad_max(const GCTAObservation &obs, const GSkyRegions &on) const
Check if effective area IRF has a radius cut.
const GEnergy & energy(void) const
Return energy.
const GCTAInstDir & dir(void) const
Return instrument direction.
const GCTAResponseIrf & cta_rsp_irf(const std::string &origin, const GObservation &obs)
Retrieve CTA IRF response from generic observation.
void load(const GFilename &filename)
Load Redistribution Matrix File.
Definition: GRmf.cpp:510
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GArf.hpp:237
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
void emax_obs(const GEnergy &emax_obs)
Set maximum energy of observations.
Definition: GPha.hpp:398
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1889
CTA On/Off observation class.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
GSkyRegion * append(const GSkyRegion &region)
Append region to container.