Justine Devin / gammalib

Commit 91fe979af5eb5ad652c27adad5fba3ef12d54195

Authored by Jürgen Knödlseder
1 parent 26dcc0cd

GObservation::npred_spec now skips and invalid energy ranges. Furthermore, energ… 

…y intergation is done over all energy ranges in an energy boundary object (see #1450).

Updated at the same time the Doxygen information.
Inline Side-by-side
Showing 2 changed files with 144 additions and 105 deletions
src/model/GModelSky.cpp
  View file @91fe979
... ... @@ -473,7 +473,6 @@ GVector GModelSky::gradients(const GPhoton& photon)
473 473 double GModelSky::eval(const GEvent& event, const GObservation& obs) const
474 474 {
475 475 // Evaluate function
476   - //double value = integrate_time(event, obs, false);
477 476 double value = obs.response()->convolve(*this, event, obs, false);
478 477  
479 478 // Return
... ... @@ -495,11 +494,10 @@ double GModelSky::eval(const GEvent& event, const GObservation& obs) const
495 494 * While the value of the sky model is returned by the method, the parameter
496 495 * gradients are set as GModelPar members.
497 496 ***************************************************************************/
498   -double GModelSky::eval_gradients(const GEvent& event,
  497 +double GModelSky::eval_gradients(const GEvent& event,
499 498 const GObservation& obs) const
500 499 {
501 500 // Evaluate function
502   - //double value = integrate_time(event, obs, true);
503 501 double value = obs.response()->convolve(*this, event, obs, true);
504 502  
505 503 // Return
... ... @@ -515,33 +513,35 @@ double GModelSky::eval_gradients(const GEvent& event,
515 513 * @param[in] obs Observation.
516 514 * @return Spatially integrated sky model.
517 515 *
518   - * @exception GException::no_response
519   - * No valid instrument response function defined.
520   - *
521 516 * Computes
522   - * \f[N"_{\rm pred} = \int_{\rm ROI}
523   - * S(\vec{p}, E, t) PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t) \,
524   - * {\rm d}\vec{p'}\f]
525   - * where
526   - * \f$S(\vec{p}, E, t)\f$ is the source model,
527   - * \f$PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t)\f$ is the point
528   - * spread function,
529   - * \f$\vec{p'}\f$ is the measured photon direction,
  517 + *
  518 + * \f[
  519 + * N_{\rm pred}(E',t') = \int_{\rm ROI}
  520 + * P(p',E',t') \, dp' \, dE'
  521 + * \f]
  522 + *
  523 + * of the event probability
  524 + *
  525 + * \f[
  526 + * P(p',E',t') = \int \int \int
  527 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  528 + * \f]
  529 + *
  530 + * where
  531 + * \f$S(p,E,t)\f$ is the source model,
  532 + * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  533 + * \f$p'\f$ is the measured photon direction,
530 534 * \f$E'\f$ is the measured photon energy,
531 535 * \f$t'\f$ is the measured photon arrival time,
532   - * \f$\vec{p}\f$ is the true photon arrival direction,
533   - * \f$E\f$ is the true photon energy,
534   - * \f$t\f$ is the true photon arrival time, and
535   - * \f$d\f$ is the instrument pointing.
536   - *
537   - * \f${\rm ROI}\f$ is the region of interest that is stored in the
538   - * GObservation::m_roi member. The integration over the ROI is performed
539   - * by the GResponse::npred() method.
  536 + * \f$p\f$ is the true photon arrival direction,
  537 + * \f$E\f$ is the true photon energy, and
  538 + * \f$t\f$ is the true photon arrival time.
540 539 *
541   - * The method takes care of any instrument dependent scale factors. These
542   - * scale factors will be applied to the predicted number of model counts.
  540 + * The method calls the GResponse::nroi() method that does the relevant
  541 + * integration.
543 542 ***************************************************************************/
544   -double GModelSky::npred(const GEnergy& obsEng, const GTime& obsTime,
  543 +double GModelSky::npred(const GEnergy& obsEng,
  544 + const GTime& obsTime,
545 545 const GObservation& obs) const
546 546 {
547 547 // Initialise result
... ...
src/obs/GObservation.cpp
  View file @91fe979
1 1 /***************************************************************************
2 2 * GObservation.cpp - Abstract observation base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -253,8 +253,9 @@ double GObservation::likelihood(const GModels& models,
253 253 * and observation identifiers matches those of the observation. Models that
254 254 * do not match will be skipped.
255 255 ***************************************************************************/
256   -double GObservation::model(const GModels& models, const GEvent& event,
257   - GVector* gradient) const
  256 +double GObservation::model(const GModels& models,
  257 + const GEvent& event,
  258 + GVector* gradient) const
258 259 {
259 260 // Initialise method variables
260 261 double model = 0.0; // Reset model value
... ... @@ -430,27 +431,35 @@ double GObservation::npred(const GModels& models, GVector* gradient) const
430 431 * Good Time Interval is invalid.
431 432 *
432 433 * Computes
433   - * \f[N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
434   - * S(\vec{p}, E, t) PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t) \,
435   - * {\rm d}\vec{p'} {\rm d}E' {\rm d}t'\f]
436   - * where
437   - * \f$S(\vec{p}, E, t)\f$ is the source model,
438   - * \f$PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t)\f$ is the point
439   - * spread function,
440   - * \f$\vec{p'}\f$ is the measured photon direction,
  434 + *
  435 + * \f[
  436 + * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
  437 + * P(p',E',t') \, dp' \, dE' \, dt'
  438 + * \f]
  439 + *
  440 + * of the event probability
  441 + *
  442 + * \f[
  443 + * P(p',E',t') = \int \int \int
  444 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  445 + * \f]
  446 + *
  447 + * where
  448 + * \f$S(p,E,t)\f$ is the source model,
  449 + * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  450 + * \f$p'\f$ is the measured photon direction,
441 451 * \f$E'\f$ is the measured photon energy,
442 452 * \f$t'\f$ is the measured photon arrival time,
443   - * \f$\vec{p}\f$ is the true photon arrival direction,
444   - * \f$E\f$ is the true photon energy,
445   - * \f$t\f$ is the true photon arrival time, and
446   - * \f$d\f$ is the instrument pointing.
447   - *
448   - * \f${\rm GTI}\f$ are the Good Time Intervals that are stored in the
449   - * GObservation::m_gti member.
450   - * Note that MET is used for the time integration interval. This, however,
451   - * is no specialisation since npred_grad_kern_spec::eval() converts the
452   - * argument back in a GTime object by assuming that the argument is in MET,
453   - * hence the correct time system will be used at the end by the method.
  453 + * \f$p\f$ is the true photon arrival direction,
  454 + * \f$E\f$ is the true photon energy, and
  455 + * \f$t\f$ is the true photon arrival time.
  456 + *
  457 + * This method performs the time integration over the Good Time Intervals
  458 + * \f${\rm GTI}\f$. Note that MET is used for the time integration interval.
  459 + * This, however, is no specialisation since npred_grad_kern_spec::eval()
  460 + * converts the argument back in a GTime object by assuming that the argument
  461 + * is in MET, hence the correct time system will be used at the end by the
  462 + * method.
454 463 ***************************************************************************/
455 464 double GObservation::npred(const GModel& model) const
456 465 {
... ... @@ -610,24 +619,36 @@ double GObservation::model_grad(const GModel& model,
610 619 * @param[in] par Model parameter.
611 620 *
612 621 * Computes
613   - * \f[\frac{{\rm d} N_{\rm pred}}{{\rm d} a_i}\f]
  622 + *
  623 + * \f[
  624 + * \frac{{\rm d} N_{\rm pred}}{{\rm d} a_i}
  625 + * \f]
  626 + *
614 627 * where
615   - * \f[N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
616   - * S(\vec{p}, E, t) PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t) \,
617   - * {\rm d}\vec{p'} {\rm d}E' {\rm d}t'\f]
  628 + *
  629 + * \f[
  630 + * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
  631 + * P(p',E',t') \, dp' \, dE' \, dt'
  632 + * \f]
  633 + *
  634 + * is the integral of the event probability
  635 + *
  636 + * \f[
  637 + * P(p',E',t') = \int \int \int
  638 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  639 + * \f]
  640 + *
618 641 * and
619 642 * \f$a_i\f$ is the model parameter \f$i\f$.
620   - * Furthermore,
621   - * \f$S(\vec{p}, E, t)\f$ is the source model,
622   - * \f$PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t)\f$ is the point
623   - * spread function,
624   - * \f$\vec{p'}\f$ is the measured photon direction,
  643 + * Furthermore
  644 + * \f$S(p,E,t)\f$ is the source model,
  645 + * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  646 + * \f$p'\f$ is the measured photon direction,
625 647 * \f$E'\f$ is the measured photon energy,
626 648 * \f$t'\f$ is the measured photon arrival time,
627   - * \f$\vec{p}\f$ is the true photon arrival direction,
628   - * \f$E\f$ is the true photon energy,
629   - * \f$t\f$ is the true photon arrival time, and
630   - * \f$d\f$ is the instrument pointing.
  649 + * \f$p\f$ is the true photon arrival direction,
  650 + * \f$E\f$ is the true photon energy, and
  651 + * \f$t\f$ is the true photon arrival time.
631 652 *
632 653 * This method uses a robust but simple difference method to estimate
633 654 * parameter gradients that have not been provided by the model. We use here
... ... @@ -832,7 +853,9 @@ void GObservation::free_members(void)
832 853 * optimisation using unbinned analysis and Poisson statistics.
833 854 * The -(log-likelihood) function is given by
834 855 *
835   - * \f$L = N_{\rm pred} - \sum_i \log e_i\f$
  856 + * \f[
  857 + * L = N_{\rm pred} - \sum_i \log e_i
  858 + * \f]
836 859 *
837 860 * where
838 861 * \f$N_{\rm pred}\f$ is the number of events predicted by the model, and
... ... @@ -949,7 +972,11 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models,
949 972 * This method evaluates the -(log-likelihood) function for parameter
950 973 * optimisation using binned analysis and Poisson statistics.
951 974 * The -(log-likelihood) function is given by
952   - * \f$L=-\sum_i n_i \log e_i - e_i\f$
  975 + *
  976 + * \f[
  977 + * L=-\sum_i n_i \log e_i - e_i
  978 + * \f]
  979 + *
953 980 * where the sum is taken over all data space bins, \f$n_i\f$ is the
954 981 * observed number of counts and \f$e_i\f$ is the model.
955 982 * This method also computes the parameter gradients
... ... @@ -1154,7 +1181,11 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
1154 1181 * This method evaluates the -(log-likelihood) function for parameter
1155 1182 * optimisation using binned analysis and Poisson statistics.
1156 1183 * The -(log-likelihood) function is given by
1157   - * \f$L = 1/2 \sum_i (n_i - e_i)^2 \sigma_i^{-2}\f$
  1184 + *
  1185 + * \f[
  1186 + * L = 1/2 \sum_i (n_i - e_i)^2 \sigma_i^{-2}
  1187 + * \f]
  1188 + *
1158 1189 * where the sum is taken over all data space bins, \f$n_i\f$ is the
1159 1190 * observed number of counts, \f$e_i\f$ is the model and \f$\sigma_i\f$
1160 1191 * is the statistical uncertainty.
... ... @@ -1281,8 +1312,6 @@ double GObservation::likelihood_gaussian_binned(const GModels& models,
1281 1312 = =
1282 1313 ==========================================================================*/
1283 1314  
1284   -
1285   -
1286 1315 /***********************************************************************//**
1287 1316 * @brief Model function evaluation for gradient computation
1288 1317 *
... ... @@ -1351,56 +1380,68 @@ double GObservation::npred_kern::eval(const double& x)
1351 1380 * Energy range is invalid.
1352 1381 *
1353 1382 * Computes
1354   - * \f[N'_{\rm pred} = \int_{E_{\rm bounds}} \int_{\rm ROI}
1355   - * S(\vec{p}, E, t) PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t) \,
1356   - * {\rm d}\vec{p'} {\rm d}E'\f]
1357   - * where
1358   - * \f$S(\vec{p}, E, t)\f$ is the source model,
1359   - * \f$PSF(\vec{p'}, E', t' | \vec{d}, \vec{p}, E, t)\f$ is the point
1360   - * spread function,
1361   - * \f$\vec{p'}\f$ is the measured photon direction,
  1383 + *
  1384 + * \f[
  1385 + * N_{\rm pred}(t') = \int_{E_{\rm bounds}} \int_{\rm ROI}
  1386 + * P(p',E',t') \, dp' \, dE'
  1387 + * \f]
  1388 + *
  1389 + * of the event probability
  1390 + *
  1391 + * \f[
  1392 + * P(p',E',t') = \int \int \int
  1393 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  1394 + * \f]
  1395 + *
  1396 + * where
  1397 + * \f$S(p,E,t)\f$ is the source model,
  1398 + * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  1399 + * \f$p'\f$ is the measured photon direction,
1362 1400 * \f$E'\f$ is the measured photon energy,
1363 1401 * \f$t'\f$ is the measured photon arrival time,
1364   - * \f$\vec{p}\f$ is the true photon arrival direction,
1365   - * \f$E\f$ is the true photon energy,
1366   - * \f$t\f$ is the true photon arrival time, and
1367   - * \f$d\f$ is the instrument pointing.
  1402 + * \f$p\f$ is the true photon arrival direction,
  1403 + * \f$E\f$ is the true photon energy, and
  1404 + * \f$t\f$ is the true photon arrival time.
1368 1405 *
1369   - * \f$E_{\rm bounds}\f$ are the energy boundaries that are stored in the
1370   - * GObservation::m_ebounds member.
1371   - *
1372   - * @todo Loop also over energy boundaries (is more general; there is no
1373   - * reason for not doing it).
  1406 + * This method performs the energy intergration over the energy boundaries
  1407 + * \f$E_{\rm bounds}\f$.
1374 1408 ***************************************************************************/
1375 1409 double GObservation::npred_spec(const GModel& model,
1376 1410 const GTime& obsTime) const
1377 1411 {
1378   - // Set integration energy interval in MeV
1379   - GEnergy e_min = events()->ebounds().emin();
1380   - GEnergy e_max = events()->ebounds().emax();
  1412 + // Initialise result
  1413 + double result = 0.0;
1381 1414  
1382   - // Get energy interval in MeV
1383   - double emin = e_min.MeV();
1384   - double emax = e_max.MeV();
  1415 + // Get energy boundaries
  1416 + GEbounds ebounds = events()->ebounds();
1385 1417  
1386   - // Throw exception if energy range is not valid
1387   - if (emax <= emin) {
1388   - throw GException::erange_invalid(G_NPRED_SPEC, emin, emax);
1389   - }
  1418 + // Loop over energy boundaries
  1419 + for (int i = 0; i < ebounds.size(); ++i) {
1390 1420  
1391   - // Setup integration function
1392   - GObservation::npred_spec_kern integrand(this, &model, &obsTime);
1393   - GIntegral integral(&integrand);
  1421 + // Get boundaries in MeV
  1422 + double emin = ebounds.emin(i).MeV();
  1423 + double emax = ebounds.emax(i).MeV();
1394 1424  
1395   - // Set integration precision
1396   - integral.eps(1.0e-6); // Needed for fluctuating bgd. model !!!
  1425 + // Continue only if valid
  1426 + if (emax > emin) {
1397 1427  
1398   - // Do Romberg integration
1399   - #if defined(G_LN_ENERGY_INT)
1400   - emin = std::log(emin);
1401   - emax = std::log(emax);
1402   - #endif
1403   - double result = integral.romberg(emin, emax);
  1428 + // Setup integration function
  1429 + GObservation::npred_spec_kern integrand(this, &model, &obsTime);
  1430 + GIntegral integral(&integrand);
  1431 +
  1432 + // Set integration precision
  1433 + integral.eps(1.0e-6); // Needed for fluctuating bgd. model !!!
  1434 +
  1435 + // Do Romberg integration
  1436 + #if defined(G_LN_ENERGY_INT)
  1437 + emin = std::log(emin);
  1438 + emax = std::log(emax);
  1439 + #endif
  1440 + result += integral.romberg(emin, emax);
  1441 +
  1442 + } // endif: energy interval was valid
  1443 +
  1444 + } // endfor: looped over energy boundaries
1404 1445  
1405 1446 // Compile option: Check for NaN/Inf
1406 1447 #if defined(G_NAN_CHECK)
... ... @@ -1408,8 +1449,6 @@ double GObservation::npred_spec(const GModel&amp; model,
1408 1449 std::cout << "*** ERROR: GObservation::npred_spec:";
1409 1450 std::cout << " NaN/Inf encountered";
1410 1451 std::cout << " (result=" << result;
1411   - std::cout << ", emin=" << emin;
1412   - std::cout << ", emax=" << emax;
1413 1452 std::cout << ")" << std::endl;
1414 1453 }
1415 1454 #endif
... ...