Marianne Lemoine / ctools

Commit 28a6964e3fc33a74658b67a67704d853517e6de4

Authored by Jürgen Knödlseder
1 parent 8af322b8

Compute model values only for bins that are within the energy boundaries of the observation.

Inline Side-by-side
Showing 1 changed file with 31 additions and 10 deletions
src/ctmodel/ctmodel.cpp
  View file @28a6964
... ... @@ -504,9 +504,8 @@ void ctmodel::get_parameters(void)
504 504 m_deadc = (*this)["deadc"].real();
505 505 m_tmin = (*this)["tmin"].real();
506 506 m_tmax = (*this)["tmax"].real();
507   - if (m_obs.models().size() == 0) {
508   - m_srcmdl = (*this)["srcmdl"].filename();
509   - }
  507 + m_emin = (*this)["emin"].real();
  508 + m_emax = (*this)["emax"].real();
510 509 }
511 510  
512 511 }
... ... @@ -592,13 +591,20 @@ void ctmodel::setup_obs(void)
592 591 GTime tmax(m_tmax);
593 592 gti.append(tmin, tmax);
594 593  
  594 + // Setup energy range covered by model
  595 + GEnergy emin(m_emin, "TeV");
  596 + GEnergy emax(m_emax, "TeV");
  597 + GEbounds ebounds(1, emin, emax);
  598 +
595 599 // Allocate CTA observation and empty event list
596 600 GCTAObservation obs;
597 601 GCTAEventList list;
598 602  
599   - // Set event list GTI and append event list to observation
600   - // (we need in fact only the GTI for the computations)
  603 + // Set event list GTI and energy boundaries and append event list
  604 + // to observation (we need in fact only the GTI for the
  605 + // computations)
601 606 list.gti(gti);
  607 + list.ebounds(ebounds);
602 608 obs.events(list);
603 609  
604 610 // Set CTA observation attributes
... ... @@ -755,8 +761,13 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
755 761 // Continue only if observation pointer is valid
756 762 if (obs != NULL) {
757 763  
  764 + // Get energy boundaries and GTI references for observation
  765 + const GEbounds& ebounds = obs->events()->ebounds();
  766 + const GGti& gti = obs->events()->gti();
  767 +
758 768 // Initialise statistics
759   - double sum = 0.0;
  769 + double sum = 0.0;
  770 + int num_outside_ebds = 0;
760 771  
761 772 // Setup cube GTIs for this observation
762 773 m_cube.gti(obs->events()->gti());
... ... @@ -764,9 +775,17 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
764 775 // Loop over all cube bins
765 776 for (int i = 0; i < m_cube.size(); ++i) {
766 777  
767   - // Get cube bin and actual value
768   - GCTAEventBin* bin = m_cube[i];
769   - double value = bin->counts();
  778 + // Get cube bin
  779 + GCTAEventBin* bin = m_cube[i];
  780 +
  781 + // Skip bin if it is outside the energy range of the observation
  782 + if (!ebounds.contains(bin->energy())) {
  783 + num_outside_ebds++;
  784 + continue;
  785 + }
  786 +
  787 + // Get actual bin value
  788 + double value = bin->counts();
770 789  
771 790 // Compute model value for cube bin
772 791 double model = m_obs.models().eval(*bin, *obs) * bin->size();
... ... @@ -781,7 +800,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
781 800 } // endfor: looped over all cube bins
782 801  
783 802 // Append GTIs of observation to list of GTIs
784   - m_gti.extend(obs->events()->gti());
  803 + m_gti.extend(gti);
785 804  
786 805 // Update GTIs
787 806 m_cube.gti(m_gti);
... ... @@ -790,6 +809,8 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
790 809 if (logTerse()) {
791 810 log << gammalib::parformat("Model events in cube");
792 811 log << sum << std::endl;
  812 + log << gammalib::parformat("Cube bins outside energy range");
  813 + log << num_outside_ebds << std::endl;
793 814 }
794 815  
795 816 // Log cube
... ...