Integrate code

Jürgen Knödlseder
1 parent d69e7faa
Showing 5 changed files with 80 additions and 35 deletions
AUTHORS
ChangeLog
NEWS
inst/cta/src/GCTAModelAeffBackground.cpp
inst/cta/src/GCTAModelIrfBackground.cpp
@@ -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)
-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)
 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
 ---------------------
 /***************************************************************************
  *       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&amp; obs,
         // Retrieve CTA observation
         const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&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<const GCTAResponseIrf*>(cta->response());
+        const GCTAResponseIrf* rsp =
+              dynamic_cast<const GCTAResponseIrf*>(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&amp; obs,
         }
  
         // Retrieve event list to access the ROI, energy boundaries and GTIs
-        const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
+        const GCTAEventList* events =
+              dynamic_cast<const GCTAEventList*>(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&amp; 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&amp; 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&amp; 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&amp; 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);
 /***************************************************************************
  *       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&amp; obs, GRan&amp; ran) co
         // Retrieve CTA observation
         const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&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<const GCTAResponseIrf*>(cta->response());
+        const GCTAResponseIrf* rsp =
+              dynamic_cast<const GCTAResponseIrf*>(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&amp; obs, GRan&amp; ran) co
         }
  
         // Retrieve event list to access the ROI, energy boundaries and GTIs
-        const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
+        const GCTAEventList* events =
+              dynamic_cast<const GCTAEventList*>(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&amp; obs, GRan&amp; 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) {