Simon Bonnefoy / gammalib

Commit 6dca3a06ad8c20bb4d704b9ce2bd222f3b0726f6

Authored by Jürgen Knödlseder
1 parent 4a6eab84

Assume constant background in the absence of an IRF template in GCTAOnOffObserva… 

…tion::compute_bgd (#2683)
Inline Side-by-side
Showing 1 changed file with 196 additions and 168 deletions
inst/cta/src/GCTAOnOffObservation.cpp
  View file @6dca3a0
... ... @@ -293,8 +293,8 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations& obs) :
293 293 // Signal first On/Off observation
294 294 bool first = true;
295 295  
296   - // Initialise exposure
297   - double exposure = 0.0;
  296 + // Initialise total exposure
  297 + double total_exposure = 0.0;
298 298  
299 299 // Allocate observation energy range
300 300 GEnergy emin_obs;
... ... @@ -323,38 +323,35 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations& obs) :
323 323 GArf arf_stacked = this->arf_stacked(onoff->arf(), emin, emax);
324 324 GRmf rmf_stacked = this->rmf_stacked(onoff->rmf(), emin, emax);
325 325  
  326 + // Get exposure for observation
  327 + double exposure = onoff->on_spec().exposure();
  328 +
326 329 // If this is the first On/Off observation then store the data to
327 330 // initialise the data definition
328 331 if (first) {
329 332  
330   - // Store PHA, ARF and RMF
331   - m_on_spec = onoff->on_spec();
332   - m_off_spec = onoff->off_spec();
333   - m_arf = arf_stacked * onoff->on_spec().exposure();
334   - m_rmf = rmf_stacked;
335   - m_ontime = onoff->ontime();
336   - m_livetime = onoff->livetime();
337   - exposure = onoff->on_spec().exposure();
  333 + // Store PHA, ARF and RMF, set exposure, ontime and livetime
  334 + m_on_spec = onoff->on_spec();
  335 + m_off_spec = onoff->off_spec();
  336 + m_arf = arf_stacked * exposure;
  337 + m_rmf = rmf_stacked;
  338 + total_exposure = exposure;
  339 + m_ontime = onoff->ontime();
  340 + m_livetime = onoff->livetime();
338 341  
339 342 // Compute number of background events/MeV from BACKRESP column
340 343 // and store result intermediately into BACKRESP column of off
341 344 // spectum
342 345 std::vector<double>& backresp = m_off_spec["BACKRESP"];
343 346 for (int i = 0; i < backresp.size(); ++i) {
344   - backresp[i] *= onoff->on_spec().exposure();
  347 + backresp[i] *= exposure;
345 348 }
346 349  
347   - // Compute background scaling factor contribution. Handle the case
348   - // where there is no background response, which may happen if no
349   - // IRF template exists. In that case we do not weight with the
350   - // background response.
  350 + // Compute background scaling factor contribution
351 351 for (int i = 0; i < m_on_spec.size(); ++i) {
352 352 double background = onoff->off_spec()["BACKRESP"][i];
353   - double exposure = onoff->on_spec().exposure();
354 353 double alpha = onoff->on_spec().backscal(i);
355   - double scale = (background > 0.0)
356   - ? alpha * background * exposure
357   - : alpha * exposure;
  354 + double scale = alpha * background * exposure;
358 355 m_on_spec.backscal(i,scale);
359 356 }
360 357  
... ... @@ -410,44 +407,36 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations&amp; obs) :
410 407 m_on_spec += onoff->on_spec();
411 408 m_off_spec += onoff->off_spec();
412 409  
413   - // Compute background scaling factor contribution. Handle the case
414   - // where there is no background response, which may happen if no
415   - // IRF template exists. In that case we do not weight with the
416   - // background response.
  410 + // Compute background scaling factor contribution
417 411 for (int i = 0; i < m_on_spec.size(); ++i) {
418 412 double background = onoff->off_spec()["BACKRESP"][i];
419 413 double alpha = onoff->on_spec().backscal(i);
420   - double scale = m_on_spec.backscal(i) +
421   - ((background > 0.0)
422   - ? alpha * background * onoff->on_spec().exposure()
423   - : alpha * onoff->on_spec().exposure());
424   - m_on_spec.backscal(i,scale);
  414 + double scale = alpha * background * exposure;
  415 + m_on_spec.backscal(i, m_on_spec.backscal(i) + scale);
425 416 }
426 417  
427 418 // Add ARF
428   - m_arf += arf_stacked * onoff->on_spec().exposure();
  419 + m_arf += arf_stacked * exposure;
