GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAOnOffObservation.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAOnOffObservation.hpp - CTA On/Off observation class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2018 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.hpp
23  * @brief CTA On/Off observation class definition
24  * @author Chia-Chun Lu & Christoph Deil & Pierrick Martin
25  */
26 
27 #ifndef GCTAONOFFOBSERVATION_HPP
28 #define GCTAONOFFOBSERVATION_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <string>
32 #include <vector>
33 #include "GPha.hpp"
34 #include "GArf.hpp"
35 #include "GRmf.hpp"
36 #include "GFunction.hpp"
37 #include "GObservation.hpp"
38 #include "GNodeArray.hpp"
39 
40 /* __ Forward declarations _______________________________________________ */
41 class GModels;
42 class GModelSpatial;
43 class GOptimizerPars;
44 class GObservations;
45 class GSkyRegions;
46 class GCTAObservation;
47 class GCTAResponse;
48 class GCTAResponseIrf;
49 
50 
51 /***********************************************************************//**
52  * @class GCTAOnOffObservation
53  *
54  * @brief CTA On/Off observation class
55  *
56  * This class defines a CTA On/Off observation. An On/Off observation is
57  * defined by two spectra, one for an On region including the source of
58  * interest, and one for an Off region including only background. The
59  * response of an On/Off observation is given by the Auxiliary Response File
60  * (ARF) and the Redistribution Matrix File (RMF).
61  *
62  * The class uses GPha objects to store the On and Off spectra, and GArf and
63  * GRmf objects to store the response information. On and Off regions are
64  * defined via GSkyRegionMap objects which are essentially finely-pixellized
65  * skymaps that allow handling arbitrarily complex shapes
66  ***************************************************************************/
68 
69 public:
70  // Constructors and destructors
72  GCTAOnOffObservation(const bool& dummy, const std::string& instrument);
74  GCTAOnOffObservation(const GPha& pha_on,
75  const GPha& pha_off,
76  const GArf& arf,
77  const GRmf& rmf);
79  const GModels& models,
80  const std::string& srcname,
81  const GEbounds& etrue,
82  const GEbounds& ereco,
83  const GSkyRegions& on,
84  const GSkyRegions& off,
85  const bool& use_model_bkg = true);
87  virtual ~GCTAOnOffObservation(void);
88 
89  // Operators
91 
92  // Implemented pure virtual methods
93  virtual void clear(void);
94  virtual GCTAOnOffObservation* clone(void) const;
95  virtual std::string classname(void) const;
96  virtual void response(const GResponse& rsp);
97  virtual const GCTAResponse* response(void) const;
98  virtual std::string instrument(void) const;
99  virtual double ontime(void) const;
100  virtual double livetime(void) const;
101  virtual double deadc(const GTime& time = GTime()) const;
102  virtual void read(const GXmlElement& xml);
103  virtual void write(GXmlElement& xml) const;
104  virtual std::string print(const GChatter& chatter = NORMAL) const;
105 
106  // Overloaded virtual methods
107  virtual double likelihood(const GModels& models,
108  GVector* gradient,
109  GMatrixSparse* curvature,
110  double* npred) const;
111  virtual int nobserved(void) const;
112 
113  // Other methods
114  void instrument(const std::string& instrument);
115  bool has_response(void) const;
116  const GPha& on_spec(void) const;
117  const GPha& off_spec(void) const;
118  const GArf& arf(void) const;
119  const GRmf& rmf(void) const;
120  GPha model_gamma(const GModels& models) const;
121  GPha model_background(const GModels& models) const;
122 
123 protected:
124  // Protected methods
125  void init_members(void);
126  void copy_members(const GCTAOnOffObservation& obs);
127  void free_members(void);
128  void set_exposure(void);
129  void check_consistency(const std::string& method) const;
130  GArf arf_stacked(const GArf& arf, const GEnergy& emin, const GEnergy& emax) const;
131  GRmf rmf_stacked(const GRmf& rmf, const GEnergy& emin, const GEnergy& emax) const;
132  void set(const GCTAObservation& obs,
133  const GModels& models,
134  const std::string& srcname,
135  const GSkyRegions& on,
136  const GSkyRegions& off,
137  const bool& use_model_bkg);
138  void compute_arf(const GCTAObservation& obs,
139  const GModelSpatial& spatial,
140  const GSkyRegions& on);
141  void compute_arf_cut(const GCTAObservation& obs,
142  const GModelSpatial& spatial,
143  const GSkyRegions& on);
144  void compute_bgd(const GCTAObservation& obs,
145  const GSkyRegions& off,
146  const GModels& models,
147  const bool& use_model_bkg);
148  void compute_alpha(const GCTAObservation& obs,
149  const GSkyRegions& on,
150  const GSkyRegions& off,
151  const GModels& models,
152  const bool& use_model_bkg);
153  void compute_rmf(const GCTAObservation& obs,
154  const GSkyRegions& on);
155  void apply_ebounds(const GCTAObservation& obs);
156  double arf_rad_max(const GCTAObservation& obs, const GSkyRegions& on) const;
157  double N_gamma(const GModels& models, const int& ibin, GVector* grad) const;
158  double N_bgd(const GModels& models, const int& ibin, GVector* grad) const;
159 
160  // Likelihood methods
161  virtual double likelihood_cstat(const GModels& models,
162  GVector* gradient,
163  GMatrixSparse* curvature,
164  double* npred) const;
165  virtual double likelihood_wstat(const GModels& models,
166  GVector* gradient,
167  GMatrixSparse* curvature,
168  double* npred) const;
169  virtual double wstat_value(const double& non,
170  const double& noff,
171  const double& alpha,
172  const double& ngam,
173  double& nonpred,
174  double& nbgd,
175  double& dlogLdsky,
176  double& d2logLdsky2) const;
177 
178  // Protected data members
179  std::string m_instrument; //!< Instrument name
180  GCTAResponse* m_response; //!< Pointer to IRFs
181  double m_ontime; //!< Ontime (seconds)
182  double m_livetime; //!< Livetime (seconds)
183  double m_deadc; //!< Deadtime correction (livetime/ontime)
184  GPha m_on_spec; //!< On counts spectrum
185  GPha m_off_spec; //!< Off counts spectrum
186  GArf m_arf; //!< Auxiliary Response Function vector
187  GRmf m_rmf; //!< Redistribution matrix
188 };
189 
190 
191 /***********************************************************************//**
192  * @brief Return class name
193  *
194  * @return String containing the class name ("GCTAOnOffObservation").
195  ***************************************************************************/
196 inline
197 std::string GCTAOnOffObservation::classname(void) const
198 {
199  return ("GCTAOnOffObservation");
200 }
201 
202 
203 /***********************************************************************//**
204  * @brief Return instrument
205  *
206  * @return Instrument.
207  ***************************************************************************/
208 inline
209 std::string GCTAOnOffObservation::instrument(void) const
210 {
211  return (m_instrument);
212 }
213 
214 
215 /***********************************************************************//**
216  * @brief Return ontime
217  *
218  * @return Ontime in seconds.
219  ***************************************************************************/
220 inline
222 {
223  return (m_ontime);
224 }
225 
226 
227 /***********************************************************************//**
228  * @brief Return livetime
229  *
230  * @return Livetime in seconds.
231  ***************************************************************************/
232 inline
234 {
235  return (m_livetime);
236 }
237 
238 
239 /***********************************************************************//**
240  * @brief Return deadtime correction factor
241  *
242  * @param[in] time Time.
243  * @return Deadtime correction factor.
244  *
245  * Returns the deadtime correction factor. Optionally, this method takes a
246  * @p time argument that takes provision for returning the deadtime
247  * correction factor as function of time.
248  *
249  * The deadtime correction factor is defined as the livetime divided by the
250  * ontime.
251  ***************************************************************************/
252 inline
253 double GCTAOnOffObservation::deadc(const GTime& time) const
254 {
255  return (m_deadc);
256 }
257 
258 
259 /***********************************************************************//**
260  * @brief Signal if observation contains response information
261  *
262  * @return True if observation contains response information.
263  ***************************************************************************/
264 inline
266 {
267  return (m_response != NULL);
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Set instrument
273  *
274  * @param[in] instrument Instrument.
275  ***************************************************************************/
276 inline
277 void GCTAOnOffObservation::instrument(const std::string& instrument)
278 {
280  return;
281 }
282 
283 
284 /***********************************************************************//**
285  * @brief Return On spectrum
286  *
287  * @return On spectrum.
288  ***************************************************************************/
289 inline
291 {
292  return (m_on_spec);
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Return Off spectrum
298  *
299  * @return Off spectrum.
300  ***************************************************************************/
301 inline
303 {
304  return (m_off_spec);
305 }
306 
307 
308 /***********************************************************************//**
309  * @brief Return Auxiliary Response File
310  *
311  * @return Auxiliary Response File.
312  ***************************************************************************/
313 inline
314 const GArf& GCTAOnOffObservation::arf(void) const
315 {
316  return (m_arf);
317 }
318 
319 
320 /***********************************************************************//**
321  * @brief Return Redistribution Matrix File
322  *
323  * @return Redistribution Matrix File.
324  ***************************************************************************/
325 inline
326 const GRmf& GCTAOnOffObservation::rmf(void) const
327 {
328  return (m_rmf);
329 }
330 
331 
332 /***********************************************************************//**
333  * @brief Return number of observed events
334  *
335  * @return Number of observed events.
336  ***************************************************************************/
337 inline
339 {
340  return (int(m_on_spec.counts()+0.5));
341 }
342 
343 #endif /* GCTAONOFFOBSERVATION_HPP */
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 m_instrument
Instrument name.
void apply_ebounds(const GCTAObservation &obs)
Apply CTA observation energy boundaries to On/Off observation.
void init_members(void)
Initialise class members.
virtual const GCTAResponse * response(void) const
Return pointer to CTA response function.
virtual ~GCTAOnOffObservation(void)
Destructor.
virtual GCTAOnOffObservation * clone(void) const
Clone instance.
Sparse matrix class interface definition.
virtual double npred(const GModels &models, GVector *gradient=NULL) const
Return total number (and optionally gradient) of predicted counts for all models. ...
virtual void clear(void)
Clear instance.
virtual double livetime(void) const
Return livetime.
double counts(void) const
Number of counts in spectrum.
Definition: GPha.cpp:703
Optimizer parameter container class.
void compute_bgd(const GCTAObservation &obs, const GSkyRegions &off, const GModels &models, const bool &use_model_bkg)
Compute background rate in Off regions.
XML element node class.
Definition: GXmlElement.hpp:47
Redistribution Matrix File class.
Definition: GRmf.hpp:55
GPha model_gamma(const GModels &models) const
Time class.
Definition: GTime.hpp:54
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...
const GArf & arf(void) const
Return Auxiliary Response File.
GCTAResponse * m_response
Pointer to IRFs.
Auxiliary Response File class.
Definition: GArf.hpp:54
GArf m_arf
Auxiliary Response Function vector.
double N_gamma(const GModels &models, const int &ibin, GVector *grad) const
XSPEC Pulse Height Analyzer class definition.
const GRmf & rmf(void) const
Return Redistribution Matrix File.
Node array class interface definition.
Model container class.
Definition: GModels.hpp:150
virtual void read(const GXmlElement &xml)
Read On/Off observation from an XML element.
GCTAOnOffObservation & operator=(const GCTAOnOffObservation &obs)
Assignment operator.
XSPEC Redistribution Matrix File class definition.
void set_exposure(void)
Set ontime, livetime and deadtime correction factor.
Energy boundaries container class.
Definition: GEbounds.hpp:60
Single parameter function abstract base class definition.
void copy_members(const GCTAOnOffObservation &obs)
Copy class members.
virtual void write(GXmlElement &xml) const
write observation to an xml element
XSPEC Auxiliary Response File 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.
bool has_response(void) const
Signal if observation contains response information.
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...
virtual std::string instrument(void) const
Return instrument.
GChatter
Definition: GTypemaps.hpp:33
const GPha & on_spec(void) const
Return On spectrum.
virtual double deadc(const GTime &time=GTime()) const
Return deadtime correction factor.
Abstract observation base class.
GRmf m_rmf
Redistribution matrix.
Abstract observation base class interface definition.
GPha model_background(const GModels &models) const
Observation container class.
CTA instrument response function class.
double m_deadc
Deadtime correction (livetime/ontime)
CTA instrument response function class.
Sky region container class.
Definition: GSkyRegions.hpp:56
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.
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.
void free_members(void)
Delete class members.
double N_bgd(const GModels &models, const int &ibin, GVector *grad) const
virtual std::string print(const GChatter &chatter=NORMAL) const
Print On/Off observation information.
GCTAOnOffObservation(void)
Void constructor.
virtual double ontime(void) const
Return ontime.
GPha m_off_spec
Off counts spectrum.
void check_consistency(const std::string &method) const
Check consistency of data members.
virtual int nobserved(void) const
Return number of observed events.
virtual double likelihood(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for On/Off analysis.
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.
GPha m_on_spec
On counts spectrum.
Abstract spatial model base class.
double m_ontime
Ontime (seconds)
Vector class.
Definition: GVector.hpp:46
virtual std::string classname(void) const
Return class name.
CTA observation class.
Abstract instrument response base class.
Definition: GResponse.hpp:67
GArf arf_stacked(const GArf &arf, const GEnergy &emin, const GEnergy &emax) const
Return ARF for stacking.
Pulse Height Analyzer class.
Definition: GPha.hpp:64
double arf_rad_max(const GCTAObservation &obs, const GSkyRegions &on) const
Check if effective area IRF has a radius cut.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
CTA On/Off observation class.