Thierry Louge / gammalib

Commit f54cd3020f284c6d6262c9fa3e05ff72faaa04cc

Authored by Jürgen Knödlseder
1 parent 92d74e59

Adjust the interface for On/Off observations

Inline Side-by-side
Showing 4 changed files with 423 additions and 495 deletions
inst/cta/include/GCTAOnOffObservation.hpp
  View file @f54cd30
... ... @@ -40,20 +40,33 @@
40 40 class GModels;
41 41 class GOptimizerPars;
42 42 class GCTAObservation;
  43 +class GBounds;
43 44  
44 45  
45 46 /***********************************************************************//**
46 47 * @class GCTAOnOffObservation
47 48 *
48 49 * @brief CTA On/Off observation class
  50 + *
  51 + * This class defines a CTA On/Off observation. An On/Off observation is
  52 + * defined by two spectra, one for an On region including the source of
  53 + * interest, and one for an Off region including only background. The
  54 + * response of an On/Off observation is given by the Auxiliary Response File
  55 + * (ARF) and the Redistribution Matrix File (RMF).
  56 + *
  57 + * The class uses GPha objects to store the On and Off spectra, and GArf and
  58 + * GRmf objects to store the response information.
49 59 ***************************************************************************/
50 60 class GCTAOnOffObservation : public GObservation {
51 61  
52 62 public:
53 63 // Constructors and destructors
54 64 GCTAOnOffObservation(void);
55   - GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
56   - const GSkyRegions& off);
  65 + GCTAOnOffObservation(const GCTAObservation& obs,
  66 + const GEbounds& etrue,
  67 + const GEbounds& ereco,
  68 + const GSkyRegions& on,
  69 + const GSkyRegions& off);
57 70 GCTAOnOffObservation(const GCTAOnOffObservation& obs);
58 71 virtual ~GCTAOnOffObservation(void);
59 72  
... ... @@ -82,28 +95,28 @@ public:
82 95 virtual int nobserved(void) const;
83 96  
84 97 // Other methods
85   - void instrument(const std::string& instrument);
86   - void on_regions(const GSkyRegions& regions);
87   - void off_regions(const GSkyRegions& regions);
88   - const GPha& on_spec(void) const;
89   - const GPha& off_spec(void) const;
90   - const GArf& arf(void) const;
91   - const GRmf& rmf(void) const;
92   - void fill(const GCTAObservation& obs);
93   - void compute_response(const GCTAObservation& obs,
94   - const GEbounds& etrue);
  98 + void instrument(const std::string& instrument);
  99 + void on_regions(const GSkyRegions& regions);
  100 + void off_regions(const GSkyRegions& regions);
  101 + const GSkyRegions& on_regions(void) const;
  102 + const GSkyRegions& off_regions(void) const;
  103 + const GPha& on_spec(void) const;
  104 + const GPha& off_spec(void) const;
  105 + const GArf& arf(void) const;
  106 + const GRmf& rmf(void) const;
95 107  
96 108 protected:
97 109 // Protected methods
98 110 void init_members(void);
99 111 void copy_members(const GCTAOnOffObservation& obs);
100 112 void free_members(void);
101   - void compute_alpha(const GCTAObservation& obs);
  113 + void set(const GCTAObservation& obs);
102 114 void compute_arf(const GCTAObservation& obs);
103 115 void compute_bgd(const GCTAObservation& obs);
104   - void compute_rmf(const GCTAObservation& obs, const GEbounds& etrue);
105   - double model_on(const GModels& models, int ibin, GVector* mod_grad) const;
106   - double model_off(const GModels& models, int ibin, GVector* mod_grad) const;
  116 + void compute_alpha(const GCTAObservation& obs);
  117 + void compute_rmf(const GCTAObservation& obs);
  118 + double N_gamma(const GModels& models, const int& ibin, GVector* grad) const;
  119 + double N_bgd(const GModels& models, const int& ibin, GVector* grad) const;
107 120  
108 121 // Protected data members
109 122 std::string m_instrument; //!< Instrument name
... ... @@ -228,9 +241,33 @@ void GCTAOnOffObservation::off_regions(const GSkyRegions&amp; regions)
228 241  
229 242  
230 243 /***********************************************************************//**
231   - * @brief Return ON spectrum
  244 + * @brief Return On regions
  245 + *
  246 + * @return On regions.
  247 + ***************************************************************************/
  248 +inline
  249 +const GSkyRegions& GCTAOnOffObservation::on_regions(void) const
  250 +{
  251 + return (m_on_regions);
  252 +}
  253 +
  254 +
  255 +/***********************************************************************//**
  256 + * @brief Return Off regions
  257 + *
  258 + * @return Off regions.
  259 + ***************************************************************************/
  260 +inline
  261 +const GSkyRegions& GCTAOnOffObservation::off_regions(void) const
  262 +{
  263 + return (m_off_regions);
  264 +}
  265 +
  266 +
  267 +/***********************************************************************//**
  268 + * @brief Return On spectrum
232 269 *
233   - * @return ON spectrum.
  270 + * @return On spectrum.
234 271 ***************************************************************************/
235 272 inline
236 273 const GPha& GCTAOnOffObservation::on_spec(void) const
... ... @@ -240,9 +277,9 @@ const GPha&amp; GCTAOnOffObservation::on_spec(void) const
240 277  
241 278  
242 279 /***********************************************************************//**
243   - * @brief Return OFF spectrum
  280 + * @brief Return Off spectrum
244 281 *
245   - * @return OFF spectrum.
  282 + * @return Off spectrum.
246 283 ***************************************************************************/
247 284 inline
248 285 const GPha& GCTAOnOffObservation::off_spec(void) const
... ...
inst/cta/pyext/GCTAOnOffObservation.i
  View file @f54cd30
... ... @@ -63,8 +63,11 @@ class GCTAOnOffObservation : public GObservation {
63 63 public:
64 64 // Constructors and destructors
65 65 GCTAOnOffObservation(void);
66   - GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
67   - const GSkyRegions& off);
  66 + GCTAOnOffObservation(const GCTAObservation& obs,
  67 + const GEbounds& etrue,
  68 + const GEbounds& ereco,
  69 + const GSkyRegions& on,
  70 + const GSkyRegions& off);
68 71 GCTAOnOffObservation(const GCTAOnOffObservation& obs);
69 72 virtual ~GCTAOnOffObservation(void);
70 73  
... ... @@ -89,16 +92,15 @@ public:
89 92 virtual int nobserved(void) const;
90 93  
91 94 // Other methods
92   - void instrument(const std::string& instrument);
93   - void on_regions(const GSkyRegions& regions);
94   - void off_regions(const GSkyRegions& regions);
95   - const GPha& on_spec(void) const;
96   - const GPha& off_spec(void) const;
97   - const GArf& arf(void) const;
98   - const GRmf& rmf(void) const;
99   - void fill(const GCTAObservation& obs);
100   - void compute_response(const GCTAObservation& obs,
101   - const GEbounds& etrue);
  95 + void instrument(const std::string& instrument);
  96 + void on_regions(const GSkyRegions& regions);
  97 + void off_regions(const GSkyRegions& regions);
  98 + const GSkyRegions& on_regions(void) const;
  99 + const GSkyRegions& off_regions(void) const;
  100 + const GPha& on_spec(void) const;
  101 + const GPha& off_spec(void) const;
  102 + const GArf& arf(void) const;
  103 + const GRmf& rmf(void) const;