429 420  
430 421 // Add number of background events/MeV from BACKRESP column
431 422 const std::vector<double>& src = onoff->off_spec()["BACKRESP"];
432 423 std::vector<double>& dst = m_off_spec["BACKRESP"];
433 424 for (int i = 0; i < dst.size(); ++i) {
434   - dst[i] += src[i] * onoff->on_spec().exposure();
  425 + dst[i] += src[i] * exposure;
435 426 }
436 427  
437 428 // Add RMF
438 429 for (int itrue = 0; itrue < m_rmf.ntrue(); ++itrue) {
439   - double arf = arf_stacked[itrue] * onoff->on_spec().exposure();
  430 + double arf = arf_stacked[itrue] * exposure;
440 431 for (int imeasured = 0; imeasured < m_rmf.nmeasured(); ++imeasured) {
441 432 m_rmf(itrue,imeasured) += rmf_stacked(itrue,imeasured) * arf;
442 433 }
443 434 }
444 435  
445   - // Add exposure
446   - exposure += onoff->on_spec().exposure();
447   -
448   - // Add ontime and livetime
449   - m_ontime += onoff->ontime();
450   - m_livetime += onoff->livetime();
  436 + // Add exposure, ontime and livetime
  437 + total_exposure += exposure;
  438 + m_ontime += onoff->ontime();
  439 + m_livetime += onoff->livetime();
451 440  
452 441 // Update observation energy range
453 442 if (emin < emin_obs) {
... ... @@ -461,9 +450,7 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations&amp; obs) :
461 450  
462 451 } // endfor: looped over observations
463 452  
464   - // Compute stacked background scaling factor. Handle the case where there
465   - // is no background response, which may happen if no IRF template exists.
466   - // In that case we simply have to devide by the exposure.
  453 + // Compute stacked background scaling factor
