Rolf Buehler / gammalib

Commit 0658d8678a894878935ca427a240916b74f4b130

Authored by Jürgen Knödlseder
1 parent bb345382

Implement logging capability into CTA cube response classes. 

This enhances the user information in ctools.
Inline Side-by-side
Showing 11 changed files with 177 additions and 155 deletions
inst/cta/include/GCTACubeBackground.hpp
  View file @0658d86
... ... @@ -65,7 +65,7 @@ public:
65 65 void clear(void);
66 66 GCTACubeBackground* clone(void) const;
67 67 std::string classname(void) const;
68   - void fill(const GObservations& obs, GLog* log = 0);
  68 + void fill(const GObservations& obs, GLog* log = NULL);
69 69 double integral(const double& logE) const;
70 70 void read(const GFits& fits);
71 71 void write(GFits& file) const;
... ...
inst/cta/include/GCTACubeExposure.hpp
  View file @0658d86
... ... @@ -74,7 +74,7 @@ public:
74 74 GCTACubeExposure* clone(void) const;
75 75 std::string classname(void) const;
76 76 void set(const GCTAObservation& obs);
77   - void fill(const GObservations& obs, GLog* log = 0);
  77 + void fill(const GObservations& obs, GLog* log = NULL);
78 78 const GSkymap& cube(void) const;
79 79 const GEbounds& ebounds(void) const;
80 80 const GGti& gti(void) const;
... ...
inst/cta/include/GCTACubePsf.hpp
  View file @0658d86
... ... @@ -82,7 +82,7 @@ public:
82 82 GCTACubePsf* clone(void) const;
83 83 std::string classname(void) const;
84 84 void set(const GCTAObservation& obs);
85   - void fill(const GObservations& obs, GLog* log = 0);
  85 + void fill(const GObservations& obs, GLog* log = NULL);
86 86 const GSkymap& map(void) const;
87 87 const GEbounds& ebounds(void) const;
88 88 const GNodeArray& deltas(void) const;
... ...
inst/cta/include/GCTAObservation.hpp
  View file @0658d86
... ... @@ -124,12 +124,12 @@ public:
124 124 void deadc(const double& deadc);
125 125 void eventfile(const std::string& filename);
126 126 const std::string& eventfile(void) const;
  127 + const std::string& eventtype(void) const;
127 128 void dispose_events(void);
128 129 const double& lo_user_thres(void) const;
129 130 const double& hi_user_thres(void) const;
130 131 void n_tels(const int& tels);
131 132 const int& n_tels(void) const;
132   - const std::string& eventtype(void) const;
133 133  
134 134 protected:
135 135 // Protected methods
... ... @@ -414,6 +414,18 @@ const std::string& GCTAObservation::eventfile(void) const
414 414  
415 415  
416 416 /***********************************************************************//**
  417 + * @brief Return type of events
  418 + *
  419 + * @return Type of events.
  420 + ***************************************************************************/
  421 +inline
  422 +const std::string& GCTAObservation::eventtype(void) const
  423 +{
  424 + return m_eventtype;
  425 +}
  426 +
  427 +
  428 +/***********************************************************************//**
417 429 * @brief Return user low energy threshold
418 430 *
419 431 * @return User low energy threshold.
... ... @@ -460,15 +472,4 @@ const int& GCTAObservation::n_tels(void) const
460 472 return m_n_tels;
461 473 }
462 474  
463   -
464   -/***********************************************************************//**
465   - * @brief Return type of events
466   - *
467   - * @return Type of events.
468   - ***************************************************************************/
469   -inline
470   -const std::string& GCTAObservation::eventtype(void) const
471   -{
472   - return m_eventtype;
473   -}
474 475 #endif /* GCTAOBSERVATION_HPP */
... ...
inst/cta/pyext/GCTACubeBackground.i
  View file @0658d86