102 104 };
103 105  
104 106  
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @f54cd30
... ... @@ -55,16 +55,15 @@ const GObservationRegistry g_onoff_obs_cta_registry(&amp;g_onoff_obs_cta_seed);
55 55 #define G_RESPONSE_GET "GCTAOnOffObservation::response()"
56 56 #define G_WRITE "GCTAOnOffObservation::write(GXmlElement&)"
57 57 #define G_READ "GCTAOnOffObservation::read(GXmlElement&)"
58   -#define G_FILL "GCTAOnOffObservation::fill(GCTAObservation&)"
59   -#define G_COMPUTE_ALPHA "GCTAOnOffObservation::compute_alpha("\
60   - "GCTAObservation&)"
61   -#define G_COMPUTE_ARF "GCTAOnOffObservation::compute_arf(GCTAObservation&)"
62   -#define G_COMPUTE_BGD "GCTAOnOffObservation::compute_bgd(GCTAObservation&)"
63   -#define G_COMPUTE_RMF "GCTAOnOffObservation::compute_rmf(GCTAObservation&,"\
64   - " GEbounds&)"
65 58 #define G_LIKELIHOOD "GCTAOnOffObservation::likelihood(GModels&, "\
66 59 "GOptimizerPars&, GMatrixSparse&, "\
67 60 "GVector&, double&, double&)"
  61 +#define G_SET "GCTAOnOffObservation::set(GCTAObservation&)"
  62 +#define G_COMPUTE_ARF "GCTAOnOffObservation::compute_arf(GCTAObservation&)"
  63 +#define G_COMPUTE_BGD "GCTAOnOffObservation::compute_bgd(GCTAObservation&)"
  64 +#define G_COMPUTE_ALPHA "GCTAOnOffObservation::compute_alpha("\
  65 + "GCTAObservation&)"
  66 +#define G_COMPUTE_RMF "GCTAOnOffObservation::compute_rmf(GCTAObservation&)"