467 454 for (int i = 0; i < m_on_spec.size(); ++i) {
468 455 double norm = m_off_spec["BACKRESP"][i];
469 456 if (norm > 0.0) {
... ... @@ -471,8 +458,7 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations&amp; obs) :
471 458 m_on_spec.backscal(i,scale);
472 459 }
473 460 else {
474   - double scale = m_on_spec.backscal(i) / exposure;
475   - m_on_spec.backscal(i,scale);
  461 + m_on_spec.backscal(i,0.0);
476 462 }
477 463 }
478 464  
... ... @@ -487,15 +473,15 @@ GCTAOnOffObservation::GCTAOnOffObservation(const GObservations&amp; obs) :
487 473 }
488 474  
489 475 // Compute ARF
490   - if (exposure > 0.0) {
491   - m_arf /= exposure;
  476 + if (total_exposure > 0.0) {
  477 + m_arf /= total_exposure;
492 478 }
493 479  
494 480 // Compute background events/MeV/sec and store them in BACKRESP column
495   - if (exposure > 0.0) {
  481 + if (total_exposure > 0.0) {
496 482 std::vector<double>& backresp = m_off_spec["BACKRESP"];
497 483 for (int i = 0; i < backresp.size(); ++i) {
498   - backresp[i] /= exposure;
  484 + backresp[i] /= total_exposure;
499 485 }
500 486 }
501 487  
... ... @@ -834,6 +820,122 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
834 820 }
835 821  
836 822  
  823 +/***********************************************************************
  824 + * @brief Compute predicted gamma-ray counts for given model
  825 + *
  826 + * @param[in] models Model container.
  827 + * @return Model PHA with predicted gamma-ray counts.
  828 + *
  829 + * Returns spectrum of predicted number of source events in the On
  830 + * regions.
  831 + ***********************************************************************/
  832 +GPha GCTAOnOffObservation::model_gamma(const GModels& models) const
  833 +{
  834 + // Get number of model parameters in model container
  835 + int npars = models.npars();
  836 +
  837 + // Create dummy gradient vectors to provide required argument to N_gamma
  838 + GVector dummy_grad(npars);
  839 +
  840 + // Initialise output spectrum
  841 + GEbounds ereco = m_on_spec.ebounds();
  842 + GPha gammas = GPha(ereco);
  843 +
  844 + // Loop over all energy bins
  845 + for (int i = 0; i < m_on_spec.size(); ++i) {
  846 + gammas.fill(ereco.emean(i), N_gamma(models, i, &dummy_grad));
  847 + }
  848 +
  849 + // Return predicted gamma-ray counts spectrum
  850 + return gammas;
  851 +}
  852 +
  853 +
  854 +/***********************************************************************
  855 + * @brief Compute predicted background counts PHA for given model
  856 + *
  857 + * @param[in] models Model container.
  858 + * @return Model background PHA.
  859 + *
  860 + * Returns spectrum of predicted number of background events in the Off
  861 + * regions. The computation method changed depending on the statistic used
  862 + * to match the spectrum used for likelihood fitting.
  863 + * To be noted: for the WSTAT statistic the model background depends on
  864 + * on the model adopted for gamma rays, that therefore need to be passed
  865 + * to this function.
  866 + ***********************************************************************/
  867 +GPha GCTAOnOffObservation::model_background(const GModels& models) const
  868 +{
  869 + // Get number of model parameters in model container
  870 + int npars = models.npars();
  871 +
  872 + // Create dummy gradient vectors to provide required argument to N_gamma
  873 + GVector dummy_grad(npars);
  874 +
  875 + // Initialise output spectrum
  876 + GEbounds ereco = m_on_spec.ebounds();
  877 + GPha bgds = GPha(ereco);
  878 +
  879 + // Extract statistic for this observation
  880 + std::string statistic = gammalib::toupper(this->statistic());
  881 +
  882 + // Loop over all energy bins
  883 + for (int i = 0; i < m_on_spec.size(); ++i) {
  884 +
  885 + // Initialise variable to store number of background counts
  886 + double nbgd = 0.0;
  887 +
  888 + // Choose background evaluation method based on fit statistics
  889 +
  890 + // CSTAT
  891 + if ((statistic == "POISSON") || (statistic == "CSTAT")) {
  892 + nbgd = N_bgd(models, i, &dummy_grad);
  893 + }
  894 +
  895 + // ... or Poisson statistic with measured background
  896 + else if (statistic == "WSTAT") {
  897 +
  898 + // Get number of On and Off counts
  899 + double non = m_on_spec[i];
  900 + double noff = m_off_spec[i];
  901 +
  902 + // Get background scaling
  903 + double alpha = m_on_spec.backscal(i);
  904 +
  905 + // Get number of gamma and background events (and corresponding
  906 + // spectral model gradients)
  907 + double ngam = N_gamma(models, i, &dummy_grad);
  908 +
  909 + // Initialise variables for likelihood computation
  910 + double nonpred = 0.0;
  911 + double dlogLdsky = 0.0;
  912 + double d2logLdsky2 = 0.0;
  913 +
  914 + // Perform likelihood profiling and derive likelihood value
  915 + // The number of background counts is updated to the profile value
  916 + double logL = wstat_value(non, noff, alpha, ngam, nonpred, nbgd,
  917 + dlogLdsky,d2logLdsky2);
  918 +
  919 + } // endelse: WSTAT
  920 +
  921 + // ... or unsupported
  922 + else {
  923 + std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
  924 + "Either specify \"POISSON\", \"CSTAT\" or "
  925 + "\"WSTAT\".";
  926 + throw GException::invalid_value(G_MODEL_BACKGROUND, msg);
  927 + }
  928 +
  929 + // Fill background spectrum
  930 + bgds.fill(ereco.emean(i), nbgd);
  931 +
  932 + } // endfor: looped over energy bins
  933 +
  934 + // Return model count spectrum
  935 + return bgds;
  936 +}
  937 +
  938 +
837 939 /***********************************************************************//**
838 940 * @brief Print On/Off observation information
839 941 *
... ... @@ -1470,13 +1572,15 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1470 1572 * \f$\Omega_{{\rm off},p}\f$ is the solid angle in units of \f${\rm sr}\f$
1471 1573 * of the pixel index \f$p\f$ in the Off region \f${\rm off}\f$.
1472 1574 *
1473   - * If @p use_irf_bkg is @c false, all `BACKRESP` elements will be set to
1474   - * zero, i.e.
  1575 + * If @p use_irf_bkg is @c false, `BACKRESP` will be computed using
1475 1576 *
1476 1577 * \f[
1477   - * {\tt BACKRESP}(E_{\rm reco}) = 0
  1578 + * {\tt BACKRESP}(E_{\rm reco}) = \sum_{\rm off} \sum_p \Omega_{{\rm off},p}
1478 1579 * \f]
1479 1580 *
  1581 + * which actually assumes that the background rate is constant over the
  1582 + * field of view.
  1583 + *
1480 1584 * @todo Integrate background rate over energy bin instead of computing the
1481 1585 * rate at the energy bin centre.
1482 1586 ***************************************************************************/
... ... @@ -1494,12 +1598,12 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1494 1598 // Initialise background rates to zero
1495 1599 std::vector<double> background(nreco, 0.0);
1496 1600  
  1601 + // Get CTA observation pointing direction
  1602 + GCTAPointing obspnt = obs.pointing();
  1603 +
1497 1604 // Continue only if the IRF background template should be used
1498 1605 if (use_irf_bkg) {
1499 1606  
1500   - // Get CTA observation pointing direction
1501   - GCTAPointing obspnt = obs.pointing();
1502   -
1503 1607 // Get CTA background
1504 1608 const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_COMPUTE_BGD, obs);
1505 1609  
... ... @@ -1545,6 +1649,46 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs,
1545 1649  
1546 1650 } // endif: IRF background template was used
1547 1651  
  1652 + // ... otherwise
  1653 + else {
  1654 +
  1655 + // Loop over regions
  1656 + for (int k = 0; k < off.size(); ++k) {
  1657 +
  1658 + // Get pointer on sky region map
  1659 + const GSkyRegionMap* off_map =
  1660 + static_cast<const GSkyRegionMap*>(off[k]);
  1661 +
  1662 + // Loop over pixels in Off region map and integrate
  1663 + // background rate
  1664 + for (int j = 0; j < off_map->nonzero_indices().size(); ++j) {
  1665 +
  1666 + // Get pixel index
  1667 + int pixidx = off_map->nonzero_indices()[j];
  1668 +
  1669 + // Get direction to pixel center
  1670 + GSkyDir pixdir = off_map->map().inx2dir(pixidx);
  1671 +
  1672 + // Translate sky direction into instrument direction
  1673 + GCTAInstDir pixinstdir = obspnt.instdir(pixdir);
  1674 +
  1675 + // Get solid angle subtended by this pixel
  1676 + double pixsolid = off_map->map().solidangle(pixidx);
  1677 +
  1678 + // Loop over energy bins
  1679 + for (int i = 0; i < nreco; ++i) {
  1680 +
  1681 + // Get background rate in events/s/MeV
  1682 + background[i] += pixsolid;
  1683 +
  1684 + } // endfor: looped over energy bins
  1685 +
  1686 + } // endfor: looped over all pixels in map
  1687 +
  1688 + } // endfor: looped over all regions
  1689 +
  1690 + } // endelse: no IRF background template was used
  1691 +
