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