Commit aa56ca42d6065efb8100c911ffe42b201d23c984
Merge branch 'master' into devel
# Conflicts: # ChangeLog # NEWS # README.md # configure.ac # doc/Doxyfile # doc/source/admin/release_history/index.rst # doc/source/conf.py # sonar-project.properties
Showing
12 changed files
with
178 additions
and
41 deletions
ChangeLog
| ... | ... | @@ -51,6 +51,16 @@ |
| 51 | 51 | Implement pre-computation cache in GModelSpatialGauss (#1300) |
| 52 | 52 | |
| 53 | 53 | |
| 54 | +2021-04-28 | |
| 55 | + | |
| 56 | + * Bug fix version 1.7.4 released | |
| 57 | + ================================ | |
| 58 | + | |
| 59 | + Add conda activation and deactivation scripts for csh (#3628) | |
| 60 | + Speed-up Monte-Carlo sampling of exponentially cut-off power | |
| 61 | + law models (#3625) | |
| 62 | + | |
| 63 | + | |
| 54 | 64 | 2020-12-02 |
| 55 | 65 | |
| 56 | 66 | * Bug fix version 1.7.3 released | ... | ... |
Makefile.am
| 1 | 1 | ############################################################################# |
| 2 | 2 | # Makefile.am for GammaLib package # |
| 3 | 3 | # ------------------------------------------------------------------------- # |
| 4 | -# Copyright (C) 2016-2020 Juergen Knoedlseder # | |
| 4 | +# Copyright (C) 2016-2021 Juergen Knoedlseder # | |
| 5 | 5 | # ------------------------------------------------------------------------- # |
| 6 | 6 | # # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify # |
| ... | ... | @@ -38,9 +38,11 @@ pkgconfig_DATA = gammalib.pc |
| 38 | 38 | |
| 39 | 39 | # Conda activate directory |
| 40 | 40 | condaactivatedir = $(sysconfdir)/conda/activate.d |
| 41 | -condaactivate_DATA = dev/00-gammalib-activate.sh | |
| 41 | +condaactivate_DATA = dev/00-gammalib-activate.sh \ | |
| 42 | + dev/00-gammalib-activate.csh | |
| 42 | 43 | condadeactivatedir = $(sysconfdir)/conda/deactivate.d |
| 43 | -condadeactivate_DATA = dev/01-gammalib-deactivate.sh | |
| 44 | +condadeactivate_DATA = dev/01-gammalib-deactivate.sh \ | |
| 45 | + dev/01-gammalib-deactivate.csh | |
| 44 | 46 | |
| 45 | 47 | # Manual pages to install |
| 46 | 48 | man_MANS = |
| ... | ... | @@ -60,6 +62,7 @@ EXTRA_DIST = README.md \ |
| 60 | 62 | doc/Doxyfile doc/source \ |
| 61 | 63 | dev/testreport.py \ |
| 62 | 64 | dev/00-gammalib-activate.sh dev/01-gammalib-deactivate.sh \ |
| 65 | + dev/00-gammalib-activate.csh dev/01-gammalib-deactivate.csh \ | |
| 63 | 66 | dev/conda.recipe/build.sh \ |
| 64 | 67 | dev/conda.recipe/conda_build_config.yaml \ |
| 65 | 68 | dev/conda.recipe/meta.yaml.in | ... | ... |
dev/00-gammalib-activate.csh
0 → 100644
dev/01-gammalib-deactivate.csh
0 → 100644
doc/source/admin/download.rst
| ... | ... | @@ -4,13 +4,14 @@ Download |
| 4 | 4 | ======== |
| 5 | 5 | |
| 6 | 6 | The latest GammaLib release is |
| 7 | -`gammalib-1.7.3 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.3.tar.gz>`_ | |
| 8 | -(2 December 2020). | |
| 7 | +`gammalib-1.7.4 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.4.tar.gz>`_ | |
| 8 | +(28 April 2021). | |
| 9 | 9 | |
| 10 | 10 | Below a list of older GammaLib releases. Please read the :ref:`release` to |
| 11 | 11 | learn more about new features and corrected bugs in a given release. |
| 12 | 12 | |
| 13 | 13 | To download an older release, click on the corresponding release number: |
| 14 | +`1.7.3 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.3.tar.gz>`_ | |
| 14 | 15 | `1.7.2 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.2.tar.gz>`_ |
| 15 | 16 | `1.7.1 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.1.tar.gz>`_ |
| 16 | 17 | `1.7.0 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.0.tar.gz>`_ | ... | ... |
doc/source/admin/release_history/1.7.4.rst
0 → 100644
| 1 | +.. _1.7.4: | |
| 2 | + | |
| 3 | +GammaLib 1.7.4 release (28 April 2021) | |
| 4 | +====================================== | |
| 5 | + | |
| 6 | +Overview | |
| 7 | +-------- | |
| 8 | + | |
| 9 | +GammaLib 1.7.4 is a bug fix release for GammaLib release 1.7. | |
| 10 | + | |
| 11 | +The following bugs were fixed: | |
| 12 | + | |
| 13 | +* [`3628 <https://cta-redmine.irap.omp.eu/issues/3628>`_] - | |
| 14 | + Add conda activation and deactivation scripts for csh | |
| 15 | +* [`3625 <https://cta-redmine.irap.omp.eu/issues/3625>`_] - | |
| 16 | + Speed-up Monte-Carlo sampling of exponentially cut-off power law models | ... | ... |
doc/source/admin/release_history/index.rst
inst/cta/src/GCTAModelIrfBackground.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | 2 | * GCTAModelIrfBackground.cpp - CTA IRF background model class * |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2014-2020 by Juergen Knoedlseder * | |
| 4 | + * copyright (C) 2014-2021 by Juergen Knoedlseder * | |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| 6 | 6 | * * |
| 7 | 7 | * This program is free software: you can redistribute it and/or modify * | ... | ... |
src/model/GModelSpectralComposite.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | 2 | * GModelSpectralComposite.cpp - Spectral composite model class * |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2016-2017 by Michael Mayer * | |
| 4 | + * copyright (C) 2016-2021 by Michael Mayer * | |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| 6 | 6 | * * |
| 7 | 7 | * This program is free software: you can redistribute it and/or modify * |
| ... | ... | @@ -54,6 +54,7 @@ const GModelSpectralRegistry g_spectral_comp_registry(&g_spectral_comp_seed); |
| 54 | 54 | /* __ Coding definitions _________________________________________________ */ |
| 55 | 55 | |
| 56 | 56 | /* __ Debug definitions __________________________________________________ */ |
| 57 | +//#define G_DEBUG_UPDATE_MC_CACHE //!< Debug MC cache update | |
| 57 | 58 | |
| 58 | 59 | |
| 59 | 60 | /*========================================================================== |
| ... | ... | @@ -698,7 +699,7 @@ void GModelSpectralComposite::copy_members(const GModelSpectralComposite& model) |
| 698 | 699 | // Clear parameters and copy the pointers from the clone |
| 699 | 700 | m_pars.clear(); |
| 700 | 701 | for (int i = 0; i < model.components(); ++i) { |
| 701 | - | |
| 702 | + | |
| 702 | 703 | // Clone spectral component |
| 703 | 704 | m_spectral.push_back(model.m_spectral[i]->clone()); |
| 704 | 705 | |
| ... | ... | @@ -808,6 +809,17 @@ void GModelSpectralComposite::update_mc_cache(const GEnergy& emin, |
| 808 | 809 | |
| 809 | 810 | } //endfor: looped over spectral components |
| 810 | 811 | |
| 812 | + // Debug option: signal update | |
| 813 | + #if defined(G_DEBUG_UPDATE_MC_CACHE) | |
| 814 | + std::cout << "GModelSpectralComposite::update_mc_cache("; | |
| 815 | + std::cout << emin.print() << "," << emax.print() << "):"; | |
| 816 | + std::cout << " flux=" << m_mc_flux; | |
| 817 | + for (int i = 0; i < m_spectral.size(); ++i) { | |
| 818 | + std::cout << " prob[" << i << "]=" << m_mc_probs[i]; | |
| 819 | + } | |
| 820 | + std::cout << std::endl; | |
| 821 | + #endif | |
| 822 | + | |
| 811 | 823 | } // endif: emin and emax have changed |
| 812 | 824 | |
| 813 | 825 | // Return | ... | ... |
src/model/GModelSpectralExpInvPlaw.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | 2 | * GModelSpectralExpInvPlaw.cpp - Exponential cut off power law model * |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2016 by Alexander Ziegler * | |
| 4 | + * copyright (C) 2016-2021 by Alexander Ziegler * | |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| 6 | 6 | * * |
| 7 | 7 | * This program is free software: you can redistribute it and/or modify * |
| ... | ... | @@ -57,6 +57,7 @@ const GModelSpectralRegistry g_spectral_einvplaw_registry1(&g_spectral_einvpla |
| 57 | 57 | /* __ Coding definitions _________________________________________________ */ |
| 58 | 58 | |
| 59 | 59 | /* __ Debug definitions __________________________________________________ */ |
| 60 | +//#define G_DEBUG_MC //!< Debug Monte-Carlo sampling | |
| 60 | 61 | |
| 61 | 62 | |
| 62 | 63 | /*========================================================================== |
| ... | ... | @@ -402,7 +403,7 @@ double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, |
| 402 | 403 | // Compile option: Check for NaN/Inf |
| 403 | 404 | #if defined(G_NAN_CHECK) |
| 404 | 405 | if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { |
| 405 | - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; | |
| 406 | + std::cout << "*** ERROR: GModelSpectralExpInvPlaw::eval"; | |
| 406 | 407 | std::cout << "(srcEng=" << srcEng; |
| 407 | 408 | std::cout << ", srcTime=" << srcTime << "):"; |
| 408 | 409 | std::cout << " NaN/Inf encountered"; |
| ... | ... | @@ -436,11 +437,11 @@ double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, |
| 436 | 437 | * The integration is done numerically. |
| 437 | 438 | ***************************************************************************/ |
| 438 | 439 | double GModelSpectralExpInvPlaw::flux(const GEnergy& emin, |
| 439 | - const GEnergy& emax) const | |
| 440 | + const GEnergy& emax) const | |
| 440 | 441 | { |
| 441 | 442 | // Initialise flux |
| 442 | 443 | double flux = 0.0; |
| 443 | - | |
| 444 | + | |
| 444 | 445 | // Compute only if integration range is valid |
| 445 | 446 | if (emin < emax) { |
| 446 | 447 | |
| ... | ... | @@ -482,11 +483,11 @@ double GModelSpectralExpInvPlaw::flux(const GEnergy& emin, |
| 482 | 483 | * The integration is done numerically. |
| 483 | 484 | ***************************************************************************/ |
| 484 | 485 | double GModelSpectralExpInvPlaw::eflux(const GEnergy& emin, |
| 485 | - const GEnergy& emax) const | |
| 486 | + const GEnergy& emax) const | |
| 486 | 487 | { |
| 487 | 488 | // Initialise flux |
| 488 | 489 | double eflux = 0.0; |
| 489 | - | |
| 490 | + | |
| 490 | 491 | // Compute only if integration range is valid |
| 491 | 492 | if (emin < emax) { |
| 492 | 493 | |
| ... | ... | @@ -531,9 +532,9 @@ double GModelSpectralExpInvPlaw::eflux(const GEnergy& emin, |
| 531 | 532 | * on an acceptance fraction that is computed from the exponential cut off. |
| 532 | 533 | ***************************************************************************/ |
| 533 | 534 | GEnergy GModelSpectralExpInvPlaw::mc(const GEnergy& emin, |
| 534 | - const GEnergy& emax, | |
| 535 | - const GTime& time, | |
| 536 | - GRan& ran) const | |
| 535 | + const GEnergy& emax, | |
| 536 | + const GTime& time, | |
| 537 | + GRan& ran) const | |
| 537 | 538 | { |
| 538 | 539 | // Throw an exception if energy range is invalid |
| 539 | 540 | if (emin >= emax) { |
| ... | ... | @@ -549,11 +550,25 @@ GEnergy GModelSpectralExpInvPlaw::mc(const GEnergy& emin, |
| 549 | 550 | |
| 550 | 551 | // Initialise energy |
| 551 | 552 | double eng; |
| 552 | - | |
| 553 | + | |
| 553 | 554 | // Initialse acceptance fraction |
| 554 | 555 | double acceptance_fraction; |
| 555 | 556 | double inv_ecut = m_lambda.value(); |
| 556 | 557 | |
| 558 | + // Compute random number generator normalisation for speed-up of the | |
| 559 | + // computations. Normally it should be sufficient to always use the | |
| 560 | + // minimum energy for the normalisation, yet for the case that the model | |
| 561 | + // grows with energy we better also test the maximum energy, and we | |
| 562 | + // just use the larger of both normalisations | |
| 563 | + double norm_emin = std::exp(-emin.MeV() * inv_ecut); | |
| 564 | + double norm_emax = std::exp(-emax.MeV() * inv_ecut); | |
| 565 | + double norm = (norm_emin > norm_emax) ? norm_emin : norm_emax; | |
| 566 | + | |
| 567 | + // Debug option: initialise number of samples | |
| 568 | + #if defined(G_DEBUG_MC) | |
| 569 | + int samples = 0; | |
| 570 | + #endif | |
| 571 | + | |
| 557 | 572 | // Use rejection method to draw a random energy. We first draw |
| 558 | 573 | // analytically from a power law, and then compare the power law |
| 559 | 574 | // at the drawn energy to the exponentially cut off function. This |
| ... | ... | @@ -561,6 +576,11 @@ GEnergy GModelSpectralExpInvPlaw::mc(const GEnergy& emin, |
| 561 | 576 | // a uniform random number is <= the acceptance fraction. |
| 562 | 577 | do { |
| 563 | 578 | |
| 579 | + // Debug option: increment number of samples | |
| 580 | + #if defined(G_DEBUG_MC) | |
| 581 | + samples++; | |
| 582 | + #endif | |
| 583 | + | |
| 564 | 584 | // Get uniform random number |
| 565 | 585 | double u = ran.uniform(); |
| 566 | 586 | |
| ... | ... | @@ -579,15 +599,23 @@ GEnergy GModelSpectralExpInvPlaw::mc(const GEnergy& emin, |
| 579 | 599 | else { |
| 580 | 600 | eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin); |
| 581 | 601 | } |
| 582 | - | |
| 602 | + | |
| 583 | 603 | // Compute acceptance fraction |
| 584 | 604 | acceptance_fraction = std::exp(-eng * inv_ecut); |
| 585 | 605 | |
| 586 | - } while (ran.uniform() > acceptance_fraction); | |
| 587 | - | |
| 606 | + } while (ran.uniform() * norm > acceptance_fraction); | |
| 607 | + | |
| 588 | 608 | // Set energy |
| 589 | 609 | energy.MeV(eng); |
| 590 | 610 | |
| 611 | + // Debug option: write result | |
| 612 | + #if defined(G_DEBUG_MC) | |
| 613 | + std::cout << "GModelSpectralExpInvPlaw::mc("; | |
| 614 | + std::cout << emin.print() << "," << emax.print() << "," << time.print() << "):"; | |
| 615 | + std::cout << " energy=" << energy.print(); | |
| 616 | + std::cout << " samples=" << samples << std::endl; | |
| 617 | + #endif | |
| 618 | + | |
| 591 | 619 | // Return energy |
| 592 | 620 | return energy; |
| 593 | 621 | } |
| ... | ... | @@ -841,7 +869,7 @@ void GModelSpectralExpInvPlaw::update_eval_cache(const GEnergy& energy) const |
| 841 | 869 | double index = m_index.value(); |
| 842 | 870 | double lambda = m_lambda.value(); |
| 843 | 871 | double pivot = m_pivot.value(); |
| 844 | - | |
| 872 | + | |
| 845 | 873 | // If the energy or one of the parameters index, lambda or pivot |
| 846 | 874 | // energy has changed then recompute the cache |
| 847 | 875 | if ((m_last_energy != energy) || | ... | ... |
src/model/GModelSpectralExpPlaw.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | 2 | * GModelSpectralExpPlaw.cpp - Exponential cut off power law model * |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2010-2016 by Juergen Knoedlseder * | |
| 4 | + * copyright (C) 2010-2021 by Juergen Knoedlseder * | |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| 6 | 6 | * * |
| 7 | 7 | * This program is free software: you can redistribute it and/or modify * |
| ... | ... | @@ -64,6 +64,7 @@ const GModelSpectralRegistry g_spectral_eplaw_registry2(&g_spectral_eplaw_seed2) |
| 64 | 64 | /* __ Coding definitions _________________________________________________ */ |
| 65 | 65 | |
| 66 | 66 | /* __ Debug definitions __________________________________________________ */ |
| 67 | +//#define G_DEBUG_MC //!< Debug Monte-Carlo sampling | |
| 67 | 68 | |
| 68 | 69 | |
| 69 | 70 | /*========================================================================== |
| ... | ... | @@ -410,7 +411,7 @@ double GModelSpectralExpPlaw::flux(const GEnergy& emin, |
| 410 | 411 | { |
| 411 | 412 | // Initialise flux |
| 412 | 413 | double flux = 0.0; |
| 413 | - | |
| 414 | + | |
| 414 | 415 | // Compute only if integration range is valid |
| 415 | 416 | if (emin < emax) { |
| 416 | 417 | |
| ... | ... | @@ -456,7 +457,7 @@ double GModelSpectralExpPlaw::eflux(const GEnergy& emin, |
| 456 | 457 | { |
| 457 | 458 | // Initialise flux |
| 458 | 459 | double eflux = 0.0; |
| 459 | - | |
| 460 | + | |
| 460 | 461 | // Compute only if integration range is valid |
| 461 | 462 | if (emin < emax) { |
| 462 | 463 | |
| ... | ... | @@ -519,11 +520,25 @@ GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, |
| 519 | 520 | |
| 520 | 521 | // Initialise energy |
| 521 | 522 | double eng; |
| 522 | - | |
| 523 | - // Initialse acceptance fraction | |
| 523 | + | |
| 524 | + // Initialise acceptance fraction | |
| 524 | 525 | double acceptance_fraction; |
| 525 | 526 | double inv_ecut = 1.0 / m_ecut.value(); |
| 526 | 527 | |
| 528 | + // Compute random number generator normalisation for speed-up of the | |
| 529 | + // computations. Normally it should be sufficient to always use the | |
| 530 | + // minimum energy for the normalisation, yet for the case that the model | |
| 531 | + // grows with energy we better also test the maximum energy, and we | |
| 532 | + // just use the larger of both normalisations | |
| 533 | + double norm_emin = std::exp(-emin.MeV() * inv_ecut); | |
| 534 | + double norm_emax = std::exp(-emax.MeV() * inv_ecut); | |
| 535 | + double norm = (norm_emin > norm_emax) ? norm_emin : norm_emax; | |
| 536 | + | |
| 537 | + // Debug option: initialise number of samples | |
| 538 | + #if defined(G_DEBUG_MC) | |
| 539 | + int samples = 0; | |
| 540 | + #endif | |
| 541 | + | |
| 527 | 542 | // Use rejection method to draw a random energy. We first draw |
| 528 | 543 | // analytically from a power law, and then compare the power law |
| 529 | 544 | // at the drawn energy to the exponentially cut off function. This |
| ... | ... | @@ -531,6 +546,11 @@ GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, |
| 531 | 546 | // a uniform random number is <= the acceptance fraction. |
| 532 | 547 | do { |
| 533 | 548 | |
| 549 | + // Debug option: increment number of samples | |
| 550 | + #if defined(G_DEBUG_MC) | |
| 551 | + samples++; | |
| 552 | + #endif | |
| 553 | + | |
| 534 | 554 | // Get uniform random number |
| 535 | 555 | double u = ran.uniform(); |
| 536 | 556 | |
| ... | ... | @@ -549,15 +569,23 @@ GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, |
| 549 | 569 | else { |
| 550 | 570 | eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin); |
| 551 | 571 | } |
| 552 | - | |
| 572 | + | |
| 553 | 573 | // Compute acceptance fraction |
| 554 | 574 | acceptance_fraction = std::exp(-eng * inv_ecut); |
| 555 | 575 | |
| 556 | - } while (ran.uniform() > acceptance_fraction); | |
| 557 | - | |
| 576 | + } while (ran.uniform() * norm > acceptance_fraction); | |
| 577 | + | |
| 558 | 578 | // Set energy |
| 559 | 579 | energy.MeV(eng); |
| 560 | 580 | |
| 581 | + // Debug option: write result | |
| 582 | + #if defined(G_DEBUG_MC) | |
| 583 | + std::cout << "GModelSpectralExpPlaw::mc("; | |
| 584 | + std::cout << emin.print() << "," << emax.print() << "," << time.print() << "):"; | |
| 585 | + std::cout << " energy=" << energy.print(); | |
| 586 | + std::cout << " samples=" << samples << std::endl; | |
| 587 | + #endif | |
| 588 | + | |
| 561 | 589 | // Return energy |
| 562 | 590 | return energy; |
| 563 | 591 | } |
| ... | ... | @@ -632,7 +660,7 @@ void GModelSpectralExpPlaw::write(GXmlElement& xml) const |
| 632 | 660 | /***********************************************************************//** |
| 633 | 661 | * @brief Print model information |
| 634 | 662 | * |
| 635 | - * @param[in] chatter Chattiness (defaults to NORMAL). | |
| 663 | + * @param[in] chatter Chattiness. | |
| 636 | 664 | * @return String containing model information. |
| 637 | 665 | ***************************************************************************/ |
| 638 | 666 | std::string GModelSpectralExpPlaw::print(const GChatter& chatter) const |
| ... | ... | @@ -810,7 +838,7 @@ void GModelSpectralExpPlaw::update_eval_cache(const GEnergy& energy) const |
| 810 | 838 | double index = m_index.value(); |
| 811 | 839 | double ecut = m_ecut.value(); |
| 812 | 840 | double pivot = m_pivot.value(); |
| 813 | - | |
| 841 | + | |
| 814 | 842 | // If the energy or one of the parameters index, cut-off or pivot |
| 815 | 843 | // energy has changed then recompute the cache |
| 816 | 844 | if ((m_last_energy != energy) || | ... | ... |
src/model/GModelSpectralSuperExpPlaw.cpp
| ... | ... | @@ -67,6 +67,7 @@ const GModelSpectralRegistry g_spectral_seplaw_registry2(&g_spectral_seplaw_ |
| 67 | 67 | /* __ Coding definitions _________________________________________________ */ |
| 68 | 68 | |
| 69 | 69 | /* __ Debug definitions __________________________________________________ */ |
| 70 | +//#define G_DEBUG_MC //!< Debug Monte-Carlo sampling | |
| 70 | 71 | |
| 71 | 72 | |
| 72 | 73 | /*========================================================================== |
| ... | ... | @@ -435,7 +436,7 @@ double GModelSpectralSuperExpPlaw::flux(const GEnergy& emin, |
| 435 | 436 | { |
| 436 | 437 | // Initialise flux |
| 437 | 438 | double flux = 0.0; |
| 438 | - | |
| 439 | + | |
| 439 | 440 | // Compute only if integration range is valid |
| 440 | 441 | if (emin < emax) { |
| 441 | 442 | |
| ... | ... | @@ -481,7 +482,7 @@ double GModelSpectralSuperExpPlaw::eflux(const GEnergy& emin, |
| 481 | 482 | { |
| 482 | 483 | // Initialise flux |
| 483 | 484 | double eflux = 0.0; |
| 484 | - | |
| 485 | + | |
| 485 | 486 | // Compute only if integration range is valid |
| 486 | 487 | if (emin < emax) { |
| 487 | 488 | |
| ... | ... | @@ -544,11 +545,25 @@ GEnergy GModelSpectralSuperExpPlaw::mc(const GEnergy& emin, |
| 544 | 545 | |
| 545 | 546 | // Initialise energy |
| 546 | 547 | double eng; |
| 547 | - | |
| 548 | - // Initialse acceptance fraction | |
| 548 | + | |
| 549 | + // Initialise acceptance fraction | |
| 549 | 550 | double acceptance_fraction; |
| 550 | 551 | double inv_ecut = 1.0 / m_ecut.value(); |
| 551 | 552 | |
| 553 | + // Compute random number generator normalisation for speed-up of the | |
| 554 | + // computations. Normally it should be sufficient to always use the | |
| 555 | + // minimum energy for the normalisation, yet for the case that the model | |
| 556 | + // grows with energy we better also test the maximum energy, and we | |
| 557 | + // just use the larger of both normalisations | |
| 558 | + double norm_emin = std::exp(- std::pow(emin.MeV() * inv_ecut, m_index2.value())); | |
| 559 | + double norm_emax = std::exp(- std::pow(emax.MeV() * inv_ecut, m_index2.value())); | |
| 560 | + double norm = (norm_emin > norm_emax) ? norm_emin : norm_emax; | |
| 561 | + | |
| 562 | + // Debug option: initialise number of samples | |
| 563 | + #if defined(G_DEBUG_MC) | |
| 564 | + int samples = 0; | |
| 565 | + #endif | |
| 566 | + | |
| 552 | 567 | // Use rejection method to draw a random energy. We first draw |
| 553 | 568 | // analytically from a power law, and then compare the power law |
| 554 | 569 | // at the drawn energy to the exponentially cut off function. This |
| ... | ... | @@ -556,6 +571,11 @@ GEnergy GModelSpectralSuperExpPlaw::mc(const GEnergy& emin, |
| 556 | 571 | // a uniform random number is <= the acceptance fraction. |
| 557 | 572 | do { |
| 558 | 573 | |
| 574 | + // Debug option: increment number of samples | |
| 575 | + #if defined(G_DEBUG_MC) | |
| 576 | + samples++; | |
| 577 | + #endif | |
| 578 | + | |
| 559 | 579 | // Get uniform random number |
| 560 | 580 | double u = ran.uniform(); |
| 561 | 581 | |
| ... | ... | @@ -574,15 +594,23 @@ GEnergy GModelSpectralSuperExpPlaw::mc(const GEnergy& emin, |
| 574 | 594 | else { |
| 575 | 595 | eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin); |
| 576 | 596 | } |
| 577 | - | |
| 597 | + | |
| 578 | 598 | // Compute acceptance fraction |
| 579 | 599 | acceptance_fraction = std::exp(- std::pow(eng * inv_ecut, m_index2.value())); |
| 580 | 600 | |
| 581 | - } while (ran.uniform() > acceptance_fraction); | |
| 582 | - | |
| 601 | + } while (ran.uniform() * norm > acceptance_fraction); | |
| 602 | + | |
| 583 | 603 | // Set energy |
| 584 | 604 | energy.MeV(eng); |
| 585 | 605 | |
| 606 | + // Debug option: write result | |
| 607 | + #if defined(G_DEBUG_MC) | |
| 608 | + std::cout << "GModelSpectralSuperExpPlaw::mc("; | |
| 609 | + std::cout << emin.print() << "," << emax.print() << "," << time.print() << "):"; | |
| 610 | + std::cout << " energy=" << energy.print(); | |
| 611 | + std::cout << " samples=" << samples << std::endl; | |
| 612 | + #endif | |
| 613 | + | |
| 586 | 614 | // Return energy |
| 587 | 615 | return energy; |
| 588 | 616 | } |
| ... | ... | @@ -750,7 +778,7 @@ void GModelSpectralSuperExpPlaw::init_members(void) |
| 750 | 778 | m_index2.name("Index2"); |
| 751 | 779 | m_index2.scale(1.0); |
| 752 | 780 | m_index2.value(1.0); // default: -2.0 |
| 753 | - m_index2.range(0.0,5.0); // range: [-10,+10] | |
| 781 | + m_index2.range(0.0,5.0); // range: [0,5] | |
| 754 | 782 | m_index2.free(); |
| 755 | 783 | m_index2.gradient(0.0); |
| 756 | 784 | m_index2.has_grad(true); |
| ... | ... | @@ -857,7 +885,7 @@ void GModelSpectralSuperExpPlaw::update_eval_cache(const GEnergy& energy) const |
| 857 | 885 | double ecut = m_ecut.value(); |
| 858 | 886 | double pivot = m_pivot.value(); |
| 859 | 887 | double index2 = m_index2.value(); |
| 860 | - | |
| 888 | + | |
| 861 | 889 | // If the energy or one of the parameters index1, index2, cut-off or pivot |
| 862 | 890 | // energy has changed then recompute the cache |
| 863 | 891 | if ((m_last_energy != energy) || | ... | ... |