Samuel Flis / gammalib

Commit b447670c8374f46535284f555a97b94a0aacf6ba

Authored by Jürgen Knödlseder
1 parent 555c1e7c

No data space selection in CTA cube response computation 

The CTA cube response classes do no longer consider Region of Interest and energy boundary selections as they are now handled by the cube weighting.
Inline Side-by-side
Showing 7 changed files with 138 additions and 305 deletions
ChangeLog
  View file @b447670
1   -2016-06-20 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-06-21 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + CTA cube response no longer selects on RoI and energies (#1789)
6 7 Add weights to CTA event cube and event bin (#1789)
7 8 Remove GCTACubeExposure, GCTACubePsf and GCTACubeBackground elogmeans() methods
8 9 GCTACubeExposure, GCTACubePsf and GCTACubeBackground use energy values (#1788)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -20 June 2016
  5 +21 June 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -408,6 +408,10 @@ multiplied to the bin size in the GCTAEventBin::size() method so that model valu
408 408 weigthed according to the bin weight. This enables the usage of only partially covered
409 409 bins, arising for example from different energy thresholds or ROIs (#1789).
410 410  
  411 +The cube response classes GCTACubeExposure, GCTACubePsf, GCTACubeEdisp and
  412 +GCTACubeBackground do no longer select on RoI and energy boundaries as this is now
  413 +handled by the weight array (#1789).
  414 +
411 415  
412 416 21. COMPTEL interface
413 417 ---------------------
... ...
README.md
  View file @b447670
1 1 GammaLib information
2 2 ====================
3   -* Version: 1.1.0.dev1 (20 June 2016)
  3 +* Version: 1.1.0.dev1 (21 June 2016)
4 4 * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
5 5  
6 6 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @b447670
... ... @@ -339,20 +339,6 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
339 339 *log << " in background cube computation." << std::endl;
340 340 }
341 341  
342   - // Extract region of interest from CTA observation
343   - GCTARoi roi = cta->roi();
344   -
345   - // Extract energy boundaries from CTA observation
346   - GEbounds obs_ebounds = cta->ebounds();
347   -
348   - // Check for RoI sanity
349   - if (!roi.is_valid()) {
350   - std::string msg = "No RoI information found in input observation "
351   - "\""+cta->name()+"\". Run ctselect to specify "
352   - "an RoI for this observation";
353   - throw GException::invalid_value(G_FILL, msg);
354   - }
355   -
356 342 // Set GTI of actual observations as the GTI of the event cube
357 343 eventcube.gti(cta->gti());
358 344  
... ... @@ -365,31 +351,20 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
365 351 // Get event bin
366 352 GCTAEventBin* bin = eventcube[i];
367 353  
368   - // Skip if energy is not contained within the energy boundaries of
369   - // the observation
370   - if (!obs_ebounds.contains(bin->energy())) {
371   - continue;
372   - }
373   -
374   - // Continue only if binned in contained in ROI
375   - if (roi.contains(*bin)) {
376   -
377   - // Compute model value for event bin. The model value is
378   - // given in counts/MeV/s/sr.
379   - double model = obs.models().eval(*bin, *cta);
380   -
381   - // Multiply by livetime to get the correct weighting for
382   - // each observation. We divide by the total livetime later
383   - // to get the background model in units of counts/MeV/s/sr.
384   - model *= livetime;
  354 + // Compute model value for event bin. The model value is
  355 + // given in counts/MeV/s/sr.
  356 + double model = obs.models().eval(*bin, *cta);
385 357  
386   - // Add existing number of counts
387   - model += bin->counts();
  358 + // Multiply by livetime to get the correct weighting for
  359 + // each observation. We divide by the total livetime later
  360 + // to get the background model in units of counts/MeV/s/sr.
  361 + model *= livetime;
388 362  
389   - // Store cumulated value (units: counts/MeV/sr)
390   - bin->counts(model);
  363 + // Add existing number of counts
  364 + model += bin->counts();
391 365  
392   - } // endif: bin was contained in RoI
  366 + // Store cumulated value (units: counts/MeV/sr)
  367 + bin->counts(model);
393 368  
394 369 } // endfor: looped over all bins
395 370  
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @b447670
... ... @@ -428,16 +428,6 @@ void GCTACubeEdisp::set(const GCTAObservation&amp; obs)
428 428 "valid event list.";
429 429 throw GException::invalid_value(G_SET, msg);
430 430 }
431   - const GCTARoi& roi = list->roi();
432   -
433   - // Check for RoI sanity
434   - if (!roi.is_valid()) {
435   - std::string msg = "No valid Region of Interest found in CTA "
436   - "observation \""+obs.name()+"\". Run ctselect "
437   - "to set the Region of Interest for this "
438   - "observation.";
439   - throw GException::invalid_value(G_SET, msg);
440   - }
441 431  
442 432 // Get CTA response
443 433 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
... ... @@ -469,46 +459,41 @@ void GCTACubeEdisp::set(const GCTAObservation&amp; obs)
469 459 // Get pixel sky direction
470 460 GSkyDir dir = m_cube.inx2dir(pixel);
471 461  
472   - // Continue only if pixel is within RoI
473   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
  462 + // Compute theta angle with respect to pointing direction
  463 + // in radians
  464 + double theta = pnt.dist(dir);
474 465  
475   - // Compute theta angle with respect to pointing direction
476   - // in radians
477   - double theta = pnt.dist(dir);
  466 + // Loop over all energy dispersion cube true energies
  467 + for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin){
478 468  
479   - // Loop over all energy dispersion cube true energies
480   - for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin){
  469 + // Get log10 of true energy in TeV
  470 + double logEsrc = m_elogmeans[iebin];
481 471  
482   - // Get log10 of true energy in TeV
483   - double logEsrc = m_elogmeans[iebin];
484   -
485   - // Loop over migration bins
486   - for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
487   -
488   - // Get migration for this migration bin
489   - double migra = m_migras[imigra];
  472 + // Loop over migration bins
  473 + for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
490 474  
491   - // Skip migration bin if it's not positive
492   - if (migra <= 0.0) {
493   - continue;
494   - }
  475 + // Get migration for this migration bin
  476 + double migra = m_migras[imigra];
495 477  
496   - // Compute log10 of reconstructed energy in TeV
497   - double logEobs = std::log10(migra) - logEsrc;
  478 + // Skip migration bin if it's not positive
  479 + if (migra <= 0.0) {
  480 + continue;
  481 + }
498 482  
499   - // Set map index
500   - int imap = offset(imigra, iebin);
  483 + // Compute log10 of reconstructed energy in TeV
  484 + double logEobs = std::log10(migra) - logEsrc;
501 485  
502   - // Set energy dispersion cube value
503   - m_cube(pixel, imap) = (*edisp)(logEobs, logEsrc,
504   - theta, 0.0,
505   - 0.0, 0.0);
  486 + // Set map index
  487 + int imap = offset(imigra, iebin);
506 488  
507   - } // endfor: looped over migration bins
  489 + // Set energy dispersion cube value
  490 + m_cube(pixel, imap) = (*edisp)(logEobs, logEsrc,
  491 + theta, 0.0,
  492 + 0.0, 0.0);
508 493  
509   - } // endfor: looped over energy bins
  494 + } // endfor: looped over migration bins