1548 1692 // Append background vector to Off spectrum
1549 1693 m_off_spec.append("BACKRESP", background);
1550 1694  
... ... @@ -2318,122 +2462,6 @@ double GCTAOnOffObservation::N_bgd(const GModels&amp; models,
2318 2462 }
2319 2463  
2320 2464  
2321   -/***********************************************************************
2322   - * @brief Compute predicted gamma-ray counts for given model
2323   - *
2324   - * @param[in] models Model container.
2325   - * @return Model PHA with predicted gamma-ray counts.
2326   - *
2327   - * Returns spectrum of predicted number of source events in the On
2328   - * regions.
2329   - ***********************************************************************/
2330   -GPha GCTAOnOffObservation::model_gamma(const GModels& models) const
2331   -{
2332   - // Get number of model parameters in model container
2333   - int npars = models.npars();
2334   -
2335   - // Create dummy gradient vectors to provide required argument to N_gamma
2336   - GVector dummy_grad(npars);
2337   -
2338   - // Initialise output spectrum
2339   - GEbounds ereco = m_on_spec.ebounds();
2340   - GPha gammas = GPha(ereco);
2341   -
2342   - // Loop over all energy bins
2343   - for (int i = 0; i < m_on_spec.size(); ++i) {
2344   - gammas.fill(ereco.emean(i), N_gamma(models, i, &dummy_grad));
2345   - }
2346   -
2347   - // Return predicted gamma-ray counts spectrum
2348   - return gammas;
2349   -}
2350   -
2351   -
2352   -/***********************************************************************
2353   - * @brief Compute predicted background counts PHA for given model
2354   - *
2355   - * @param[in] models Model container.
2356   - * @return Model background PHA.
2357   - *
2358   - * Returns spectrum of predicted number of background events in the Off
2359   - * regions. The computation method changed depending on the statistic used
2360   - * to match the spectrum used for likelihood fitting.
2361   - * To be noted: for the WSTAT statistic the model background depends on
2362   - * on the model adopted for gamma rays, that therefore need to be passed
2363   - * to this function.
2364   - ***********************************************************************/
2365   -GPha GCTAOnOffObservation::model_background(const GModels& models) const
2366   -{
2367   - // Get number of model parameters in model container
2368   - int npars = models.npars();
2369   -
2370   - // Create dummy gradient vectors to provide required argument to N_gamma
2371   - GVector dummy_grad(npars);
2372   -
2373   - // Initialise output spectrum
2374   - GEbounds ereco = m_on_spec.ebounds();
2375   - GPha bgds = GPha(ereco);
2376   -
2377   - // Extract statistic for this observation
2378   - std::string statistic = gammalib::toupper(this->statistic());
2379   -
2380   - // Loop over all energy bins
2381   - for (int i = 0; i < m_on_spec.size(); ++i) {
2382   -
2383   - // Initialise variable to store number of background counts
2384   - double nbgd = 0.0;
2385   -
2386   - // Choose background evaluation method based on fit statistics
2387   -
2388   - // CSTAT
2389   - if ((statistic == "POISSON") || (statistic == "CSTAT")) {
2390   - nbgd = N_bgd(models, i, &dummy_grad);
2391   - }
2392   -
2393   - // ... or Poisson statistic with measured background
2394   - else if (statistic == "WSTAT") {
2395   -
2396   - // Get number of On and Off counts
2397   - double non = m_on_spec[i];
2398   - double noff = m_off_spec[i];
2399   -
2400   - // Get background scaling
2401   - double alpha = m_on_spec.backscal(i);
2402   -
2403   - // Get number of gamma and background events (and corresponding
2404   - // spectral model gradients)
2405   - double ngam = N_gamma(models, i, &dummy_grad);
2406   -
2407   - // Initialise variables for likelihood computation
2408   - double nonpred = 0.0;
2409   - double dlogLdsky = 0.0;
2410   - double d2logLdsky2 = 0.0;
2411   -
2412   - // Perform likelihood profiling and derive likelihood value
2413   - // The number of background counts is updated to the profile value
2414   - double logL = wstat_value(non, noff, alpha, ngam, nonpred, nbgd,
2415   - dlogLdsky,d2logLdsky2);
2416   -
2417   - } // endelse: WSTAT
2418   -
2419   - // ... or unsupported
2420   - else {
2421   - std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
2422   - "Either specify \"POISSON\", \"CSTAT\" or "
2423   - "\"WSTAT\".";
2424   - throw GException::invalid_value(G_MODEL_BACKGROUND, msg);
2425   - }
2426   -
2427   - // Fill background spectrum
2428   - bgds.fill(ereco.emean(i), nbgd);
2429   -
2430   - } // endfor: looped over energy bins
2431   -
2432   - // Return model count spectrum
2433   - return bgds;
2434   -}
2435   -
2436   -
2437 2465 /***********************************************************************//**
2438 2466 * @brief Evaluate log-likelihood function for On/Off analysis in the
2439 2467 * case of Poisson signal with modeled Poisson background
... ...