... ... @@ -52,7 +52,7 @@ public:
52 52 void clear(void);
53 53 GCTACubeBackground* clone(void) const;
54 54 std::string classname(void) const;
55   - void fill(const GObservations& obs, GLog* log=0);
  55 + void fill(const GObservations& obs, GLog* log = NULL);
56 56 double integral(const double& logE) const;
57 57 void read(const GFits& fits);
58 58 void write(GFits& file) const;
... ...
inst/cta/pyext/GCTACubeExposure.i
  View file @0658d86
... ... @@ -65,7 +65,7 @@ public:
65 65 GCTACubeExposure* clone(void) const;
66 66 std::string classname(void) const;
67 67 void set(const GCTAObservation& obs);
68   - void fill(const GObservations& obs, GLog* log=0);
  68 + void fill(const GObservations& obs, GLog* log = NULL);
69 69 const GSkymap& cube(void) const;
70 70 const GEbounds& ebounds(void) const;
71 71 const GGti& gti(void) const;
... ...
inst/cta/pyext/GCTACubePsf.i
  View file @0658d86
... ... @@ -74,7 +74,7 @@ public:
74 74 GCTACubePsf* clone(void) const;
75 75 std::string classname(void) const;
76 76 void set(const GCTAObservation& obs);
77   - void fill(const GObservations& obs, GLog* log=0);
  77 + void fill(const GObservations& obs, GLog* log = NULL);
78 78 const GSkymap& map(void) const;
79 79 const GEbounds& ebounds(void) const;
80 80 const GNodeArray& deltas(void) const;
... ...
inst/cta/pyext/GCTAObservation.i
  View file @0658d86
... ... @@ -96,12 +96,12 @@ public:
96 96 void deadc(const double& deadc);
97 97 void eventfile(const std::string& filename);
98 98 const std::string& eventfile(void) const;
  99 + const std::string& eventtype(void) const;
99 100 void dispose_events(void);
100 101 const double& lo_user_thres(void) const;
101 102 const double& hi_user_thres(void) const;
102 103 void n_tels(const int& tels);
103 104 const int& n_tels(void) const;
104   - const std::string& eventtype(void) const;
105 105 };
106 106  
107 107  
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @0658d86
... ... @@ -264,6 +264,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
264 264 * @brief Fill background cube from observation container
265 265 *
266 266 * @param[in] obs Observation container.
  267 + * @param[in] log Pointer to logger (optional).