510 495  
511   - } // endif: pixel was within RoI
  496 + } // endfor: looped over energy bins
512 497  
513 498 } // endfor: looped over all pixels
514 499  
... ... @@ -576,23 +561,6 @@ void GCTACubeEdisp::fill(const GObservations&amp; obs, GLog* log)
576 561 continue;
577 562 }
578 563  
579   - // Extract region of interest from CTA observation
580   - GCTARoi roi = cta->roi();
581   -
582   - // Extract energy boundaries from CTA observation
583   - GEbounds obs_ebounds = cta->ebounds();
584   - GEnergy obs_emin = obs_ebounds.emin() - g_energy_margin;
585   - GEnergy obs_emax = obs_ebounds.emax() + g_energy_margin;
586   -
587   - // Check for RoI sanity
588   - if (!roi.is_valid()) {
589   - std::string msg = "No valid Region of Interest found in CTA "
590   - "observation \""+cta->name()+"\". Run ctselect "
591   - "to set the Region of Interest for this "
592   - "observation.";
593   - throw GException::invalid_value(G_FILL, msg);
594   - }
595   -
596 564 // Get CTA response
597 565 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
598 566 if (rsp == NULL) {
... ... @@ -632,60 +600,47 @@ void GCTACubeEdisp::fill(const GObservations&amp; obs, GLog* log)
632 600 // Get pixel sky direction
633 601 GSkyDir dir = m_cube.inx2dir(pixel);
634 602  
635   - // Continue only if pixel is within RoI
636   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
637   -
638   - // Compute theta angle with respect to pointing
639   - // direction in radians
640   - double theta = pnt.dist(dir);
641   -
642   - // Loop over all energy bins
643   - for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin) {
644   -
645   - // Skip if energy is not within observation energy
646   - // boundaries
647   - if (m_energies[iebin] < obs_emin ||
648   - m_energies[iebin] > obs_emax) {
649   - continue;
650   - }
651   -
652   - // Get log10 of true energy in TeV
653   - double logEsrc = m_elogmeans[iebin];
  603 + // Compute theta angle with respect to pointing direction in radians
  604 + double theta = pnt.dist(dir);
654 605  
655   - // Compute exposure weight
656   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
657   - cta->livetime();
  606 + // Loop over all energy bins
  607 + for (int iebin = 0; iebin < m_elogmeans.size(); ++iebin) {
658 608  
659   - // Accumulate weights
660   - exposure(pixel, iebin) += weight;
  609 + // Get log10 of true energy in TeV
  610 + double logEsrc = m_elogmeans[iebin];
661 611  
662   - // Loop over delta values
663   - for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
  612 + // Compute exposure weight
  613 + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
  614 + cta->livetime();
664 615  
665   - // Compute delta in radians
666   - double migra = m_migras[imigra];
  616 + // Accumulate weights
  617 + exposure(pixel, iebin) += weight;
667 618  
668   - // Skip migra bin if zero
669   - if (migra <= 0.0) {
670   - continue;
671   - }
  619 + // Loop over delta values
  620 + for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
672 621  
673   - // Compute log10 of reconstructed energy in TeV
674   - double logEobs = std::log10(migra) - logEsrc;
  622 + // Compute delta in radians
  623 + double migra = m_migras[imigra];
675 624  
676   - // Set map index
677   - int imap = offset(imigra, iebin);
  625 + // Skip migra bin if zero
  626 + if (migra <= 0.0) {
  627 + continue;
  628 + }
678 629  
679   - // Add energy dispersion cube value
680   - m_cube(pixel, imap) += (*edisp)(logEobs, logEsrc,
681   - theta, 0.0,
682   - 0.0, 0.0) * weight;
  630 + // Compute log10 of reconstructed energy in TeV
  631 + double logEobs = std::log10(migra) - logEsrc;
683 632  
684   - } // endfor: looped over migration bins
  633 + // Set map index
  634 + int imap = offset(imigra, iebin);
685 635  
686   - } // endfor: looped over energy bins
  636 + // Add energy dispersion cube value
  637 + m_cube(pixel, imap) += (*edisp)(logEobs, logEsrc,
  638 + theta, 0.0,
  639 + 0.0, 0.0) * weight;
