Leonardo Di Venere / gammalib

Commit 2a959989ed64e1b802f8e7dd096aeb23e9220c52

Authored by Jürgen Knödlseder
1 parent b94ea338

CTA background rates are now per ontime (#2758)

Inline Side-by-side
Showing 13 changed files with 146 additions and 214 deletions
ChangeLog
  View file @2a95998
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + CTA background rates are now per ontime (#2758)
6 7 Remove diffuse response cube caching (#2722)
7 8 Revise GCTAEdispCube::operator() interface and disable stacked response
8 9 cutting in true energy (#2721)
... ...
... ... @@ -315,6 +315,10 @@ None
315 315  
316 316 20. CTA interface
317 317 -----------------
  318 +CTA background rates are now per ontime and not livetime. This makes the IRF
  319 +background consistent with the cube background, and is also consistent with
  320 +the open gamma data format (#2758).
  321 +
318 322 The classes GCTACubeSourcePoint and GCTACubeSourceDiffuse were removed since
319 323 the response caching is now performed at the GCTAResponse level, hence no
320 324 dedicated response cache classes are needed anymore (#2722).
... ...
doc/source/admin/release_history/1.6.rst
  View file @2a95998
... ... @@ -123,6 +123,8 @@ New features
123 123 Clean up
124 124 --------
125 125  
  126 +* [`2758 <https://cta-redmine.irap.omp.eu/issues/2758>`_] -
  127 + CTA background rates are now per ontime
126 128 * [`2716 <https://cta-redmine.irap.omp.eu/issues/2716>`_] -
127 129 Rename ``GCTAResponseIrf::lo_save_thres()`` and
128 130 ``GCTAResponseIrf::hi_save_thres()`` to ``GCTAResponseIrf::lo_safe_thres()``
... ...
doc/source/users/user_manual/modules/cta_overview.rst
  View file @2a95998
... ... @@ -70,7 +70,7 @@ Function by :doxy:`GCTAResponseIrf` or as a precomputed response for
70 70 stacked analysis by :doxy:`GCTAResponseCube`.
71 71 For stacked analysis the response is composed of the exposure
72 72 (implemented by :doxy:`GCTACubeExposure`), the livetime averaged point
73   -spread function (implemented by :doxy:`GCTACubePsf`) and the livetime
  73 +spread function (implemented by :doxy:`GCTACubePsf`) and the ontime
74 74 averaged background rate (implemented by :doxy:`GCTACubeBackground`).
75 75 For all other analysis, the factorized IRFs are used directly.
76 76 Dependent on the response format, the components of the factorization
... ...
inst/cta/src/GCTABackground3D.cpp
  View file @2a95998
... ... @@ -172,15 +172,17 @@ GCTABackground3D&amp; GCTABackground3D::operator=(const GCTABackground3D&amp; bgd)
172 172  
173 173  
174 174 /***********************************************************************//**
175   - * @brief Return background rate in units of events/s/MeV/sr
  175 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  176 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
176 177 *
177 178 * @param[in] logE Log10 of the true photon energy (TeV).
178 179 * @param[in] detx Tangential X coordinate in nominal system (radians).
179 180 * @param[in] dety Tangential Y coordinate in nominal system (radians).
  181 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
180 182 *
181   - * Returns the background rate in units of events/s/MeV/sr for a given energy
182   - * and detector coordinates. The operator interpolates linearly in DETX and
183   - * DETY, and logarithmically in energy.
  183 + * Returns the background rate for a given energy and detector coordinates.
  184 + * The rate is given per ontime. The operator interpolates linearly in
  185 + * ``DETX`` and ``DETY``, and logarithmically in energy.
184 186 *
185 187 * The operator assures that the background rate never becomes negative. For
186 188 * invalid background models or detector coordinates outside the range
... ... @@ -285,8 +287,10 @@ GCTABackground3D* GCTABackground3D::clone(void) const
285 287 * ENERG_HI - Energy upper bin boundaries
286 288 * BKG - Background template (or BGD as legacy column name)
287 289 *
288   - * The data are stored in the m_background member. The DETX and DETY axes
289   - * will be set to radians, the energy axis will be set to log10.
  290 + * The data are stored in the m_background member. The ``DETX`` and ``DETY``
  291 + * axes will be set to radians, the energy axis will be set to log10. The
  292 + * data are expected to be given in unit of events per MeV, second and
  293 + * steradian, where the time is the real ontime.
290 294 *
291 295 * This method ignores all column names that are not the mandatory column
292 296 * names in the FITS @p table.
... ... @@ -310,7 +314,9 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
310 314 }
311 315 }
312 316  
313   - // Read background table from reduced FITS table
  317 + // Read background table from reduced FITS table. The background table
  318 + // is expected to provide the number of events per MeV, second and
  319 + // steradian, where the time is the real ontime.
314 320 m_background.read(*ptr);
315 321  
316 322 // Delete reduced FITS table
... ... @@ -539,15 +545,17 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy&amp; energy,
539 545  
540 546  
541 547 /***********************************************************************//**
542   - * @brief Returns background count rate integrated over energy interval
  548 + * @brief Returns background rate integrated over energy interval in units
  549 + * of events s\f$^{-1}\f$ sr\f$^{-1}\f$
543 550 *
544 551 * @param[in] dir Instrument direction.
545 552 * @param[in] emin Minimum energy of energy interval.
546 553 * @param[in] emax Maximum energy of energy interval.
547   - * @return Integrated background count rate (counts/src/sr).
  554 + * @return Integrated background count rate (events s\f$^{-1}\f$ sr\f$^{-1}\f$).
548 555 *
549 556 * Returns the background count rate for a given instrument direction that
550   - * is integrated over a specified energy interval.
  557 + * is integrated over a specified energy interval. The rate is given per
  558 + * ontime.
551 559 *
552 560 * If the energy interval is not positive, a zero background rate is
553 561 * returned.
... ... @@ -568,10 +576,6 @@ double GCTABackground3D::rate_ebin(const GCTAInstDir&amp; dir,
568 576 (dir.dety() <= m_dety_max) &&
569 577 (emax > emin)) {
570 578  
571   - // Retrieve references to node arrays
572   - //const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx);
573   - //const GNodeArray& detys = m_background.axis_nodes(m_inx_dety);
574   -
575 579 // Initialise first and second node
576 580 double x1 = emin.MeV();
577 581 double x2 = 0.0;
... ... @@ -1046,7 +1050,7 @@ void GCTABackground3D::init_mc_cache(void) const
1046 1050 #if defined(G_DEBUG_MC_INIT)
1047 1051 std::cout << "Energy=" << energies[i];
1048 1052 std::cout << " Rate_total=" << total_flux;
1049   - std::cout << " events/s/MeV" << std::endl;
  1053 + std::cout << " events/(s MeV)" << std::endl;
1050 1054 #endif
1051 1055  
1052 1056 } // endfor: looped over energy bins
... ...
inst/cta/src/GCTABackgroundPerfTable.cpp
  View file @2a95998
... ... @@ -161,15 +161,16 @@ GCTABackgroundPerfTable&amp; GCTABackgroundPerfTable::operator=(const GCTABackground
161 161  
162 162  
163 163 /***********************************************************************//**
164   - * @brief Return background rate in units of events/s/MeV/sr
  164 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  165 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
165 166 *
166 167 * @param[in] logE Log10 of the true photon energy (TeV).
167 168 * @param[in] detx Tangential coordinate in nominal system (rad).
168 169 * @param[in] dety Tangential coordinate in nominal system (rad).
169   - * @return Background rate (events/s/MeV/sr)
  170 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
170 171 *
171   - * Returns the background rate in units of events/s/MeV/sr for a given energy
172   - * and detector coordinates. The method assures that the background rate
  172 + * Returns the background rate for a given energy and detector coordinates.
  173 + * The rate is given per ontime. The method assures that the background rate
173 174 * never becomes negative.
174 175 *
175 176 * If the performance table contains less than 2 nodes an empty value is
... ... @@ -377,15 +378,17 @@ GCTAInstDir GCTABackgroundPerfTable::mc(const GEnergy&amp; energy,
377 378  
378 379  
379 380 /***********************************************************************//**
380   - * @brief Returns background count rate integrated over energy interval
  381 + * @brief Returns background rate integrated over energy interval in units
  382 + * of events s\f$^{-1}\f$ sr\f$^{-1}\f$
381 383 *
382 384 * @param[in] dir Instrument direction.
383 385 * @param[in] emin Minimum energy of energy interval.
384 386 * @param[in] emax Maximum energy of energy interval.
385   - * @return Integrated background count rate (counts/src/sr).
  387 + * @return Integrated background count rate (events s\f$^{-1}\f$ sr\f$^{-1}\f$).
386 388 *
387 389 * Returns the background count rate for a given instrument direction that
388   - * is integrated over a specified energy interval.
  390 + * is integrated over a specified energy interval. The rate is given per
  391 + * ontime.
389 392 *
390 393 * If the energy interval is not positive, a zero background rate is
391 394 * returned.
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @2a95998
... ... @@ -240,13 +240,15 @@ GCTACubeBackground&amp; GCTACubeBackground::operator=(const GCTACubeBackground&amp; bgd)
240 240  
241 241  
242 242 /***********************************************************************//**
243   - * @brief Return background rate in units of events/s/MeV/sr
  243 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  244 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
244 245 *
245 246 * @param[in] dir Reconstructed event position where rate should be computed
246 247 * @param[in] energy Energy at which the rate should be computed
  248 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
247 249 *
248   - * Returns the background rate in units of events/s/MeV/sr for a given energy
249   - * and sky coordinate. The method assures that the background rate
  250 + * Returns the background rate for a giben instrument direction and energy.
  251 + * The rate is given per ontime. The method assures that the background rate
250 252 * never becomes negative.
251 253 *
252 254 * The method interpolates logarithmically in the energy direction.
... ... @@ -334,7 +336,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
334 336 * @exception GException::invalid_value
335 337 * No event list found in CTA observations.
336 338 *
337   - * Set the background cube by computing the livetime weighted background rate
  339 + * Set the background cube by computing the ontime weighted background rate
338 340 * for all CTA observations in an observation container. The cube pixel
339 341 * values are computed as the sum over the background rates.
340 342 ***************************************************************************/
... ... @@ -354,8 +356,8 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
354 356 // Initialise event cubes for model evaluate and storage
355 357 GCTAEventCube eventcube = GCTAEventCube(m_cube, m_ebounds, gti);
356 358  
357   - // Initialise total livetime
358   - double total_livetime = 0.0;
  359 + // Initialise total ontime
  360 + double total_ontime = 0.0;
359 361  
360 362 // Loop over all observations in container
361 363 for (int i = 0; i < obs.size(); ++i) {
... ... @@ -388,17 +390,17 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
388 390 continue;
389 391 }
390 392  
391   - // Get observation livetime
392   - double livetime = cta->livetime();
  393 + // Get observation ontime
  394 + double ontime = cta->ontime();
393 395  
394   - // Skip observation if livetime is zero
395   - if (livetime == 0.0) {
  396 + // Skip observation if ontime is zero
  397 + if (ontime == 0.0) {
396 398 if (log != NULL) {
397 399 *log << "Skipping unbinned ";
398 400 *log << cta->instrument();
399 401 *log << " observation ";
400 402 *log << "\"" << cta->name() << "\"";
401   - *log << " (id=" << cta->id() << ") due to zero livetime";
  403 + *log << " (id=" << cta->id() << ") due to zero ontime";
402 404 *log << std::endl;
403 405 }
404 406 continue;
... ... @@ -499,15 +501,15 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
499 501 double x2 = m_ebounds.emax(iebin).MeV();
500 502 double model = gammalib::plaw_integral(x1, f1, x2, f2) / (x2 - x1);
501 503  
502   - // Multiply by livetime to get the correct weighting for
503   - // each observation. We divide by the total livetime later
504   - // to get the background model in units of counts/MeV/s/sr.
505   - model *= livetime;
  504 + // Multiply by ontime to get the correct weighting for each
  505 + // observation. We divide by the total ontime later to get
  506 + // the background model in units of counts/(MeV s sr).
  507 + model *= ontime;
506 508  
507 509 // Add existing number of counts
508 510 model += bin->counts();
509 511  
510   - // Store cumulated value (units: counts/MeV/sr)
  512 + // Store cumulated value (units: counts/(MeV sr))
511 513 bin->counts(model);
512 514  
513 515 } // endif: added background model value to cube
... ... @@ -519,26 +521,26 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
519 521  
520 522 } // endfor: looped over spatial pixels
521 523  
522   - // Accumulate livetime
523   - total_livetime += livetime;
  524 + // Accumulate ontime
  525 + total_ontime += ontime;
524 526  
525 527 } // endfor: looped over all observations
526 528  
527   - // Re-normalize cube to get units of counts/MeV/s/sr
528   - if (total_livetime > 0.0) {
  529 + // Re-normalize cube to get units of counts/(MeV s sr)
  530 + if (total_ontime > 0.0) {
529 531  
530 532 // Loop over all bins in background cube and divide the content
531   - // by the total livetime.
  533 + // by the total ontime.
532 534 for (int i = 0; i < eventcube.size(); ++i) {
533 535 GCTAEventBin* bin = eventcube[i];
534   - double rate = bin->counts() / total_livetime;
  536 + double rate = bin->counts() / total_ontime;
535 537 bin->counts(rate);
536 538 }
537 539  
538 540 // Set background cube values from event cube
539 541 m_cube = eventcube.counts();
540 542  
541   - } // endif: livetime was positive
  543 + } // endif: ontime was positive
542 544  
543 545 // Return
544 546 return;
... ... @@ -547,10 +549,12 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
547 549  
548 550  
549 551 /***********************************************************************//**
550   - * @brief Compute spatially integrated background rate for a given energy
  552 + * @brief Return spatially integrated background rate in units of
  553 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
551 554 *
552 555 * @param[in] logE Logarithm (base 10) of energy in TeV
553   - * @return Spatially integrated background rate in units of counts/MeV/s
  556 + * @return Spatially integrated background rate
  557 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
554 558 *
555 559 * Spatially integrates the background cube at a given energy. This method
556 560 * performs an interpolation between the energy maps. The integration is
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @2a95998
... ... @@ -351,17 +351,19 @@ void GCTAModelAeffBackground::temporal(const GModelTemporal* temporal)
351 351  
352 352  
353 353 /***********************************************************************//**
354   - * @brief Evaluate function
  354 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  355 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
355 356 *
356 357 * @param[in] event Observed event.
357 358 * @param[in] obs Observation.
358 359 * @param[in] gradients Compute gradients?
359   - * @return Function value.
  360 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
  361 + *
  362 + * The method returns a real rate, defined by the number of counts per MeV,
  363 + * steradian and ontime.
360 364 *
361 365 * If the @p gradients flag is true the method will also set the parameter
362 366 * gradients of the model parameters.
363   - *
364   - * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
365 367 ***************************************************************************/
366 368 double GCTAModelAeffBackground::eval(const GEvent& event,
367 369 const GObservation& obs,
... ... @@ -392,16 +394,12 @@ double GCTAModelAeffBackground::eval(const GEvent&amp; event,
392 394 // Compute value
393 395 double value = spat * spec * temp;
394 396  
395   - // Apply deadtime correction
396   - double deadc = obs.deadc(event.time());
397   - value *= deadc;
398   -
399 397 // Optionally compute partial derivatives
400 398 if (gradients) {
401 399  
402 400 // Multiply factors to spectral gradients
403 401 if (spectral() != NULL) {
404   - double fact = spat * temp * deadc;
  402 + double fact = spat * temp;
405 403 if (fact != 1.0) {
406 404 for (int i = 0; i < spectral()->size(); ++i)
407 405 (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
... ... @@ -410,7 +408,7 @@ double GCTAModelAeffBackground::eval(const GEvent&amp; event,
410 408  
411 409 // Multiply factors to temporal gradients
412 410 if (temporal() != NULL) {
413   - double fact = spat * spec * deadc;
  411 + double fact = spat * spec;
414 412 if (fact != 1.0) {
415 413 for (int i = 0; i < temporal()->size(); ++i)
416 414 (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
... ... @@ -425,20 +423,21 @@ double GCTAModelAeffBackground::eval(const GEvent&amp; event,
425 423  
426 424  
427 425 /***********************************************************************//**
428   - * @brief Return spatially integrated background model
  426 + * @brief Return spatially integrated background rate in units of
  427 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
429 428 *
430 429 * @param[in] obsEng Measured event energy.
431 430 * @param[in] obsTime Measured event time.
432 431 * @param[in] obs Observation.
433   - * @return Spatially integrated model.
  432 + * @return Spatially integrated background rate
  433 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
434 434 *
435 435 * @exception GException::invalid_value
436 436 * Pointing direction differs from RoI centre.
437 437 *
438 438 * Spatially integrates the effective area background model for a given
439   - * measured event energy and event time. This method also applies a deadtime
440   - * correction factor, so that the normalization of the model is a real rate
441   - * (counts/MeV/s).
  439 + * measured event energy and event time. The method returns a real rate,
  440 + * defined as the number of counts per MeV and ontime.
442 441 ***************************************************************************/
443 442 double GCTAModelAeffBackground::npred(const GEnergy& obsEng,
444 443 const GTime& obsTime,
... ... @@ -532,9 +531,6 @@ double GCTAModelAeffBackground::npred(const GEnergy&amp; obsEng,
532 531 npred *= spectral()->eval(obsEng, obsTime);
533 532 npred *= temporal()->eval(obsTime);
534 533  
535   - // Apply deadtime correction
536   - npred *= obs.deadc(obsTime);
537   -
538 534 // Return Npred
539 535 return npred;
540 536 }
... ... @@ -555,10 +551,6 @@ double GCTAModelAeffBackground::npred(const GEnergy&amp; obsEng,
555 551 * number generator of type GRan which is passed by reference, hence the
556 552 * state of the random number generator will be changed by the method.
557 553 *
558   - * The method also applies a deadtime correction using a Monte Carlo process,
559   - * taking into account temporal deadtime variations. For this purpose, the
560   - * method makes use of the time dependent GObservation::deadc method.
561   - *
562 554 * For each event in the returned event list, the sky direction, the nominal
563 555 * coordinates (DETX and DETY), the energy and the time will be set.
564 556 ***************************************************************************/
... ... @@ -650,24 +642,12 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
650 642 #if defined(G_DUMP_MC)
651 643 std::cout << " Interval " << itime;
652 644 std::cout << " times=" << n_events << std::endl;
653   - int n_killed_by_deadtime = 0;
654 645 int n_killed_by_roi = 0;
655 646 #endif
656 647  
657 648 // Loop over events
658 649 for (int i = 0; i < n_events; ++i) {
659 650  
660   - // Apply deadtime correction
661   - double deadc = obs.deadc(times[i]);
662   - if (deadc < 1.0) {
663   - if (ran.uniform() > deadc) {
664   - #if defined(G_DUMP_MC)
665   - n_killed_by_deadtime++;
666   - #endif
667   - continue;
668   - }
669   - }
670   -
671 651 // Get Monte Carlo event energy from spectral model
672 652 GEnergy energy = spectral.mc(ebounds.emin(ieng),
673 653 ebounds.emax(ieng),
... ... @@ -759,8 +739,6 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
759 739  
760 740 // Debug option: provide statisics
761 741 #if defined(G_DUMP_MC)
762   - std::cout << " Killed by deadtime=";
763   - std::cout << n_killed_by_deadtime << std::endl;
764 742 std::cout << " Killed by ROI=";
765 743 std::cout << n_killed_by_roi << std::endl;
766 744 #endif
... ...
inst/cta/src/GCTAModelBackground.cpp
  View file @2a95998
... ... @@ -379,17 +379,17 @@ void GCTAModelBackground::temporal(const GModelTemporal* temporal)
379 379  
380 380  
381 381 /***********************************************************************//**
382   - * @brief Evaluate function
  382 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  383 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
383 384 *
384 385 * @param[in] event Observed event.
385 386 * @param[in] obs Observation.
386 387 * @param[in] gradients Compute gradients?
387   - * @return Function value.
  388 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
388 389 *
389 390 * Evaluates the background model which is a factorization of a spatial,
390   - * spectral and temporal model component. This method also applies a deadtime
391   - * correction factor, so that the normalization of the model is a real rate
392   - * (counts/exposure time).
  391 + * spectral and temporal model component. The method returns a real rate,
  392 + * defined as the number of counts per MeV, steradian and ontime.
393 393 *
394 394 * If the @p gradients flag is true the method will also set the parameter
395 395 * gradients of the model parameters.
... ... @@ -422,16 +422,12 @@ double GCTAModelBackground::eval(const GEvent&amp; event,
422 422 // Compute value
423 423 double value = spat * spec * temp;
424 424  
425   - // Apply deadtime correction
426   - double deadc = obs.deadc(event.time());
427   - value *= deadc;
428   -
429 425 // Optionally compute partial derivatives
430 426 if (gradients) {
431 427  
432 428 // Multiply factors to radial gradients
433 429 if (spatial() != NULL) {
434   - double fact = spec * temp * deadc;
  430 + double fact = spec * temp;
435 431 if (fact != 1.0) {
436 432 for (int i = 0; i < spatial()->size(); ++i)
437 433 (*spatial())[i].factor_gradient((*spatial())[i].factor_gradient() * fact);
... ... @@ -440,7 +436,7 @@ double GCTAModelBackground::eval(const GEvent&amp; event,
440 436  
441 437 // Multiply factors to spectral gradients
442 438 if (spectral() != NULL) {
443   - double fact = spat * temp * deadc;
  439 + double fact = spat * temp;
444 440 if (fact != 1.0) {
445 441 for (int i = 0; i < spectral()->size(); ++i)
446 442 (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
... ... @@ -449,7 +445,7 @@ double GCTAModelBackground::eval(const GEvent&amp; event,
449 445  
450 446 // Multiply factors to temporal gradients
451 447 if (temporal() != NULL) {
452   - double fact = spat * spec * deadc;
  448 + double fact = spat * spec;
453 449 if (fact != 1.0) {
454 450 for (int i = 0; i < temporal()->size(); ++i)
455 451 (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
... ... @@ -464,20 +460,21 @@ double GCTAModelBackground::eval(const GEvent&amp; event,
464 460  
465 461  
466 462 /***********************************************************************//**
467   - * @brief Return spatially integrated background model
  463 + * @brief Return spatially integrated background rate in units of
  464 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
468 465 *
469 466 * @param[in] energy Measured event energy.
470 467 * @param[in] time Measured event time.
471 468 * @param[in] obs Observation.
472   - * @return Integral of spatial model component.
  469 + * @return Spatially integrated background rate
  470 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
473 471 *
474 472 * @exception GException::invalid_value
475 473 * Pointing direction differs from RoI centre.
476 474 *
477 475 * Spatially integrates the background model for a given measured event
478   - * energy and event time. This method also applies a deadtime correction
479   - * factor, so that the normalization of the model is a real rate
480   - * (counts/exposure time).
  476 + * energy and event time. The method returns a real rate, defined as the
  477 + * number of counts per MeV and ontime.
481 478 ***************************************************************************/
482 479 double GCTAModelBackground::npred(const GEnergy& energy,
483 480 const GTime& time,
... ... @@ -511,9 +508,6 @@ double GCTAModelBackground::npred(const GEnergy&amp; energy,
511 508 npred *= spectral()->eval(energy, time);
512 509 npred *= temporal()->eval(time);
513 510  
514   - // Apply deadtime correction
515   - npred *= obs.deadc(time);
516   -
517 511 // Return Npred
518 512 return npred;
519 513 }
... ... @@ -534,10 +528,6 @@ double GCTAModelBackground::npred(const GEnergy&amp; energy,
534 528 * number generator of type GRan which is passed by reference, hence the
535 529 * state of the random number generator will be changed by the method.
536 530 *
537   - * The method also applies a deadtime correction using a Monte Carlo process,
538   - * taking into account temporal deadtime variations. For this purpose, the
539   - * method makes use of the time dependent GObservation::deadc method.
540   - *
541 531 * For each event in the returned event list, the sky direction, the nominal
542 532 * coordinates (DETX and DETY), the energy and the time will be set.
543 533 *
... ... @@ -624,24 +614,12 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
624 614 #if defined(G_DUMP_MC)
625 615 std::cout << " Interval " << itime;
626 616 std::cout << " events=" << n_events << std::endl;
627   - int n_removed_by_deadtime = 0;
628   - int n_trial_outside_roi = 0;
  617 + int n_trial_outside_roi = 0;
629 618 #endif
630 619  
631 620 // Loop over events
632 621 for (int i = 0; i < n_events; ++i) {
633 622  
634   - // Apply deadtime correction
635   - double deadc = obs.deadc(times[i]);
636   - if (deadc < 1.0) {
637   - if (ran.uniform() > deadc) {
638   - #if defined(G_DUMP_MC)
639   - n_removed_by_deadtime++;
640   - #endif
641   - continue;
642   - }
643   - }
644   -
645 623 // Get Monte Carlo event energy from spectral model
646 624 GEnergy energy = spectral->mc(ebounds.emin(ieng),
647 625 ebounds.emax(ieng),
... ... @@ -684,8 +662,6 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
684 662  
685 663 // Debug option: provide statisics
686 664 #if defined(G_DUMP_MC)
687   - std::cout << " Events removed due to deadtime=";
688   - std::cout << n_removed_by_deadtime << std::endl;
689 665 std::cout << " Trials outside RoI=";
690 666 std::cout << n_trial_outside_roi << std::endl;
691 667 #endif
... ...
inst/cta/src/GCTAModelCubeBackground.cpp
  View file @2a95998
... ... @@ -326,12 +326,16 @@ void GCTAModelCubeBackground::temporal(const GModelTemporal* temporal)
326 326  
327 327  
328 328 /***********************************************************************//**
329   - * @brief Evaluate function
  329 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  330 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
330 331 *
331 332 * @param[in] event Observed event.
332 333 * @param[in] obs Observation.
333 334 * @param[in] gradients Compute gradients?
334   - * @return Function value.
  335 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
  336 + *
  337 + * Evaluates the background model. The method returns a real rate, defined
  338 + * as the number of counts per MeV, steradian and ontime.
335 339 *
336 340 * If the @p gradients flag is true the method will also set the parameter
337 341 * gradients of the model parameters.
... ... @@ -369,7 +373,7 @@ double GCTAModelCubeBackground::eval(const GEvent&amp; event,
369 373 double fact = spat * temp;
370 374 if (fact != 1.0) {
371 375 for (int i = 0; i < spectral()->size(); ++i)
372   - (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  376 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
373 377 }
374 378 }
375 379  
... ... @@ -378,7 +382,7 @@ double GCTAModelCubeBackground::eval(const GEvent&amp; event,
378 382 double fact = spat * spec;
379 383 if (fact != 1.0) {
380 384 for (int i = 0; i < temporal()->size(); ++i)
381   - (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  385 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
382 386 }
383 387 }
384 388  
... ... @@ -390,17 +394,18 @@ double GCTAModelCubeBackground::eval(const GEvent&amp; event,
390 394  
391 395  
392 396 /***********************************************************************//**
393   - * @brief Return spatially integrated background model
  397 + * @brief Return spatially integrated background rate in units of
  398 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
394 399 *
395 400 * @param[in] obsEng Measured event energy.
396 401 * @param[in] obsTime Measured event time.
397 402 * @param[in] obs Observation.
398   - * @return Spatially integrated model.
  403 + * @return Spatially integrated background rate
  404 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
399 405 *
400 406 * Spatially integrates the cube background model for a given measured event
401   - * energy and event time. This method also applies a deadtime correction
402   - * factor, so that the normalization of the model is a real rate
403   - * (counts/MeV/s).
  407 + * energy and event time. The method returns a real rate,
  408 + * defined as the number of counts per MeV and ontime.
404 409 ***************************************************************************/
405 410 double GCTAModelCubeBackground::npred(const GEnergy& obsEng,
406 411 const GTime& obsTime,
... ... @@ -479,9 +484,6 @@ double GCTAModelCubeBackground::npred(const GEnergy&amp; obsEng,
479 484 npred *= spectral()->eval(obsEng, obsTime);
480 485 npred *= temporal()->eval(obsTime);
481 486  
482   - // Apply deadtime correction
483   - npred *= obs.deadc(obsTime);
484   -
485 487 // Return Npred
486 488 return npred;
487 489 }
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @2a95998
... ... @@ -323,17 +323,19 @@ void GCTAModelIrfBackground::temporal(const GModelTemporal* temporal)
323 323  
324 324  
325 325 /***********************************************************************//**
326   - * @brief Evaluate function
  326 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  327 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
327 328 *
328 329 * @param[in] event Observed event.
329 330 * @param[in] obs Observation.
330 331 * @param[in] gradients Compute gradients?
331   - * @return Function value.
  332 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
  333 + *
  334 + * The method returns a real rate, defined by the number of counts per MeV,
  335 + * steradian and ontime.
332 336 *
333 337 * If the @p gradients flag is true the method will also set the parameter
334 338 * gradients of the model parameters.
335   - *
336   - * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
337 339 ***************************************************************************/
338 340 double GCTAModelIrfBackground::eval(const GEvent& event,
339 341 const GObservation& obs,
... ... @@ -360,33 +362,25 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
360 362 // Compute value
361 363 double value = spat * spec * temp;
362 364  
363   - // Apply deadtime correction. This assumes that the IRF background rates
364   - // are stored per livetime, and by multiplying these rates with the
365   - // deadtime correction factor thay are given per ontime.
366   - /*
367   - double deadc = obs.deadc(event.time());
368   - value *= deadc;
369   - */
370   -
371 365 // If gradients were requested then multiply factors to spectral and
372 366 // temporal gradients
373 367 if (gradients) {
374 368  
375 369 // Multiply factors to spectral gradients
376 370 if (spectral() != NULL) {
377   - double fact = spat * temp; // * deadc;
  371 + double fact = spat * temp;
378 372 if (fact != 1.0) {
379 373 for (int i = 0; i < spectral()->size(); ++i)
380   - (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  374 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
381 375 }
382 376 }
383 377  
384 378 // Multiply factors to temporal gradients
385 379 if (temporal() != NULL) {
386   - double fact = spat * spec; // * deadc;
  380 + double fact = spat * spec;
387 381 if (fact != 1.0) {
388 382 for (int i = 0; i < temporal()->size(); ++i)
389   - (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  383 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
390 384 }
391 385 }
392 386  
... ... @@ -398,20 +392,21 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
398 392  
399 393  
400 394 /***********************************************************************//**
401   - * @brief Return spatially integrated background model
  395 + * @brief Return spatially integrated background rate in units of
  396 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
402 397 *
403 398 * @param[in] obsEng Measured event energy.
404 399 * @param[in] obsTime Measured event time.
405 400 * @param[in] obs Observation.
406   - * @return Spatially integrated model.
  401 + * @return Spatially integrated background rate
  402 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
407 403 *
408 404 * @exception GException::invalid_value
409 405 * Pointing direction differs from RoI centre.
410 406 *
411 407 * Spatially integrates the instrumental background model for a given
412   - * measured event energy and event time. This method also applies a deadtime
413   - * correction factor, so that the normalization of the model is a real rate
414   - * (counts/MeV/s).
  408 + * measured event energy and event time. The method returns a real rate,
  409 + * defined as the number of counts per MeV and ontime.
415 410 ***************************************************************************/
416 411 double GCTAModelIrfBackground::npred(const GEnergy& obsEng,
417 412 const GTime& obsTime,
... ... @@ -531,9 +526,6 @@ double GCTAModelIrfBackground::npred(const GEnergy&amp; obsEng,
531 526 npred *= spectral()->eval(obsEng, obsTime);
532 527 npred *= temporal()->eval(obsTime);
533 528  
534   - // Apply deadtime correction
535   - npred *= obs.deadc(obsTime);
536   -
537 529 // Return Npred
538 530 return npred;
539 531 }
... ... @@ -554,10 +546,6 @@ double GCTAModelIrfBackground::npred(const GEnergy&amp; obsEng,
554 546 * number generator of type GRan which is passed by reference, hence the
555 547 * state of the random number generator will be changed by the method.
556 548 *
557   - * The method also applies a deadtime correction using a Monte Carlo process,
558   - * taking into account temporal deadtime variations. For this purpose, the
559   - * method makes use of the time dependent GObservation::deadc method.
560   - *
561 549 * For each event in the returned event list, the sky direction, the nominal
562 550 * coordinates (DETX and DETY), the energy and the time will be set.
563 551 ***************************************************************************/
... ... @@ -637,24 +625,12 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
637 625 #if defined(G_DUMP_MC)
638 626 std::cout << " Interval " << itime;
639 627 std::cout << " events=" << n_events << std::endl;
640   - int n_killed_by_deadtime = 0;
641   - int n_killed_by_roi = 0;
  628 + int n_killed_by_roi = 0;
642 629 #endif
643 630  
644 631 // Loop over events
645 632 for (int i = 0; i < n_events; ++i) {
646 633  
647   - // Apply deadtime correction
648   - double deadc = obs.deadc(times[i]);
649   - if (deadc < 1.0) {
650   - if (ran.uniform() > deadc) {
651   - #if defined(G_DUMP_MC)
652   - n_killed_by_deadtime++;
653   - #endif
654   - continue;
655   - }
656   - }
657   -
658 634 // Get Monte Carlo event energy from spectral model
659 635 GEnergy energy = spectral.mc(ebounds.emin(ieng),
660 636 ebounds.emax(ieng),
... ... @@ -693,8 +669,6 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
693 669  
694 670 // Debug option: provide statisics
695 671 #if defined(G_DUMP_MC)
696   - std::cout << " Killed by deadtime=";
697   - std::cout << n_killed_by_deadtime << std::endl;
698 672 std::cout << " Killed by ROI=";
699 673 std::cout << n_killed_by_roi << std::endl;
700 674 #endif
... ...
inst/cta/src/GCTAModelRadialAcceptance.cpp
  View file @2a95998
... ... @@ -384,17 +384,17 @@ void GCTAModelRadialAcceptance::temporal(const GModelTemporal* temporal)
384 384  
385 385  
386 386 /***********************************************************************//**
387   - * @brief Evaluate function
  387 + * @brief Return background rate in units of events MeV\f$^{-1}\f$
  388 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
388 389 *
389 390 * @param[in] event Observed event.
390 391 * @param[in] obs Observation.
391 392 * @param[in] gradients Compute gradients?
392   - * @return Function value.
  393 + * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
393 394 *
394 395 * Evaluates tha CTA radial acceptance model which is a factorization of a
395   - * spatial, spectral and temporal model component. This method also applies
396   - * a deadtime correction factor, so that the normalization of the model is
397   - * a real rate (counts/exposure time).
  396 + * spatial, spectral and temporal model component. The method returns a
  397 + * real rate, defined by the number of counts per MeV, steradian and ontime.
398 398 *
399 399 * If the @p gradients flag is true the method will also set the parameter
400 400 * gradients of the model parameters.
... ... @@ -432,37 +432,33 @@ double GCTAModelRadialAcceptance::eval(const GEvent&amp; event,
432 432 // Compute value
433 433 double value = rad * spec * temp;
434 434  
435   - // Apply deadtime correction
436   - double deadc = obs.deadc(event.time());
437   - value *= deadc;
438   -
439 435 // Optionally compute partial derivatives
440 436 if (gradients) {
441 437  
442 438 // Multiply factors to radial gradients
443 439 if (radial() != NULL) {
444   - double fact = spec * temp * deadc;
  440 + double fact = spec * temp;
445 441 if (fact != 1.0) {
446 442 for (int i = 0; i < radial()->size(); ++i)
447   - (*radial())[i].factor_gradient((*radial())[i].factor_gradient() * fact );
  443 + (*radial())[i].factor_gradient((*radial())[i].factor_gradient() * fact);
448 444 }
449 445 }
450 446  
451 447 // Multiply factors to spectral gradients
452 448 if (spectral() != NULL) {
453   - double fact = rad * temp * deadc;
  449 + double fact = rad * temp;
454 450 if (fact != 1.0) {
455 451 for (int i = 0; i < spectral()->size(); ++i)
456   - (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  452 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
457 453 }
458 454 }
459 455  
460 456 // Multiply factors to temporal gradients
461 457 if (temporal() != NULL) {
462   - double fact = rad * spec * deadc;
  458 + double fact = rad * spec;
463 459 if (fact != 1.0) {
464 460 for (int i = 0; i < temporal()->size(); ++i)
465   - (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  461 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
466 462 }
467 463 }
468 464  
... ... @@ -474,11 +470,14 @@ double GCTAModelRadialAcceptance::eval(const GEvent&amp; event,
474 470  
475 471  
476 472 /***********************************************************************//**
477   - * @brief Return spatially integrated data model
  473 + * @brief Return spatially integrated background rate in units of
  474 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
478 475 *
479 476 * @param[in] obsEng Measured event energy.
480 477 * @param[in] obsTime Measured event time.
481 478 * @param[in] obs Observation.
  479 + * @return Spatially integrated background rate
  480 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
482 481 *
483 482 * @exception GException::no_list
484 483 * No valid CTA event list found in observation
... ... @@ -486,8 +485,8 @@ double GCTAModelRadialAcceptance::eval(const GEvent&amp; event,
486 485 * No valid CTA pointing found in observation
487 486 *
488 487 * Spatially integrates the data model for a given measured event energy and
489   - * event time. This method also applies a deadtime correction factor, so that
490   - * the normalization of the model is a real rate (counts/exposure time).
  488 + * event time. The method returns a real rate, defined as the number of
  489 + * counts per MeV and ontime.
491 490 ***************************************************************************/
492 491 double GCTAModelRadialAcceptance::npred(const GEnergy& obsEng,
493 492 const GTime& obsTime,
... ... @@ -526,9 +525,6 @@ double GCTAModelRadialAcceptance::npred(const GEnergy&amp; obsEng,
526 525 npred *= spectral()->eval(obsEng, obsTime);
527 526 npred *= temporal()->eval(obsTime);
528 527  
529   - // Apply deadtime correction
530   - npred *= obs.deadc(obsTime);
531   -
532 528 } // endif: model was valid
533 529  
534 530 // Return
... ... @@ -551,10 +547,6 @@ double GCTAModelRadialAcceptance::npred(const GEnergy&amp; obsEng,
551 547 * argument that is passed to the method. The method also requires a random
552 548 * number generator of type GRan which is passed by reference, hence the
553 549 * state of the random number generator will be changed by the method.
554   - *
555   - * The method also applies a deadtime correction using a Monte Carlo process,
556   - * taking into account temporal deadtime variations. For this purpose, the
557   - * method makes use of the time dependent GObservation::deadc method.
558 550 ***************************************************************************/
559 551 GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation& obs,
560 552 GRan& ran) const
... ... @@ -614,14 +606,6 @@ GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation&amp; obs,
614 606 // Loop over events
615 607 for (int i = 0; i < n_events; ++i) {
616 608  
617   - // Apply deadtime correction
618   - double deadc = obs.deadc(times[i]);
619   - if (deadc < 1.0) {
620   - if (ran.uniform() > deadc) {
621   - continue;
622   - }
623   - }
624   -
625 609 // Set event energy
626 610 GEnergy energy = spectral()->mc(obs.events()->ebounds().emin(ieng),
627 611 obs.events()->ebounds().emax(ieng),
... ... @@ -781,7 +765,7 @@ void GCTAModelRadialAcceptance::write(GXmlElement&amp; xml) const
781 765 /***********************************************************************//**
782 766 * @brief Print model information
783 767 *
784   - * @param[in] chatter Chattiness (defaults to NORMAL).
  768 + * @param[in] chatter Chattiness.
785 769 * @return String containing model information.
786 770 ***************************************************************************/
787 771 std::string GCTAModelRadialAcceptance::print(const GChatter& chatter) const
... ...
inst/cta/test/test_CTA.cpp
  View file @2a95998
... ... @@ -2640,7 +2640,7 @@ void TestGCTAOptimize::test_unbinned_optimizer(void)
2640 2640 300000, 0,
2641 2641 1, 0,
2642 2642 2.94185, 0.0379307,
2643   - 6.43398e-05, 1.80715e-06,
  2643 + 6.11228e-05, 1.71680e-06,
2644 2644 -1.82866, 0.0155452,
2645 2645 1000000, 0,
2646 2646 1, 0};
... ... @@ -2684,7 +2684,7 @@ void TestGCTAOptimize::test_binned_optimizer(void)
2684 2684 300000, 0,
2685 2685 1, 0,
2686 2686 2.95662, 0.0704027,
2687   - 6.42985e-05, 1.96478e-06,
  2687 + 6.10835e-05, 1.86654e-06,
2688 2688 -1.82084, 0.0163749,
2689 2689 1000000, 0,
2690 2690 1, 0};
... ...