267 268 *
268 269 * @exception GException::invalid_value
269 270 * No event list found in CTA observations.
... ... @@ -289,60 +290,63 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log)
289 290 // Get observation and continue only if it is a CTA observation
290 291 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
291 292  
292   - if (cta != NULL) {
293   -
294   - // Skip observation if we don't have an unbinned observation
295   - if (cta->eventtype() != "EventList") {
296   -
297   - // Log that we skip the this observation
298   - if (log != NULL) {
299   - *log << "Warning: Skipping binned observation ";
300   - *log << "\"" << cta->name() << "\"" <<std::endl;
301   - }
302   - continue;
  293 + // Skip observation if it's not CTA
  294 + if (cta == NULL) {
  295 + if (log != NULL) {
  296 + *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
  297 + *log << "\"" << obs[i]->name() << "\"" <<std::endl;
303 298 }
  299 + continue;
  300 + }
  301 +
  302 + // Skip observation if we don't have an unbinned observation
  303 + if (cta->eventtype() != "EventList") {
  304 + if (log != NULL) {
  305 + *log << "Warning: Skipping binned observation ";
  306 + *log << "\"" << cta->name() << "\"" <<std::endl;
  307 + }
  308 + continue;
  309 + }
304 310  
305   - // Extract region of interest from CTA observation
306   - GCTARoi roi = cta->roi();
307   -
308   - // Set GTI of actual observations as the GTI of the event cube
309   - eventcube.gti(cta->gti());
  311 + // Extract region of interest from CTA observation
  312 + GCTARoi roi = cta->roi();
310 313  
311   - // Get observation livetime
312   - double livetime = cta->livetime();
  314 + // Set GTI of actual observations as the GTI of the event cube
  315 + eventcube.gti(cta->gti());
313 316  
314   - // Loop over all bins in background cube
315   - for (int i = 0; i < eventcube.size(); ++i) {
  317 + // Get observation livetime
  318 + double livetime = cta->livetime();
316 319  
317   - // Get event bin
318   - GCTAEventBin* bin = eventcube[i];
  320 + // Loop over all bins in background cube
  321 + for (int i = 0; i < eventcube.size(); ++i) {
319 322  
320   - // Continue only if binned in contained in ROI
321   - if (roi.contains(*bin)) {
  323 + // Get event bin
  324 + GCTAEventBin* bin = eventcube[i];
322 325  
323   - // Compute model value for event bin. The model value is
324   - // given in counts/MeV/s/sr.
325   - double model = obs.models().eval(*bin, *cta);
  326 + // Continue only if binned in contained in ROI
  327 + if (roi.contains(*bin)) {
326 328  
327   - // Multiply by livetime to get the correct weighting for
328   - // each observation. We divide by the total livetime later
329   - // to get the background model in units of counts/MeV/s/sr.
330   - model *= livetime;
  329 + // Compute model value for event bin. The model value is
  330 + // given in counts/MeV/s/sr.
  331 + double model = obs.models().eval(*bin, *cta);
331 332  
332   - // Add existing number of counts
333   - model += bin->counts();
  333 + // Multiply by livetime to get the correct weighting for
  334 + // each observation. We divide by the total livetime later
  335 + // to get the background model in units of counts/MeV/s/sr.
  336 + model *= livetime;
334 337  
335   - // Store cumulated value (units: counts/MeV/sr)
336   - bin->counts(model);
  338 + // Add existing number of counts
  339 + model += bin->counts();
337 340  
338   - } // endif: bin was contained in RoI
  341 + // Store cumulated value (units: counts/MeV/sr)
  342 + bin->counts(model);
339 343  
340   - } // endfor: looped over all bins
  344 + } // endif: bin was contained in RoI
341 345  
342   - // Accumulate livetime
343   - total_livetime += livetime;
  346 + } // endfor: looped over all bins
344 347  
345   - } // endif: cta observation was vaild
  348 + // Accumulate livetime
  349 + total_livetime += livetime;
346 350  
347 351 } // endfor: looped over all observations
348 352  
... ...
inst/cta/src/GCTACubeExposure.cpp
  View file @0658d86
... ... @@ -364,6 +364,7 @@ void GCTACubeExposure::set(const GCTAObservation&amp; obs)
364 364 * @brief Fill exposure cube from observation container
365 365 *
366 366 * @param[in] obs Observation container.
  367 + * @param[in] log Pointer to logger (optional).