687 640  
688   - } // endif: pixel was within RoI
  641 + } // endfor: looped over migration bins
  642 +
  643 + } // endfor: looped over energy bins
689 644  
690 645 } // endfor: looped over all pixels
691 646  
... ...
inst/cta/src/GCTACubeExposure.cpp
  View file @b447670
... ... @@ -300,18 +300,17 @@ GCTACubeExposure* GCTACubeExposure::clone(void) const
300 300  
301 301  
302 302 /***********************************************************************//**
303   - * @brief Set exposure cube from one CTA observation
  303 + * @brief Set exposure cube for one CTA observation
304 304 *
305 305 * @param[in] obs CTA observation.
306 306 *
307   - * Set the exposure cube from a single CTA observations. The cube pixel
308   - * values are computed as product of the effective area and the livetime.
  307 + * Set the exposure cube for a single CTA observation. The cube pixel values
  308 + * are computed as product of the effective area and the livetime.
309 309 *
310 310 * @todo: Throw an exception if response is not valid
311 311 ***************************************************************************/
312 312 void GCTACubeExposure::set(const GCTAObservation& obs)
313 313 {
314   -
315 314 // Only continue if we have an unbinned observation
316 315 if (obs.eventtype() == "EventList") {
317 316  
... ... @@ -320,17 +319,6 @@ void GCTACubeExposure::set(const GCTAObservation&amp; obs)
320 319 m_livetime = 0.0;
321 320 m_cube = 0.0;
322 321  
323   - // Extract region of interest from CTA observation
324   - GCTARoi roi = obs.roi();
325   -
326   - // Check for RoI sanity
327   - if (!roi.is_valid()) {
328   - std::string msg = "No RoI information found in input observation "
329   - "\""+obs.name()+"\". Run ctselect to specify "
330   - "an RoI for this observation";
331   - throw GException::invalid_value(G_SET, msg);
332   - }
333   -
334 322 // Get references on CTA response and pointing direction
335 323 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
336 324 const GSkyDir& pnt = obs.pointing().dir();
... ... @@ -344,26 +332,21 @@ void GCTACubeExposure::set(const GCTAObservation&amp; obs)
344 332 // Get pixel sky direction
345 333 GSkyDir dir = m_cube.inx2dir(pixel);
346 334  
347   - // Continue only if pixel is within RoI
348   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
349   -
350   - // Compute theta angle with respect to pointing direction
351   - // in radians
352   - double theta = pnt.dist(dir);
353   -
354   - // Loop over all exposure cube energies
355   - for (int iebin = 0; iebin < m_energies.size(); ++iebin){
  335 + // Compute theta angle with respect to pointing direction
  336 + // in radians
  337 + double theta = pnt.dist(dir);
356 338  
357   - // Get logE/TeV
358   - double logE = m_energies[iebin].log10TeV();
  339 + // Loop over all exposure cube energies
  340 + for (int iebin = 0; iebin < m_energies.size(); ++iebin){
359 341  
360   - // Set exposure cube (effective area * livetime)
361   - m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
362   - obs.livetime();
  342 + // Get logE/TeV
  343 + double logE = m_energies[iebin].log10TeV();
363 344  
364   - } // endfor: looped over energy bins
  345 + // Set exposure cube (effective area * livetime)
  346 + m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  347 + obs.livetime();
365 348  
366   - } // endif: pixel was within RoI
  349 + } // endfor: looped over energy bins
367 350  
368 351 } // endfor: looped over all pixels
369 352  
... ... @@ -430,20 +413,6 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
430 413 continue;
431 414 }
432 415  
433   - // Extract region of interest from CTA observation
434   - GCTARoi roi = cta->roi();
435   -
436   - // Extract energy boundaries from CTA observation
437   - GEbounds obs_ebounds = cta->ebounds();
438   -
439   - // Check for RoI sanity
440   - if (!roi.is_valid()) {
441   - std::string msg = "No RoI information found in input observation "
442   - "\""+cta->name()+"\". Run ctselect to specify "
443   - "an RoI for this observation";
444   - throw GException::invalid_value(G_FILL, msg);
445   - }
446   -
447 416 // Get references on CTA response and pointing direction
448 417 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
449 418 const GSkyDir& pnt = cta->pointing().dir();
... ... @@ -476,36 +445,20 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
476 445 // Get pixel sky direction
477 446 GSkyDir dir = m_cube.inx2dir(pixel);
478 447  
479   - // Continue only if pixel is within RoI
480   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
481   -
482   - // Compute theta angle with respect to pointing
483   - // direction in radians
484   - double theta = pnt.dist(dir);
  448 + // Compute theta angle with respect to pointing direction in radians
  449 + double theta = pnt.dist(dir);
