From 030a746c04abc8d91e82dc4717288f9fdc5d1ec2 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 7 Apr 2017 20:25:01 +0200 Subject: [PATCH] Integrate code --- AUTHORS | 1 + ChangeLog | 3 ++- NEWS | 7 ++++++- inst/cta/src/GCTAModelAeffBackground.cpp | 87 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------- inst/cta/src/GCTAModelIrfBackground.cpp | 17 ++++++++++++----- 5 files changed, 80 insertions(+), 35 deletions(-) diff --git a/AUTHORS b/AUTHORS index dc30cac..033feed 100644 --- a/AUTHORS +++ b/AUTHORS @@ -4,6 +4,7 @@ Joshua Cardenzana (jcardenzana@irap.omp.eu) Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) Alex Donath (Axel.Donath@mpi-hd.mpg.de) +Stefan Eschbach (stefan.eschbach@fau.de) Florent Forest (florent.forest@isae.fr) Lucie Gerard (lucie.gerard@desy.de) Tarek Hassan (thassan@ifae.es) diff --git a/ChangeLog b/ChangeLog index 4ef0990..261056f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-03-27 +2017-04-07 * Version 1.3.0 released ========================= + Avoid hang-up of GCTAModelAeffBackground::mc() method (#1988) Enhance precision of GCTAEdisp2D::compute_ebounds_src Add conversion methods between GEbounds and GEnergies (#1977) Add GModelSpectralSmoothBrokenPlaw class (#1948) diff --git a/NEWS b/NEWS index 8e49fff..3bcdb1f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.3.0.dev1 -27 March 2017 +7 April 2017 1. Introduction @@ -160,6 +160,11 @@ method. In addition, the validity range of the table is now stored in members of the class, and a non-zero energy dispersion is only returned for valid elements of the response table. +The GCTAModelAeffBackground::mc() could get stuck in the rejection loop if the +effective area returned is always zero. A counter has been added that leads to +the exit of the loop in case that 1000 subsequent zero values are returned +(#1988). + 21. COMPTEL interface --------------------- diff --git a/inst/cta/src/GCTAModelAeffBackground.cpp b/inst/cta/src/GCTAModelAeffBackground.cpp index dbf18f2..ba139b7 100644 --- a/inst/cta/src/GCTAModelAeffBackground.cpp +++ b/inst/cta/src/GCTAModelAeffBackground.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelAeffBackground.cpp - CTA Aeff background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2015-2016 by Michael Mayer * + * copyright (C) 2015-2017 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -509,13 +509,14 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, // Retrieve CTA observation const GCTAObservation* cta = dynamic_cast(&obs); if (cta == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); + std::string msg = "Specified observation is not a CTA " + "observation.\n" + obs.print(); throw GException::invalid_argument(G_MC, msg); } // Get pointer on CTA IRF response - const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); + const GCTAResponseIrf* rsp = + dynamic_cast(cta->response()); if (rsp == NULL) { std::string msg = "Specified observation does not contain" " an IRF response.\n" + obs.print(); @@ -534,7 +535,8 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, } // Retrieve event list to access the ROI, energy boundaries and GTIs - const GCTAEventList* events = dynamic_cast(obs.events()); + const GCTAEventList* events = + dynamic_cast(obs.events()); if (events == NULL) { std::string msg = "No CTA event list found in observation.\n" + obs.print(); @@ -559,19 +561,21 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, // Set up spectral model to draw random energies from. Here we use // a fixed energy sampling for an instance of GModelSpectralNodes. // This is analogous to to the GCTAModelIrfBackground::mc method. - GEbounds spectral_ebounds = GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true); + // We make sure that only non-negative nodes get appended. + GEbounds spectral_ebounds = + GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true); GModelSpectralNodes spectral; for (int i = 0; i < spectral_ebounds.size(); ++i) { - GEnergy energy = spectral_ebounds.elogmean(i); - double intensity = aeff_integral(obs, energy.log10TeV()); - double norm = m_spectral->eval(energy, events->tstart()); - if(norm * intensity > 0){ - spectral.append(energy, norm * intensity); - } + GEnergy energy = spectral_ebounds.elogmean(i); + double intensity = aeff_integral(obs, energy.log10TeV()); + double norm = m_spectral->eval(energy, events->tstart()); + double arg = norm * intensity; + if (arg > 0.0) { + spectral.append(energy, arg); + } } - - - // Loop over all energy boundaries + + // Loop over all energy bins for (int ieng = 0; ieng < ebounds.size(); ++ieng) { // Compute the background rate in model within the energy @@ -586,6 +590,11 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, std::cout << "rate=" << rate << " cts/s)" << std::endl; #endif + // If the rate is not positive then skip this energy bins + if (rate <= 0.0) { + continue; + } + // Loop over all good time intervals for (int itime = 0; itime < gti.size(); ++itime) { @@ -633,17 +642,24 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, ran); // Get maximum effective area for rejection method - double max_aeff = aeff->max(energy.log10TeV(), pnt.zenith(), pnt.azimuth(), false); + double max_aeff = aeff->max(energy.log10TeV(), pnt.zenith(), + pnt.azimuth(), false); + + // Skip event if the maximum effective area is not positive + if (max_aeff <= 0.0) { + continue; + } // Initialise randomised coordinates double offset = 0.0; double phi = 0.0; - // Initialise acceptance fraction for rejection method + // Initialise acceptance fraction and counter of zeros for + // rejection method double acceptance_fraction = 0.0; + int zeros = 0; // Start rejection method loop - int numberoftrys=0; do { // Throw random offset and azimuth angle in @@ -653,24 +669,39 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, phi = ran.uniform() * gammalib::twopi; // Compute function value at this offset angle - double value = (*aeff)(energy.log10TeV(), offset, phi, pnt.zenith(), pnt.azimuth(), false); - - // std::cout << "value: " << value << ", max_aeff: " << max_aeff << std::endl; - if(numberoftrys==1000&&value==0){ - break; - } + double value = (*aeff)(energy.log10TeV(), offset, phi, + pnt.zenith(), pnt.azimuth(), + false); + + // If the value is not positive then increment the + // zeros counter and fall through. The counter assures + // that this loop does not lock up. + if (value <= 0.0) { + zeros++; + continue; + } + + // Value is non-zero so reset the zeros counter + zeros = 0; + // Compute acceptance fraction acceptance_fraction = value / max_aeff; - - numberoftrys++; - } while (ran.uniform() > acceptance_fraction); + } while ((ran.uniform() > acceptance_fraction) && + (zeros < 1000)); + + // If the zeros counter is non-zero then the loop was + // exited due to exhaustion and the event is skipped + if (zeros > 0) { + continue; + } // Convert CTA pointing direction in instrument system GCTAInstDir mc_dir(pnt.dir()); // Rotate pointing direction by offset and azimuth angle - mc_dir.dir().rotate_deg(phi * gammalib::rad2deg, offset * gammalib::rad2deg); + mc_dir.dir().rotate_deg(phi * gammalib::rad2deg, + offset * gammalib::rad2deg); // Compute DETX and DETY coordinates double detx(0.0); diff --git a/inst/cta/src/GCTAModelIrfBackground.cpp b/inst/cta/src/GCTAModelIrfBackground.cpp index 7bb4f35..25c1e29 100644 --- a/inst/cta/src/GCTAModelIrfBackground.cpp +++ b/inst/cta/src/GCTAModelIrfBackground.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelIrfBackground.cpp - CTA IRF background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -551,13 +551,14 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation& obs, GRan& ran) co // Retrieve CTA observation const GCTAObservation* cta = dynamic_cast(&obs); if (cta == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); + std::string msg = "Specified observation is not a CTA " + "observation.\n" + obs.print(); throw GException::invalid_argument(G_MC, msg); } // Get pointer on CTA IRF response - const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); + const GCTAResponseIrf* rsp = + dynamic_cast(cta->response()); if (rsp == NULL) { std::string msg = "Specified observation does not contain" " an IRF response.\n" + obs.print(); @@ -576,7 +577,8 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation& obs, GRan& ran) co } // Retrieve event list to access the ROI, energy boundaries and GTIs - const GCTAEventList* events = dynamic_cast(obs.events()); + const GCTAEventList* events = + dynamic_cast(obs.events()); if (events == NULL) { std::string msg = "No CTA event list found in observation.\n" + obs.print(); @@ -618,6 +620,11 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation& obs, GRan& ran) co std::cout << "rate=" << rate << " cts/s)" << std::endl; #endif + // If the rate is not positive then skip this energy bins + if (rate <= 0.0) { + continue; + } + // Loop over all good time intervals for (int itime = 0; itime < gti.size(); ++itime) { -- libgit2 0.21.4