367 368 *
368 369 * @exception GException::invalid_value
369 370 * No event list found in CTA observations.
... ... @@ -384,66 +385,74 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
384 385  
385 386 // Get observation and continue only if it is a CTA observation
386 387 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
387   - if (cta != NULL) {
388 388  
389   - // Skip observation if we don't have an unbinned observation
390   - if (cta->eventtype() != "EventList") {
391   -
392   - // Log that we skip the this observation
393   - if (log != NULL) {
394   - *log << "Warning: Skipping binned observation ";
395   - *log << "\"" << cta->name() << "\"" <<std::endl;
396   - }
397   - continue;
  389 + // Skip observation if it's not CTA
  390 + if (cta == NULL) {
  391 + if (log != NULL) {
  392 + *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
  393 + *log << "\"" << obs[i]->name() << "\"" <<std::endl;
398 394 }
  395 + continue;
  396 + }
399 397  
400   - // Extract region of interest from CTA observation
401   - GCTARoi roi = cta->roi();
  398 + // Skip observation if we don't have an unbinned observation
  399 + if (cta->eventtype() != "EventList") {
  400 + if (log != NULL) {
  401 + *log << "Warning: Skipping binned observation ";
  402 + *log << "\"" << cta->name() << "\"" <<std::endl;
  403 + }
  404 + continue;
  405 + }
402 406  
403   - // Get references on CTA response and pointing direction
404   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
405   - const GSkyDir& pnt = cta->pointing().dir();
  407 + // Extract region of interest from CTA observation
  408 + GCTARoi roi = cta->roi();
406 409  
407   - // Continue only if response is valid
408   - if (rsp != NULL) {
  410 + // Get references on CTA response and pointing direction
  411 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  412 + const GSkyDir& pnt = cta->pointing().dir();
  413 +
  414 + // Skip observation if we don't have an unbinned observation
  415 + if (rsp == NULL) {
  416 + if (log != NULL) {
  417 + *log << "Warning: Observation \"" << cta->name();
  418 + *log << "\" contains no IRF response." <<std::endl;
  419 + }
  420 + continue;
  421 + }
409 422  
410   - // Loop over all pixels in sky map
411   - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
  423 + // Loop over all pixels in sky map
  424 + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
412 425  
413   - // Get pixel sky direction
414   - GSkyDir dir = m_cube.inx2dir(pixel);
  426 + // Get pixel sky direction
  427 + GSkyDir dir = m_cube.inx2dir(pixel);
415 428  
416   - // Continue only if pixel is within RoI
417   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
  429 + // Continue only if pixel is within RoI
  430 + if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
418 431  
419   - // Compute theta angle with respect to pointing
420   - // direction in radians
421   - double theta = pnt.dist(dir);
422   -
423   - // Loop over all exposure cube energy bins
424   - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
  432 + // Compute theta angle with respect to pointing
  433 + // direction in radians
  434 + double theta = pnt.dist(dir);
425 435  
426   - // Get logE/TeV
427   - double logE = m_ebounds.elogmean(iebin).log10TeV();
  436 + // Loop over all exposure cube energy bins
  437 + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
428 438  
429   - // Add to exposure cube (effective area * livetime)
430   - m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
431   - cta->livetime();
  439 + // Get logE/TeV
  440 + double logE = m_ebounds.elogmean(iebin).log10TeV();
432 441  
433   - } // endfor: looped over energy bins
  442 + // Add to exposure cube (effective area * livetime)
  443 + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  444 + cta->livetime();
434 445  
435   - } // endif: pixel within RoI
  446 + } // endfor: looped over energy bins
436 447  
437   - } // endfor: looped over all pixels
  448 + } // endif: pixel within RoI
438 449  
439   - // Append GTIs and increment livetime
440   - m_gti.extend(cta->gti());
441   - m_livetime += cta->livetime();
442   -
443   - } // endif: response was valid
444   -
445   - } // endif: observation was a CTA observation
  450 + } // endfor: looped over all pixels
446 451  
  452 + // Append GTIs and increment livetime
  453 + m_gti.extend(cta->gti());
  454 + m_livetime += cta->livetime();
  455 +
447 456 } // endfor: looped over observations
448 457  
449 458 // Return
... ...
inst/cta/src/GCTACubePsf.cpp
  View file @0658d86