485 450  
486   - // Loop over all exposure cube energy bins
487   - for (int iebin = 0; iebin < m_energies.size(); ++iebin){
488   -
489   - // Only add exposure if energy is inside observation energy
490   - // boundaries
491   - // TODO: The exposure cube energies are true energies, while
492   - // the observation energy boundaries are reconstructed
493   - // energies. In principle one would need to use the
494   - // energy dispersion information to do this correctly.
495   - if (obs_ebounds.contains(m_energies[iebin])) {
496   -
497   - // Get logE/TeV
498   - double logE = m_energies[iebin].log10TeV();
  451 + // Loop over all exposure cube energy bins
  452 + for (int iebin = 0; iebin < m_energies.size(); ++iebin){
499 453  
500   - // Add to exposure cube (effective area * livetime)
501   - m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
502   - cta->livetime();
  454 + // Get logE/TeV
  455 + double logE = m_energies[iebin].log10TeV();
503 456  
504   - } // endif: energy bin was inside observation energy boundaries
505   -
506   - } // endfor: looped over energy bins
  457 + // Add to exposure cube (effective area * livetime)
  458 + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  459 + cta->livetime();
507 460  
508   - } // endif: pixel within RoI
  461 + } // endfor: looped over energy bins
509 462  
510 463 } // endfor: looped over all pixels
511 464  
... ...
inst/cta/src/GCTACubePsf.cpp
  View file @b447670