68 67  
69 68 /* __ Constants __________________________________________________________ */
70 69 const double minmod = 1.0e-100; //!< Minimum model value
... ... @@ -86,6 +85,8 @@ const double minerr = 1.0e-100; //!&lt; Minimum statistical error
86 85  
87 86 /***********************************************************************//**
88 87 * @brief Void constructor
  88 + *
  89 + * Constructs empty On/Off observation.
89 90 ***************************************************************************/
90 91 GCTAOnOffObservation::GCTAOnOffObservation(void)
91 92 {
... ... @@ -100,7 +101,7 @@ GCTAOnOffObservation::GCTAOnOffObservation(void)
100 101 /***********************************************************************//**
101 102 * @brief Copy constructor
102 103 *
103   - * @param[in] obs CTA On/Off observation.
  104 + * @param[in] obs On/Off observation.
104 105 ***************************************************************************/
105 106 GCTAOnOffObservation::GCTAOnOffObservation(const GCTAOnOffObservation& obs)
106 107 {
... ... @@ -116,15 +117,25 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GCTAOnOffObservation&amp; obs)
116 117  
117 118  
118 119 /***********************************************************************//**
119   - * @brief Parameter constructor
  120 + * @brief CTA observation constructor
120 121 *
121   - * @param[in] ereco Reconstructed energy bins.
  122 + * @param[in] obs CTA observation.
  123 + * @param[in] etrue True energy boundaries.
  124 + * @param[in] ereco Reconstructed energy boundaries.
122 125 * @param[in] on On regions.
123 126 * @param[in] off Off regions.
  127 + *
  128 + * Constructs On/Off observation a CTA observation by filling the On and Off
  129 + * spectra and computing the Auxiliary Response File (ARF) and Redistribution
  130 + * Matrix File (RMF). The method requires the specification of the true and
  131 + * reconstructed energy boundaries, as well as the definition of On and Off
  132 + * regions.
124 133 ***************************************************************************/
125   -GCTAOnOffObservation::GCTAOnOffObservation(const GEbounds& ereco,
126   - const GSkyRegions& on,
127   - const GSkyRegions& off)
  134 +GCTAOnOffObservation::GCTAOnOffObservation(const GCTAObservation& obs,
  135 + const GEbounds& etrue,
  136 + const GEbounds& ereco,
  137 + const GSkyRegions& on,
  138 + const GSkyRegions& off)
128 139 {
129 140 // Initialise private
130 141 init_members();
... ... @@ -133,10 +144,17 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GEbounds&amp; ereco,
133 144 m_on_spec = GPha(ereco);
134 145 m_off_spec = GPha(ereco);
135 146  
  147 + // Initialise response information
  148 + m_arf = GArf(ereco);
  149 + m_rmf = GRmf(etrue, ereco);
  150 +
136 151 // Store regions
137 152 m_on_regions = on;
138 153 m_off_regions = off;
139 154  
  155 + // Set On/Off observation from CTA observation
  156 + set(obs);
  157 +
140 158 // Return
141 159 return;
142 160 }
... ... @@ -164,10 +182,10 @@ GCTAOnOffObservation::~GCTAOnOffObservation(void)
164 182 /***********************************************************************//**
165 183 * @brief Assignment operator
166 184 *
167   - * @param[in] obs CTA On/Off observation.
168   - * @return CTA On/Off observation.
  185 + * @param[in] obs On/Off observation.
  186 + * @return On/Off observation.
169 187 *
170   - * Assigns one CTA On/Off observation to another On/Off observation object.
  188 + * Assigns one On/Off observation to another On/Off observation object.
171 189 ***************************************************************************/
172 190 GCTAOnOffObservation& GCTAOnOffObservation::operator=(const GCTAOnOffObservation& obs)
173 191 {
... ... @@ -237,7 +255,7 @@ GCTAOnOffObservation* GCTAOnOffObservation::clone(void) const
237 255 * @exception GException::invalid_argument
238 256 * Invalid response class specified.
239 257 *
240   - * Sets the response function for the observation.
  258 + * Sets the response function for the On/Off observation.
241 259 ***************************************************************************/
242 260 void GCTAOnOffObservation::response(const GResponse& rsp)
243 261 {
... ... @@ -291,17 +309,12 @@ const GCTAResponse* GCTAOnOffObservation::response(void) const
291 309  
292 310  
293 311 /***********************************************************************//**
294   - * @brief Read ON/OFFobservation from an xml element
  312 + * @brief Read On/Off observation from an XML element
295 313 *
296 314 * @param[in] xml XML element.
297 315 *
298   - * @exception GException::xml_invalid_parnum
299   - * Invalid number of parameters found in XML element.
300   - * @exception GException::xml_invalid_parnames
301   - * Invalid parameter names found in XML element.
302   - *
303   - * Reads information for a CTA On/Off observation from an XML element.
304   - * The expected format of the XML element is
  316 + * Reads information for an On/Off observation from an XML element. The
  317 + * expected format of the XML element is
305 318 *
306 319 * <observation name="..." id="..." instrument="...">
307 320 * <parameter name="Pha_on" file="..."/>
... ... @@ -321,109 +334,21 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
321 334 // Extract instrument name
322 335 m_instrument = xml.attribute("instrument");
323 336  
324   - // Determine number of parameter nodes in XML element
325   - int npars = xml.elements("parameter");
326   -
327   - // Verify that XML element has exactly 6 parameters
328   - if (xml.elements() != 6 || npars != 6) {
329   - throw GException::xml_invalid_parnum(G_READ, xml,
330   - "CTA observation requires exactly 6 parameters.");
331   - }
332   -
333   - // Extract parameters
334   - int npar[] = {0, 0, 0, 0, 0, 0};
335   - for (int i = 0; i < npars; ++i) {
336   -
337   - // Get parameter element
338   - const GXmlElement* par = xml.element("parameter", i);
339   -
340   - // Handle Pha_on
341   - if (par->attribute("name") == "Pha_on") {
342   -
343   - // Read pha_on file name
344   - std::string filename = par->attribute("file");
345   -
346   - // Load pha_on
347   - m_on_spec.load(filename);
348   -
349   - // Increment parameter counter
350   - npar[0]++;
351   - }
352   -
353   - // Handle Pha_off
354   - else if (par->attribute("name") == "Pha_off") {
355   -
356   - // Read countsmap file name
357   - std::string filename = par->attribute("file");
358   -
359   - // Load pha_off
360   - m_off_spec.load(filename);
361   -
362   - // Increment parameter counter
363   - npar[1]++;
364   - }
365   -
366   - // Handle on regions
367   - else if(par->attribute("name") == "Regions_on") {
368   -
369   - // Get filename
370   - std::string filename = par->attribute("file");
371   -
372   - // Load on regions
373   - m_on_regions.load(filename);
374   -
375   - // Increase number of parameters
376   - npar[2]++;
377   - }
378   -
379   - // Handle off regions
380   - else if(par->attribute("name") == "Regions_off") {
381   -
382   - // Get filename
383   - std::string filename = par->attribute("file");
384   -
385   - // load off regions
386   - m_on_regions.load(filename);
387   -
388   - // Increase number of parameters
389   - npar[3]++;
390   - }
391   -
392   - // Handle Arf
393   - else if(par->attribute("name") == "Arf") {
394   -
395   - // Get filename
396   - std::string filename = par->attribute("file");
397   -
398   - // Load arf
399   - m_arf.load(filename);
400   -
401   - // Increase number of parameters
402   - npar[4]++;
403   - }
404   -
405   - // Handle Rmf
406   - else if(par->attribute("name") == "Rmf") {
407   -
408   - // Get filename
409   - std::string filename = par->attribute("file");
410   -
411   - // Load Rmf
412   - m_rmf.load(filename);
413   -
414   - // Increase number of parameters
415   - npar[5]++;
416   - }
417   -
418   - } // endfor: looped over all parameters
419   -
420   - // Verify that all parameters were found
421   - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 ||
422   - npar[3] != 1 || npar[4] != 1 || npar[5] != 1) {
423   - throw GException::xml_invalid_parnames(G_READ, xml,
424   - "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\""
425   - ", \"Regions_off\",\"Arf\" and \"Rmf\" parameters.");
426   - }
  337 + // Get parameters
  338 + std::string pha_on = gammalib::xml_get_attr(G_READ, xml, "Pha_on", "file");
  339 + std::string pha_off = gammalib::xml_get_attr(G_READ, xml, "Pha_off", "file");
  340 + std::string reg_on = gammalib::xml_get_attr(G_READ, xml, "Regions_on", "file");
  341 + std::string reg_off = gammalib::xml_get_attr(G_READ, xml, "Regions_off", "file");
  342 + std::string arf = gammalib::xml_get_attr(G_READ, xml, "Arf", "file");
  343 + std::string rmf = gammalib::xml_get_attr(G_READ, xml, "Rmf", "file");
  344 +
  345 + // Load files
  346 + m_on_spec.load(pha_on);
  347 + m_off_spec.load(pha_off);
  348 + m_on_regions.load(reg_on);
  349 + m_off_regions.load(reg_off);
  350 + m_arf.load(arf);
  351 + m_rmf.load(rmf);
427 352  
428 353 // Return
429 354 return;
... ... @@ -435,12 +360,7 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
435 360 *
436 361 * @param[in] xml XML element.
437 362 *
438   - * @exception GException::xml_invalid_parnum
439   - * Invalid number of parameters found in XML element.
440   - * @exception GException::xml_invalid_parnames
441   - * Invalid parameter names found in XML element.
442   - *
443   - * Writes information for a CTA On/Off observation into an XML element. The
  363 + * Writes information for an On/Off observation into an XML element. The
444 364 * expected format of the XML element is
445 365 *
446 366 * <observation name="..." id="..." instrument="...">
... ... @@ -452,150 +372,36 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
452 372 * <parameter name="Rmf" file="..."/>
453 373 * </observation>
454 374 *
  375 + * The actual files described in the XML elements are not written.
455 376 ***************************************************************************/
456 377 void GCTAOnOffObservation::write(GXmlElement& xml) const
457 378 {
458   - // If XML element has 0 nodes then append 6 parameter nodes
459   - if (xml.elements() == 0) {
460   - xml.append(GXmlElement("parameter name=\"Pha_on\""));
461   - xml.append(GXmlElement("parameter name=\"Pha_off\""));
462   - xml.append(GXmlElement("parameter name=\"Regions_on\""));
463   - xml.append(GXmlElement("parameter name=\"Regions_off\""));
464   - xml.append(GXmlElement("parameter name=\"Arf\""));
465   - xml.append(GXmlElement("parameter name=\"Rmf\""));
466   - }
467   -
468   - // Verify that XML element has exactly 4 parameters
469   - if (xml.elements() != 6 || xml.elements("parameter") != 6) {
470   - throw GException::xml_invalid_parnum(G_WRITE, xml,
471   - "CTAOnOffObservation requires exactly 6 parameters.");
472   - }
  379 + // Allocate XML element pointer
  380 + GXmlElement* par;
473 381  
474   - // Set or update parameter attributes
475   - int npar[] = {0, 0, 0, 0, 0, 0};
476   - for (int i = 0; i < 6; ++i) {
  382 + // Set Pha_on parameter
  383 + par = gammalib::xml_need_par(G_WRITE, xml, "Pha_on");
  384 + par->attribute("file", m_on_spec.filename());
477 385  
478   - // Get parameter element
479   - GXmlElement* par = xml.element("parameter", i);
480   -
481   - // Handle on counts
482   - if (par->attribute("name") == "Pha_on") {
483   - par->attribute("file", m_on_spec.filename());
484   - npar[0]++;
485   - }
486   -
487   - // Handle off counts
488   - else if (par->attribute("name") == "Pha_off") {
489   - par->attribute("file", m_off_spec.filename());
490   - npar[1]++;
491   - }
492   -
493   - // Handle on regions
494   - else if (par->attribute("name") == "Regions_on") {
495   - par->attribute("file", m_on_regions.filename());
496   - npar[2]++;
497   - }
  386 + // Set Pha_off parameter
  387 + par = gammalib::xml_need_par(G_WRITE, xml, "Pha_off");
  388 + par->attribute("file", m_off_spec.filename());
498 389  
499   - // Handle off regions
500   - else if (par->attribute("name") == "Regions_off") {
501   - par->attribute("file", m_off_regions.filename());
502   - npar[3]++;
503   - }
  390 + // Set Regions_on parameter
  391 + par = gammalib::xml_need_par(G_WRITE, xml, "Regions_on");
  392 + par->attribute("file", m_on_regions.filename());
504 393  
505   - // Handle effective area
506   - else if (par->attribute("name") == "Arf") {
507   - par->attribute("file", m_arf.filename());
508   - npar[4]++;
509   - }
  394 + // Set Regions_off parameter
  395 + par = gammalib::xml_need_par(G_WRITE, xml, "Regions_off");
  396 + par->attribute("file", m_off_regions.filename());
510 397  
511   - // Handle energy resolution
512   - else if (par->attribute("name") == "Rmf") {
513   - par->attribute("file", m_rmf.filename());
514   - npar[5]++;
515   - }
516   - }
  398 + // Set Arf parameter
  399 + par = gammalib::xml_need_par(G_WRITE, xml, "Arf");
  400 + par->attribute("file", m_arf.filename());
517 401  
518   - // Verify that all required parameters are present
519   - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 ||
520   - npar[3] != 1 || npar[4] != 1 || npar[5] != 1) {
521   - throw GException::xml_invalid_parnames(G_WRITE, xml,
522   - "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\","
523   - " \"Regions_off\", \"Arf\", and \"Rmf\" parameters.");
524   - }
525   -
526   - // Return
527   - return;
528   -}
529   -
530   -
531   -/***********************************************************************//**
532   - * @brief Fill events in the On and Off spectra
533   - *
534   - * @param[in] obs CTA observation.
535   - *
536   - * @exception GException::invalid_value
537   - * No CTA event list found in CTA observation.
538   - *
539   - * Fills events from an observation into the On and Off spectra and cumulates
540   - * their livetime. This method adds to any existing content in the On and
541   - * Off spectra.
542   - ***************************************************************************/
543   -void GCTAOnOffObservation::fill(const GCTAObservation& obs)
544   -{
545   - // Get CTA event list pointer
546   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
547   - if (events == NULL) {
548   - std::string msg = "No event list found in CTA observation \""+
549   - obs.name()+"\" (ID="+obs.id()+"). ON/OFF observation "
550   - "can only be filled from event list.";
551   - throw GException::invalid_value(G_FILL, msg);
552   - }
553   -
554   - // Loop over all events
555   - for (int i = 0; i < events->size(); ++i) {
556   -
557   - // Get measured event direction
558   - const GCTAEventAtom* atom = (*events)[i];
559   - GSkyDir dir = atom->dir().dir();
560   -
561   - // Fill in spectrum according to region containment
562   - if (m_on_regions.contains(dir)) {
563   - m_on_spec.fill(atom->energy());
564   - }
565   - if (m_off_regions.contains(dir)) {
566   - m_off_spec.fill(atom->energy());
567   - }
568   -
569   - } // endfor: looped over all events
570   -
571   - // Store the livetime as exposures of the spectra
572   - m_on_spec.exposure(obs.livetime());
573   - m_off_spec.exposure(obs.livetime());
574   -
575   - // Store the ontime, livetime and deadtime correction in the observation
576   - m_ontime = obs.ontime();
577   - m_livetime = obs.livetime();
578   - m_deadc = obs.deadc();
579   -
580   - // Return
581   - return;
582   -}
583   -
584   -
585   -/***********************************************************************//**
586   - * @brief Compute response for On/Off observation
587   - *
588   - * @param[in] obs CTA observation.
589   - * @param[in] etrue True energy boundaries.
590   - ***************************************************************************/
591   -void GCTAOnOffObservation::compute_response(const GCTAObservation& obs,
592   - const GEbounds& etrue)
593   -{
594   - // Compute response components
595   - compute_arf(obs);
596   - compute_bgd(obs);
597   - compute_rmf(obs, etrue);
598   - compute_alpha(obs);
  402 + // Set Rmf parameter
  403 + par = gammalib::xml_need_par(G_WRITE, xml, "Rmf");
  404 + par->attribute("file", m_rmf.filename());
599 405  
600 406 // Return
601 407 return;
... ... @@ -622,8 +428,8 @@ void GCTAOnOffObservation::compute_response(const GCTAObservation&amp; obs,
622 428 * \f$N_{\gamma}\f$ and background events \f$N_{\rm bgd}\f$ are
623 429 * computed from the spectral models of the relevant components in the
624 430 * model container (spatial and temporal components are ignored so far).
625   - * See the model_on() and model_off() methods for details about the
626   - * model computations.
  431 + * See the N_gamma() and N_bgd() methods for details about the model
  432 + * computations.
627 433 *
628 434 * The log-likelihood is given by
629 435 *
... ... @@ -758,8 +564,8 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
758 564  
759 565 // Get number of gamma and background events (and corresponding
760 566 // spectral model gradients)
761   - double ngam = model_on(models, i, &sky_grad);
762   - double nbgd = model_off(models, i, &bgd_grad);
  567 + double ngam = N_gamma(models, i, &sky_grad);
  568 + double nbgd = N_bgd(models, i, &bgd_grad);
763 569  
764 570 // Skip bin if model is too small (avoids -Inf or NaN gradients)
765 571 double nonpred = ngam + alpha[i] * nbgd;
... ... @@ -1028,119 +834,62 @@ void GCTAOnOffObservation::free_members(void)
1028 834  
1029 835  
1030 836 /***********************************************************************//**
1031   - * @brief Compute vector of alpha parameters
  837 + * @brief Set On/Off observation from a CTA observation
1032 838 *
1033 839 * @param[in] obs CTA observation.
1034 840 *
1035 841 * @exception GException::invalid_value
1036   - * No CTA response found in CTA observation.
1037   - *
1038   - * Compute the alpha parameters for all energy bins. The alpha parameter
1039   - * gives the ratio between the On and Off region background acceptance
1040   - * multiplied by the ratio between On and Off region solid angles.
  842 + * No CTA event list found in CTA observation.
1041 843 *
1042   - * @todo Implement for general sky regions, not only circles
1043   - * @todo Integrate over the sky region instead of taking the values at the
1044   - * sky region centre
  844 + * Sets an On/Off observation from a CTA observation by filling the events
  845 + * that fall in the On and Off regions into the PHA spectra and by computing
  846 + * the corresponding ARF and RMF response functions.
1045 847 ***************************************************************************/
1046   -void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs)
  848 +void GCTAOnOffObservation::set(const GCTAObservation& obs)
1047 849 {
1048   - // Get energy boundaries from On spectrum
1049   - GEbounds ereco = m_on_spec.ebounds();
1050   -
1051   - // Continue only if there are spectral bins
1052   - int nreco = ereco.size();
1053   - if (nreco > 0) {
  850 + // Get CTA event list pointer
  851 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  852 + if (events == NULL) {
  853 + std::string msg = "No event list found in CTA observation \""+
  854 + obs.name()+"\" (ID="+obs.id()+"). ON/OFF observation "
  855 + "can only be filled from event list.";
  856 + throw GException::invalid_value(G_SET, msg);
  857 + }
1054 858  
1055   - // Initialise On/Off exposure ratios
1056   - std::vector<double> alpha(nreco, 0.0);
1057   -
1058   - // Get CTA response pointer. Throw an exception if no response is found
1059   - const GCTAResponseIrf* response =
1060   - dynamic_cast<const GCTAResponseIrf*>(obs.response());
1061   - if (response == NULL) {
1062   - std::string msg = "Response in CTA observation \""+obs.name()+"\" "
1063   - "(ID="+obs.id()+") is not of the GCTAResponseIrf "
1064   - "type.";
1065   - throw GException::invalid_value(G_COMPUTE_ALPHA, msg);
1066   - }
1067   -
1068   - // Get CTA observation pointing direction, zenith, and azimuth
1069   - GCTAPointing obspnt = obs.pointing();
1070   - GSkyDir obsdir = obspnt.dir();
  859 + // Loop over all events
  860 + for (int i = 0; i < events->size(); ++i) {
1071 861  
1072   - // Loop over reconstructed energies
1073   - for (int i = 0; i < nreco; ++i) {
  862 + // Get measured event direction
  863 + const GCTAEventAtom* atom = (*events)[i];
  864 + GSkyDir dir = atom->dir().dir();
  865 +
  866 + // Fill in spectrum according to region containment
  867 + if (m_on_regions.contains(dir)) {
  868 + m_on_spec.fill(atom->energy());
  869 + }
  870 + if (m_off_regions.contains(dir)) {
  871 + m_off_spec.fill(atom->energy());
  872 + }
1074 873  
1075   - // Get mean log10 energy in TeV of bin
1076   - double logEreco = ereco.elogmean(i).log10TeV();
1077   -
1078   - // Initialise effective area totals
1079   - double aon(0.0);
1080   - double aoff(0.0);
1081   -
1082   - // Loop over ON sky regions
1083   - for (int m = 0; m < m_on_regions.size(); ++m) {
1084   -
1085   - // If region is of type GSkyRegionCircle
1086   - const GSkyRegionCircle* skyreg =
1087   - dynamic_cast<const GSkyRegionCircle*>(m_on_regions[m]);
1088   - if (skyreg != NULL) {
1089   -
1090   - // Get centre of circular region
1091   - GSkyDir centreg = skyreg->centre();
1092   -
1093   - // Compute instrument direction from region centre
1094   - GCTAInstDir instdir(obspnt.instdir(centreg));
  874 + } // endfor: looped over all events
1095 875  
1096   - // Add up acceptance
1097   - aon += (*response->background())(logEreco,
1098   - instdir.detx(),
1099   - instdir.dety()) *
1100   - skyreg->solidangle();
1101   -
1102   - } // endif: sky region is of type GSkyRegionCircle
1103   -
1104   - } // Looped over ON regions
1105   -
1106   - // Loop over OFF sky regions
1107   - for (int m = 0; m < m_off_regions.size(); ++m) {
1108   -
1109   - // If region is of type GSkyRegionCircle
1110   - const GSkyRegionCircle* skyreg =
1111   - dynamic_cast<const GSkyRegionCircle*>(m_off_regions[m]);
1112   - if (skyreg != NULL) {
1113   -
1114   - // Get centre of circular region
1115   - GSkyDir centreg = skyreg->centre();
1116   -
1117   - // Compute instrument direction from region centre
1118   - GCTAInstDir instdir(obspnt.instdir(centreg));
  876 + // Store the livetime as exposures of the spectra
  877 + m_on_spec.exposure(obs.livetime());
  878 + m_off_spec.exposure(obs.livetime());
1119 879  
1120   - // Add up acceptance
1121   - aoff += (*response->background())(logEreco,
1122   - instdir.detx(),
1123   - instdir.dety()) *
1124   - skyreg->solidangle();
1125   -
1126   - } // endif: sky region is of type GSkyRegionCircle
1127   -
1128   - } // endfor: looped over OFF regions
1129   -
1130   - // Compute alpha for this energy bin
1131   - if (aoff > 0.0) {
1132   - alpha[i] = aon/aoff;
1133   - }
1134   -
1135   - } // endfor: looped over reconstructed energies
  880 + // Store the ontime, livetime and deadtime correction in the observation
  881 + m_ontime = obs.ontime();
  882 + m_livetime = obs.livetime();
  883 + m_deadc = obs.deadc();
1136 884  
1137   - // Append alpha vector to ARF
1138   - m_arf.append("ALPHA", alpha);
  885 + // Compute response components
  886 + compute_arf(obs);
  887 + compute_bgd(obs);
  888 + compute_alpha(obs);
  889 + compute_rmf(obs);
1139 890  
1140   - } // endif: there were energy bins
1141   -
1142   - // Return
1143   - return;
  891 + // Return
  892 + return;
1144 893 }
1145 894  
1146 895  
... ... @@ -1158,11 +907,11 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs)
1158 907 ***************************************************************************/
1159 908 void GCTAOnOffObservation::compute_arf(const GCTAObservation& obs)
1160 909 {
1161   - // Get energy boundaries from on spectrum
1162   - GEbounds ereco = m_on_spec.ebounds();
1163   -
1164   - // Continue only if there are spectral bins
1165   - int nreco = ereco.size();
  910 + // Get reconstructed energy boundaries from on ARF
  911 + GEbounds ereco = m_arf.ebounds();
  912 + int nreco = ereco.size();
  913 +
  914 + // Continue only if there are ARF bins
1166 915 if (nreco > 0) {
1167 916  
1168 917 // Get CTA response pointer. Throw an exception if no response is found
... ... @@ -1181,9 +930,6 @@ void GCTAOnOffObservation::compute_arf(const GCTAObservation&amp; obs)
1181 930 double zenith = obspnt.zenith();
1182 931 double azimuth = obspnt.azimuth();
1183 932  
1184   - // Initialize ARF
1185   - m_arf = GArf(ereco);
1186   -
1187 933 // Loop over reconstructed energies
1188 934 for (int i = 0; i < nreco; ++i) {
1189 935  
... ... @@ -1202,11 +948,11 @@ void GCTAOnOffObservation::compute_arf(const GCTAObservation&amp; obs)
1202 948 if (skyreg != NULL) {
1203 949  
1204 950 // Get centre of circular region
1205   - GSkyDir centreg(skyreg->centre());
  951 + GSkyDir centre(skyreg->centre());
1206 952  
1207 953 // Compute position of region centre in instrument coordinates
1208   - double theta = obsdir.dist(centreg);
1209   - double phi = obsdir.posang(centreg);
  954 + double theta = obsdir.dist(centre);
  955 + double phi = obsdir.posang(centre);
1210 956  
1211 957 // Add up effective area
1212 958 m_arf[i] += response->aeff(theta,
... ... @@ -1229,33 +975,35 @@ void GCTAOnOffObservation::compute_arf(const GCTAObservation&amp; obs)
1229 975  
1230 976  
1231 977 /***********************************************************************//**
1232   - * @brief Compute background rate in Off observations
  978 + * @brief Compute background rate in Off region
1233 979 *
1234 980 * @param[in] obs CTA observation.
1235 981 *
1236 982 * @exception GException::invalid_argument
1237   - * Specified observation does not contain relevant response or
1238   - * background information
  983 + * Observation does not contain relevant response or background
  984 + * information
  985 + *
  986 + * Compute the background rate in units of events/s/MeV in all Off regions
  987 + * and stores the background rate as additional column with name
  988 + * "BACKGROUND" in the Auxiliary Response File (ARF).
1239 989 *
1240 990 * @todo Should implement integration of background model over region.
1241 991 * @todo Allow general region, not only circles.
1242 992 ***************************************************************************/
1243 993 void GCTAOnOffObservation::compute_bgd(const GCTAObservation& obs)
1244 994 {
1245   - // Get energy boundaries from on spectrum
1246   - GEbounds ereco = m_on_spec.ebounds();
  995 + // Get reconstructed energy boundaries from on ARF
  996 + GEbounds ereco = m_arf.ebounds();
  997 + int nreco = ereco.size();
1247 998  
1248   - // Continue only if there are spectral bins
1249   - int nreco = ereco.size();
  999 + // Continue only if there are ARF bins
1250 1000 if (nreco > 0) {
1251 1001  
1252   - // Initialise background rayes
  1002 + // Initialise background rates to zero
1253 1003 std::vector<double> background(nreco, 0.0);
1254 1004  
1255   - // Get CTA observation pointing direction, zenith, and azimuth
1256   - GCTAPointing obspnt(obs.pointing());
1257   - double zenith = obspnt.zenith();
1258   - double azimuth = obspnt.azimuth();
  1005 + // Get CTA observation pointing direction
  1006 + GCTAPointing obspnt = obs.pointing();
1259 1007  
1260 1008 // Get pointer on CTA IRF response
1261 1009 const GCTAResponseIrf* rsp =
... ... @@ -1266,7 +1014,7 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs)
1266 1014 throw GException::invalid_argument(G_COMPUTE_BGD, msg);
1267 1015 }
1268 1016  
1269   - // Retrieve pointer to CTA background
  1017 + // Get pointer to CTA background
1270 1018 const GCTABackground* bgd = rsp->background();
1271 1019 if (bgd == NULL) {
1272 1020 std::string msg = "Specified observation contains no "
... ... @@ -1280,24 +1028,24 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs)
1280 1028 // Get region centre direction, fall through if region is not
1281 1029 // of type GSkyRegionCircle
1282 1030 const GSkyRegionCircle* skyreg =
1283   - dynamic_cast<const GSkyRegionCircle*>(m_off_regions[m]);
  1031 + dynamic_cast<const GSkyRegionCircle*>(m_off_regions[m]);
1284 1032 if (skyreg == NULL) {
1285 1033 continue;
1286 1034 }
1287 1035  
1288 1036 // Get centre of circular region
1289   - GSkyDir centreg = skyreg->centre();
  1037 + GSkyDir centre = skyreg->centre();
1290 1038  
1291   - // Set direction to centre of the region in instrument coordinates
1292   - GCTAInstDir offdir = obspnt.instdir(centreg);
  1039 + // Get instrument direction of region centre
  1040 + GCTAInstDir offdir = obspnt.instdir(centre);
1293 1041  
1294 1042 // Loop over energy bins
1295 1043 for (int i = 0; i < nreco; ++i) {
1296 1044  
1297   - // Get log energy bin mean
  1045 + // Get log10(E/TeV) of mean reconstructed bin energy
1298 1046 double logEreco = m_on_spec.ebounds().elogmean(i).log10TeV();
1299 1047  
1300   - // Get rate in events/s/MeV
  1048 + // Get background rate in events/s/MeV
1301 1049 background[i] += (*bgd)(logEreco,
1302 1050 offdir.detx(),
1303 1051 offdir.dety()) * skyreg->solidangle();
... ... @@ -1317,27 +1065,146 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs)
1317 1065  
1318 1066  
1319 1067 /***********************************************************************//**
  1068 + * @brief Compute vector of alpha parameters
  1069 + *
  1070 + * @param[in] obs CTA observation.
  1071 + *
  1072 + * @exception GException::invalid_value
  1073 + * No CTA response found in CTA observation.
  1074 + *
  1075 + * Compute the alpha parameters for all energy bins. The alpha parameter
  1076 + * gives the ratio between the On and Off region background acceptance
  1077 + * multiplied by the ratio between On and Off region solid angles.
  1078 + *
  1079 + * @todo Implement for general sky regions, not only circles
  1080 + * @todo Integrate over the sky region instead of taking the values at the
  1081 + * sky region centre
  1082 + ***************************************************************************/
  1083 +void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs)
  1084 +{
  1085 + // Get reconstructed energy boundaries from on ARF
  1086 + GEbounds ereco = m_arf.ebounds();
  1087 + int nreco = ereco.size();
  1088 +
  1089 + // Continue only if there are ARF bins
  1090 + if (nreco > 0) {
  1091 +
  1092 + // Initialise On/Off exposure ratios
  1093 + std::vector<double> alpha(nreco, 0.0);
  1094 +
  1095 + // Get CTA response pointer. Throw an exception if no response is found
  1096 + const GCTAResponseIrf* response =
  1097 + dynamic_cast<const GCTAResponseIrf*>(obs.response());
  1098 + if (response == NULL) {
  1099 + std::string msg = "Response in CTA observation \""+obs.name()+"\" "
  1100 + "(ID="+obs.id()+") is not of the GCTAResponseIrf "
  1101 + "type.";
  1102 + throw GException::invalid_value(G_COMPUTE_ALPHA, msg);
  1103 + }
  1104 +
  1105 + // Get CTA observation pointing direction, zenith, and azimuth
  1106 + GCTAPointing obspnt = obs.pointing();
  1107 + GSkyDir obsdir = obspnt.dir();
  1108 +
  1109 + // Loop over reconstructed energies
  1110 + for (int i = 0; i < nreco; ++i) {
  1111 +
  1112 + // Get mean log10 energy in TeV of bin
  1113 + double logEreco = ereco.elogmean(i).log10TeV();
  1114 +
  1115 + // Initialise effective area totals
  1116 + double aon(0.0);
  1117 + double aoff(0.0);
  1118 +
  1119 + // Loop over ON sky regions
  1120 + for (int m = 0; m < m_on_regions.size(); ++m) {
  1121 +
  1122 + // If region is of type GSkyRegionCircle
  1123 + const GSkyRegionCircle* skyreg =
  1124 + dynamic_cast<const GSkyRegionCircle*>(m_on_regions[m]);
  1125 + if (skyreg != NULL) {
  1126 +
  1127 + // Get centre of circular region
  1128 + GSkyDir centreg = skyreg->centre();
  1129 +
  1130 + // Compute instrument direction from region centre
  1131 + GCTAInstDir instdir(obspnt.instdir(centreg));
  1132 +
  1133 + // Add up acceptance
  1134 + aon += (*response->background())(logEreco,
  1135 + instdir.detx(),
  1136 + instdir.dety()) *
  1137 + skyreg->solidangle();
  1138 +
  1139 + } // endif: sky region is of type GSkyRegionCircle
  1140 +
  1141 + } // Looped over ON regions
  1142 +
  1143 + // Loop over OFF sky regions
  1144 + for (int m = 0; m < m_off_regions.size(); ++m) {
  1145 +
  1146 + // If region is of type GSkyRegionCircle
  1147 + const GSkyRegionCircle* skyreg =
  1148 + dynamic_cast<const GSkyRegionCircle*>(m_off_regions[m]);
  1149 + if (skyreg != NULL) {
  1150 +
  1151 + // Get centre of circular region
  1152 + GSkyDir centreg = skyreg->centre();
  1153 +
  1154 + // Compute instrument direction from region centre
  1155 + GCTAInstDir instdir(obspnt.instdir(centreg));
  1156 +
  1157 + // Add up acceptance
  1158 + aoff += (*response->background())(logEreco,
  1159 + instdir.detx(),
  1160 + instdir.dety()) *
  1161 + skyreg->solidangle();
  1162 +
  1163 + } // endif: sky region is of type GSkyRegionCircle
  1164 +
  1165 + } // endfor: looped over OFF regions
  1166 +
  1167 + // Compute alpha for this energy bin
  1168 + if (aoff > 0.0) {
  1169 + alpha[i] = aon/aoff;
  1170 + }
  1171 +
  1172 + } // endfor: looped over reconstructed energies
  1173 +
  1174 + // Append alpha vector to ARF
  1175 + m_arf.append("ALPHA", alpha);
  1176 +
  1177 + } // endif: there were energy bins
  1178 +
  1179 + // Return
  1180 + return;
  1181 +}
  1182 +
  1183 +
  1184 +/***********************************************************************//**
1320 1185 * @brief Compute RMF of On/Off observation
1321 1186 *
1322 1187 * @param[in] obs CTA observation.
1323   - * @param[in] etrue True energy boundaries.
1324 1188 *
1325 1189 * @exception GException::invalid_value
1326 1190 * Observation does not contain IRF response
1327 1191 *
  1192 + * Compute the energy redistribution matrix for an On/Off observation. The
  1193 + * method requires that the RMF energy axes have been defined before.
  1194 + *
1328 1195 * @todo Check if we really should use reconstructed energy in the Aeff
1329 1196 * access.
1330 1197 * @todo Generalise to any kind of sky region
1331 1198 ***************************************************************************/
1332   -void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs,
1333   - const GEbounds& etrue)
  1199 +void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs)
1334 1200 {
1335   - // Get reconstructed energy boundaries from on spectrum
1336   - GEbounds ereco = m_on_spec.ebounds();
1337   -
1338   - // Continue only if there are spectral bins
1339   - int ntrue = etrue.size();
1340   - int nreco = ereco.size();
  1201 + // Get true and reconstructed energy boundaries from RMF
  1202 + GEbounds etrue = m_rmf.etrue();
  1203 + GEbounds ereco = m_rmf.emeasured();
  1204 + int ntrue = etrue.size();
  1205 + int nreco = ereco.size();
  1206 +
  1207 + // Continue only if there are RMF bins
1341 1208 if (ntrue > 0 && nreco > 0) {
1342 1209  
1343 1210 // Get CTA response pointer
... ... @@ -1356,9 +1223,6 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
1356 1223 double zenith = obspnt.zenith();
1357 1224 double azimuth = obspnt.azimuth();
1358 1225  
1359   - // Initialize RMF
1360   - m_rmf = GRmf(etrue, ereco);
1361   -
1362 1226 // Loop over reconstructed energy
1363 1227 for (int ireco = 0; ireco < nreco; ++ireco) {
1364 1228  
... ... @@ -1367,6 +1231,9 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
1367 1231  
1368 1232 // Loop over true energy
1369 1233 for (int itrue = 0; itrue < ntrue; ++itrue) {
  1234 +
  1235 + // Initialise RMF element
  1236 + m_rmf(itrue, ireco) = 0.0;
1370 1237  
1371 1238 // Compute true energy and initialise weight for this bin
1372 1239 double logEtrue = etrue.elogmean(itrue).log10TeV();
... ... @@ -1424,33 +1291,46 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
1424 1291  
1425 1292  
1426 1293 /***********************************************************************
1427   - * @brief Evaluate model for sky contribution and fill model gradients
  1294 + * @brief Compute \f$N_{\gamma}\f$ value and model parameter gradients
1428 1295 *
1429 1296 * @param[in] models Model container.
1430   - * @param[in,out] ibin Energy bin number.
1431   - * @param[in,out] mod_grad Model gradient array.
  1297 + * @param[in] ibin Energy bin number.
  1298 + * @param[in,out] grad Model gradient vector.
1432 1299 *
1433   - * Computes the number of expected gamma-ray events in an On region for a
1434   - * given energy bin, and returns it. Also computes the gradients of the
1435   - * sky model spectral components. The method assumes that
1436   - * parameters are stored in the order spatial-spectral-temporal.
  1300 + * Returns the predicted number of source events \f$N_{\gamma}\f$
  1301 + * in the On regions for a given energy bin. The method computes also
  1302 + *
  1303 + * \f[
  1304 + * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
  1305 + * \f]
  1306 + *
  1307 + * which are the gradients in the predicted number of source events with
  1308 + * respect to all model parameters.
  1309 + *
  1310 + * The method assumes that parameters are stored in the order
  1311 + * spatial-spectral-temporal.
1437 1312 *
1438 1313 * @todo I think this method only works for point sources. What happens
1439 1314 * for an extended source?
1440 1315 ***********************************************************************/
1441   -double GCTAOnOffObservation::model_on(const GModels& models,
1442   - int ibin,
1443   - GVector* mod_grad) const
  1316 +double GCTAOnOffObservation::N_gamma(const GModels& models,
  1317 + const int& ibin,
  1318 + GVector* grad) const
1444 1319 {
1445 1320 // Get total number of model parameters
1446 1321 int npars = models.npars();
1447   -
1448   - // Initialize variables (vector has 0.0 values)
  1322 +
  1323 + // Initialize results
1449 1324 double value = 0.0;
1450   - int ipar = 0;
  1325 + for (int i = 0; i < npars; ++i) {
  1326 + (*grad)[i] = 0.0;
  1327 + }
1451 1328  
1452 1329 // Continue only if bin number is in range and there are model parameters
1453   - if ((ibin < m_on_spec.size()) && (npars > 0)) {
  1330 + if ((ibin >= 0) && (ibin < m_on_spec.size()) && (npars > 0)) {
  1331 +
  1332 + // Initialise parameter index
  1333 + int ipar = 0;
1454 1334  
1455 1335 // Get energy bin bounds
1456 1336 const GEnergy emin = m_on_spec.ebounds().emin(ibin);
... ... @@ -1460,8 +1340,8 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
1460 1340  
1461 1341 // Compute normalisation factors
1462 1342 double exposure = m_on_spec.exposure();
1463   - double norm_flux = m_arf[ibin] * exposure;
1464   - double norm_grad = norm_flux * ewidth;
  1343 + double norm_flux = m_arf[ibin] * exposure; // cm2 s
  1344 + double norm_grad = norm_flux * ewidth; // cm2 s MeV
1465 1345  
1466 1346 // Loop over models
1467 1347 for (int j = 0; j < models.size(); ++j) {
... ... @@ -1500,7 +1380,7 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
1500 1380 // computing the flux over the energy bin in ph/cm2/s
1501 1381 // and multiplying this flux by the effective area (cm2)
1502 1382 // and the livetime (s)
1503   - value += spectral->flux(emin,emax) * norm_flux;
  1383 + value += spectral->flux(emin, emax) * norm_flux;
1504 1384  
1505 1385 // Determine the model gradients at the current energy
1506 1386 spectral->eval_gradients(emean);
... ... @@ -1509,7 +1389,7 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
1509 1389 for (int k = 0; k < spectral->size(); ++k, ++ipar) {
1510 1390 GModelPar& par = (*spectral)[k];
1511 1391 if (par.is_free() && ipar < npars) {
1512   - (*mod_grad)[ipar] = par.factor_gradient() * norm_grad;
  1392 + (*grad)[ipar] = par.factor_gradient() * norm_grad;
1513 1393 }
1514 1394 }
1515 1395  
... ... @@ -1531,47 +1411,58 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
1531 1411  
1532 1412  
1533 1413 /***********************************************************************
1534   - * @brief Evaluate model for bgd contribution and fill model gradients
  1414 + * @brief Compute \f$N_{\rm bgd}\f$ value and model parameter gradients
1535 1415 *
1536 1416 * @param[in] models Model container.
1537   - * @param[in,out] ibin Energy bin number.
1538   - * @param[in,out] mod_grad Model gradient array.
1539   - * @return Number of expected background events.
  1417 + * @param[in] ibin Energy bin index.
  1418 + * @param[in,out] grad Model gradient vector.
  1419 + * @return Predicted number of background events in Off regions.
1540 1420 *
1541   - * Returns the number of expected background events in the Off spectrum
1542   - * for a given energy bin. The method computes also the gradients of the
1543   - * background model spectral components.
  1421 + * Returns the predicted number of background events \f$N_{\rm bgd}\f$
  1422 + * in the Off regions for a given energy bin. The method computes also
  1423 + *
  1424 + * \f[
  1425 + * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
  1426 + * \f]
  1427 + *
  1428 + * which are the gradients in the predicted number of background events
  1429 + * with respect to all model parameters.
1544 1430 *
1545 1431 * The method assumes that the model parameters are stored in the order
1546   - * spatial-spectral-temporal.
  1432 + * spectral-temporal.
1547 1433 *
1548 1434 * @todo So far only handles the GCTAModelIrfBackground background model;
1549 1435 * Should be more generic, but there is maybe for the moment no
1550 1436 * other choice than implementing something for all possible
1551 1437 * background models; that's not very satisfactory ...
1552 1438 ***********************************************************************/
1553   -double GCTAOnOffObservation::model_off(const GModels& models,
1554   - int ibin,
1555   - GVector* mod_grad) const
  1439 +double GCTAOnOffObservation::N_bgd(const GModels& models,
  1440 + const int& ibin,
  1441 + GVector* grad) const
1556 1442 {
1557 1443 // Get total number of model parameters
1558 1444 int npars = models.npars();
1559 1445  
1560   - // Initialize variables
1561   - double value = 0.0;
1562   - int ipar = 0;
1563   -
1564   - // Continue only if bin number is in range and if there are parameters
1565   - if ((ibin < m_off_spec.size()) && (npars > 0)) {
  1446 + // Initialize results
  1447 + double value = 0.0;
  1448 + for (int i = 0; i < npars; ++i) {
  1449 + (*grad)[i] = 0.0;
  1450 + }
  1451 +
  1452 + // Continue only if bin number is valid and if there are model parameters
  1453 + if ((ibin >= 0) && (ibin < m_off_spec.size()) && (npars > 0)) {
1566 1454  
1567   - // Get reference to background rates
  1455 + // Initialise parameter index
  1456 + int ipar = 0;
  1457 +
  1458 + // Get reference to background rates (events/MeV/s)
1568 1459 const std::vector<double>& background = m_arf["BACKGROUND"];
1569 1460  
1570 1461 // Get energy bin mean and width
1571 1462 const GEnergy emean = m_off_spec.ebounds().elogmean(ibin);
1572 1463 const double ewidth = m_off_spec.ebounds().ewidth(ibin).MeV();
1573 1464  
1574   - // Compute normalisation factor
  1465 + // Compute normalisation factor (events)
1575 1466 double exposure = m_off_spec.exposure();
1576 1467 double norm = background[ibin] * exposure * ewidth;
1577 1468  
... ... @@ -1591,7 +1482,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
1591 1482 continue;
1592 1483 }
1593 1484  
1594   - // Fall through if model is not a IRF background component
  1485 + // Fall through if model is not an IRF background component
1595 1486 const GCTAModelIrfBackground* bgd =
1596 1487 dynamic_cast<const GCTAModelIrfBackground*>(mptr);
1597 1488 if (bgd == NULL) {
... ... @@ -1604,16 +1495,17 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
1604 1495 if (spectral != NULL) {
1605 1496  
1606 1497 // Determine the number of background events in model by
1607   - // computing the rate at the energy bin in events/s/MeV and
1608   - // by multiplying the rate with the off livetime (s), the
1609   - // energy bin width (MeV) and the spectral normalisation
  1498 + // computing the model normalization at the mean bin energy bin
  1499 + // and multiplying the normalisation with the number of
  1500 + // background events
1610 1501 value += spectral->eval_gradients(emean) * norm;
1611 1502  
1612   - // Loop over spectral model parameters
  1503 + // Compute the parameter gradients for all spectral model
  1504 + // parameters
1613 1505 for (int k = 0; k < spectral->size(); ++k, ++ipar) {
1614   - GModelPar& par=(*spectral)[k];
  1506 + GModelPar& par = (*spectral)[k];
1615 1507 if (par.is_free() && ipar < npars) {
1616   - (*mod_grad)[ipar] = par.factor_gradient() * norm;
  1508 + (*grad)[ipar] += par.factor_gradient() * norm;
1617 1509 }
1618 1510 }
1619 1511  
... ...
inst/cta/test/test_CTA.py
  View file @f54cd30
... ... @@ -265,15 +265,12 @@ class Test(gammalib.GPythonTestSuite):
265 265 ereco = gammalib.GEbounds(10, gammalib.GEnergy(0.1, 'TeV'),
266 266 gammalib.GEnergy(10.0, 'TeV'))
267 267  
268   - # Create On/Off observations by filling all events found in
269   - # the observation container and computing the response
  268 + # Create On/Off observations from CTA observations
270 269 filename = self.caldb + '/../data/irf_unbinned.xml'
271 270 inobs = gammalib.GObservations(filename)
272 271 outobs = gammalib.GObservations()
273 272 for run in inobs:
274   - onoff = gammalib.GCTAOnOffObservation(ereco, on, off)
275   - onoff.fill(run)
276   - onoff.compute_response(run, etrue)
  273 + onoff = gammalib.GCTAOnOffObservation(run, etrue, ereco, on, off)
277 274 outobs.append(onoff)
278 275  
279 276 # Load model container and attach it to the observations
... ...