... ... @@ -450,79 +450,87 @@ void GCTACubePsf::fill(const GObservations&amp; obs, GLog* log)
450 450  
451 451 // Get observation and continue only if it is a CTA observation
452 452 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
453   - if (cta != NULL) {
454 453  
455   - // Skip observation if we don't have an unbinned observation
456   - if (cta->eventtype() != "EventList") {
457   -
458   - // Log that we skip the this observation
459   - if (log != NULL) {
460   - *log << "Warning: Skipping binned observation ";
461   - *log << "\"" << cta->name() << "\"" << std::endl;
462   - }
463   - continue;
  454 + // Skip observation if it's not CTA
  455 + if (cta == NULL) {
  456 + if (log != NULL) {
  457 + *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
  458 + *log << "\"" << obs[i]->name() << "\"" <<std::endl;
464 459 }
  460 + continue;
  461 + }
465 462  
466   - // Extract region of interest from CTA observation
467   - GCTARoi roi = cta->roi();
  463 + // Skip observation if we don't have an unbinned observation
  464 + if (cta->eventtype() != "EventList") {
  465 + if (log != NULL) {
  466 + *log << "Warning: Skipping binned observation ";
  467 + *log << "\"" << cta->name() << "\"" <<std::endl;
  468 + }
  469 + continue;
  470 + }
468 471  
469   - // Get references on CTA response and pointing direction
470   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
471   - const GSkyDir& pnt = cta->pointing().dir();
  472 + // Extract region of interest from CTA observation
  473 + GCTARoi roi = cta->roi();
472 474  
473   - // Continue only if response is valid
474   - if (rsp != NULL) {
  475 + // Get references on CTA response and pointing direction
  476 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  477 + const GSkyDir& pnt = cta->pointing().dir();
  478 +
  479 + // Skip observation if we don't have an unbinned observation
  480 + if (rsp == NULL) {
  481 + if (log != NULL) {
  482 + *log << "Warning: Observation \"" << cta->name();
  483 + *log << "\" contains no IRF response." <<std::endl;
  484 + }
  485 + continue;
  486 + }
475 487  
476   - // Loop over all pixels in sky map
477   - for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
  488 + // Loop over all pixels in sky map
  489 + for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
478 490  
479   - // Get pixel sky direction
480   - GSkyDir dir = m_cube.inx2dir(pixel);
  491 + // Get pixel sky direction
  492 + GSkyDir dir = m_cube.inx2dir(pixel);
481 493  
482   - // Continue only if pixel is within RoI
483   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
  494 + // Continue only if pixel is within RoI
  495 + if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
484 496  
485   - // Compute theta angle with respect to pointing
486   - // direction in radians
487   - double theta = pnt.dist(dir);
488   -
489   - // Loop over all energy bins
490   - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
491   -
492   - // Get logE/TeV
493   - double logE = m_ebounds.elogmean(iebin).log10TeV();
  497 + // Compute theta angle with respect to pointing
  498 + // direction in radians
  499 + double theta = pnt.dist(dir);
494 500  
495   - // Compute exposure weight
496   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
497   - cta->livetime();
  501 + // Loop over all energy bins
  502 + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
498 503  
499   - // Accumulate weights
500   - exposure(pixel, iebin) += weight;
  504 + // Get logE/TeV
  505 + double logE = m_ebounds.elogmean(iebin).log10TeV();
501 506  
502   - // Loop over delta values
503   - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
  507 + // Compute exposure weight
  508 + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  509 + cta->livetime();
504 510  
505   - // Compute delta in radians
506   - double delta = m_deltas[idelta] * gammalib::deg2rad;
  511 + // Accumulate weights
  512 + exposure(pixel, iebin) += weight;
507 513  
508   - // Set map index
509   - int imap = offset(idelta, iebin);
  514 + // Loop over delta values
  515 + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
510 516  
511   - // Add on PSF cube
512   - m_cube(pixel, imap) +=
513   - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
  517 + // Compute delta in radians
  518 + double delta = m_deltas[idelta] * gammalib::deg2rad;
514 519  
515   - } // endfor: looped over delta bins
  520 + // Set map index
  521 + int imap = offset(idelta, iebin);
516 522  
517   - } // endfor: looped over energy bins
  523 + // Add on PSF cube
  524 + m_cube(pixel, imap) +=
  525 + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
518 526  
519   - } // endif: pixel was within RoI
  527 + } // endfor: looped over delta bins
520 528  
521   - } // endfor: looped over all pixels
  529 + } // endfor: looped over energy bins
522 530  
523   - } // endif: response was valid
  531 + } // endif: pixel was within RoI
524 532  
525   - } // endif: observation was a CTA observation
  533 + } // endfor: looped over all pixels
526 534  
527 535 } // endfor: looped over observations
528 536  
... ...