... ... @@ -51,7 +51,6 @@
51 51 /* __ Debug definitions __________________________________________________ */
52 52  
53 53 /* __ Constants __________________________________________________________ */
54   -const GEnergy g_energy_margin(1.0e-12, "TeV");
55 54  
56 55  
57 56 /*==========================================================================
... ... @@ -378,25 +377,6 @@ void GCTACubePsf::set(const GCTAObservation&amp; obs)
378 377 // Only continue if we have an unbinned observation
379 378 if (obs.eventtype() == "EventList") {
380 379  
381   - // Extract region of interest from CTA observation
382   - const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
383   - if (list == NULL) {
384   - std::string msg = "CTA Observation does not contain an event "
385   - "list. Event list information is needed to "
386   - "retrieve the Region of Interest for each "
387   - "CTA observation.";
388   - throw GException::invalid_value(G_SET, msg);
389   - }
390   - const GCTARoi& roi = list->roi();
391   -
392   - // Check for RoI sanity
393   - if (!roi.is_valid()) {
394   - std::string msg = "No RoI information found in input observation "
395   - "\""+obs.name()+"\". Run ctselect to specify "
396   - "an RoI for this observation";
397   - throw GException::invalid_value(G_SET, msg);
398   - }
399   -
400 380 // Get references on CTA response and pointing direction
401 381 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
402 382 const GSkyDir& pnt = obs.pointing().dir();
... ... @@ -410,36 +390,31 @@ void GCTACubePsf::set(const GCTAObservation&amp; obs)
410 390 // Get pixel sky direction
411 391 GSkyDir dir = m_cube.inx2dir(pixel);
412 392  
413   - // Continue only if pixel is within RoI
414   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
415   -
416   - // Compute theta angle with respect to pointing direction
417   - // in radians
418   - double theta = pnt.dist(dir);
419   -
420   - // Loop over all exposure cube energy bins
421   - for (int iebin = 0; iebin < m_energies.size(); ++iebin){
  393 + // Compute theta angle with respect to pointing direction
  394 + // in radians
  395 + double theta = pnt.dist(dir);
422 396  
423   - // Get logE/TeV
424   - double logE = m_energies[iebin].log10TeV();
  397 + // Loop over all exposure cube energy bins
  398 + for (int iebin = 0; iebin < m_energies.size(); ++iebin){
425 399  
426   - // Loop over delta values
427   - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
428   -
429   - // Compute delta in radians
430   - double delta = m_deltas[idelta] * gammalib::deg2rad;
  400 + // Get logE/TeV
  401 + double logE = m_energies[iebin].log10TeV();
431 402  
432   - // Set map index
433   - int imap = offset(idelta, iebin);
  403 + // Loop over delta values
  404 + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
434 405  
435   - // Set PSF cube
436   - m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE);
  406 + // Compute delta in radians
  407 + double delta = m_deltas[idelta] * gammalib::deg2rad;
437 408  
438   - } // endfor: looped over delta bins
  409 + // Set map index
  410 + int imap = offset(idelta, iebin);
439 411  
440   - } // endfor: looped over energy bins
  412 + // Set PSF cube
  413 + m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE);
441 414  
442   - } // endif: pixel was within RoI
  415 + } // endfor: looped over delta bins
  416 +
  417 + } // endfor: looped over energy bins
443 418  
444 419 } // endfor: looped over all pixels
445 420  
... ... @@ -501,20 +476,6 @@ void GCTACubePsf::fill(const GObservations&amp; obs, GLog* log)
501 476 continue;
502 477 }
503 478  
504   - // Extract region of interest from CTA observation
505   - GCTARoi roi = cta->roi();
506   -
507   - // Extract energy boundaries from CTA observation
508   - GEbounds obs_ebounds = cta->ebounds();
509   -
510   - // Check for RoI sanity
511   - if (!roi.is_valid()) {
512   - std::string msg = "No RoI information found in input observation "
513   - "\""+cta->name()+"\". Run ctselect to specify "
514   - "an RoI for this observation";
515   - throw GException::invalid_value(G_FILL, msg);
516   - }
517   -
518 479 // Get references on CTA response and pointing direction
519 480 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
520 481 const GSkyDir& pnt = cta->pointing().dir();
... ... @@ -547,54 +508,38 @@ void GCTACubePsf::fill(const GObservations&amp; obs, GLog* log)
547 508 // Get pixel sky direction
548 509 GSkyDir dir = m_cube.inx2dir(pixel);
549 510  
550   - // Continue only if pixel is within RoI
551   - if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
552   -
553   - // Compute theta angle with respect to pointing
554   - // direction in radians
555   - double theta = pnt.dist(dir);
556   -
557   - // Loop over all energy bins
558   - for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
559   -
560   - // Only add exposure if energy is inside observation energy
561   - // boundaries
562   - // TODO: The exposure cube energies are true energies, while
563   - // the observation energy boundaries are reconstructed
564   - // energies. In principle one would need to use the
565   - // energy dispersion information to do this correctly.
566   - if (!obs_ebounds.contains(m_energies[iebin])) {
567   - continue;
568   - }
  511 + // Compute theta angle with respect to pointing direction in radians
  512 + double theta = pnt.dist(dir);
569 513  
570   - // Get logE/TeV
571   - double logE = m_energies[iebin].log10TeV();
  514 + // Loop over all energy bins
  515 + for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
572 516  
573   - // Compute exposure weight
574   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
575   - cta->livetime();
  517 + // Get logE/TeV
  518 + double logE = m_energies[iebin].log10TeV();
576 519  
577   - // Accumulate weights
578   - exposure(pixel, iebin) += weight;
  520 + // Compute exposure weight
  521 + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
  522 + cta->livetime();
579 523  
580   - // Loop over delta values
581   - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
  524 + // Accumulate weights
  525 + exposure(pixel, iebin) += weight;
582 526  
583   - // Compute delta in radians
584   - double delta = m_deltas[idelta] * gammalib::deg2rad;
  527 + // Loop over delta values
  528 + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
585 529  
586   - // Set map index
587   - int imap = offset(idelta, iebin);
  530 + // Compute delta in radians
  531 + double delta = m_deltas[idelta] * gammalib::deg2rad;
588 532  
589   - // Add on PSF cube
590   - m_cube(pixel, imap) +=
591   - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
  533 + // Set map index
  534 + int imap = offset(idelta, iebin);
592 535  
593   - } // endfor: looped over delta bins
  536 + // Add on PSF cube
  537 + m_cube(pixel, imap) +=
  538 + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
594 539  
595   - } // endfor: looped over energy bins
  540 + } // endfor: looped over delta bins
596 541  
597   - } // endif: pixel was within RoI
  542 + } // endfor: looped over energy bins
598 543  
599 544 } // endfor: looped over all pixels
600 545  
... ...