Add maxrate parameter to throw an exception in case that the source event rate is too high (#807).

Jürgen Knödlseder
1 parent 86d4c073
Showing 5 changed files with 97 additions and 19 deletions
doc/source/reference_manual/ctobssim.rst
pyext/ctobssim.i
src/ctobssim/ctobssim.cpp
src/ctobssim/ctobssim.hpp
src/ctobssim/ctobssim.par
@@ -9,9 +9,9 @@ Simulates CTA events.
 Synopsis
 --------
-Generate photon events from astrophysical sources and process those photons
-according to the specified instrument response functions and generate
-instrumental background events.
+Generate photon events from astrophysical sources and background events 
+from an instrumental background model.
+
 General parameters
 ------------------
@@ -24,7 +24,7 @@ General parameters
 ``inmodel [file]``
     XML file that describes the astrophysical sources and the instrumental
-    background.
+    background model.
 ``outevents [file]``
     Output event list or observation definition XML file.
@@ -39,7 +39,10 @@ General parameters
     Instrumental response function.
 ``(seed = 1) [integer]``
-    Integer seed value to be used for Monte Carlo simulations.
+    Integer seed value to be used for Monte Carlo simulations. Keep this 
+    parameter at the same value for repeatable simulations, or increment 
+    this value for subsequent runs if non-repeatable simulations are
+    required.
 ``ra [real]``
     Right Ascension of CTA pointing (J2000, in degrees).
@@ -57,10 +60,10 @@ General parameters
     Mission elapsed stop time of observation (in seconds).
 ``emin [real]``
-    Lower energy limit of events (in TeV).
+    Lower energy limit of simulated events (in TeV).
 ``emax [real]``
-    Upper energy limit of events (in TeV).
+    Upper energy limit of simulated events (in TeV).
 ``(edisp = no) [boolean]``
     Apply energy dispersion?
@@ -68,6 +71,13 @@ General parameters
 ``(deadc = 0.95) [real]``
     Average deadtime correction factor.
+``(maxrate = 1.0e6) [real]``
+    Maximum photon rate for source models. Source models that exceed this
+    maximum photon rate will lead to an exception as very likely the
+    specified model normalisation is too large (probably due to the
+    a misinterpretation of units). Note that ctools specifies intensity
+    units per MeV.
+
 Standard parameters
 -------------------
 /***************************************************************************
  *                  ctobssim - Observation simulator tool                  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -48,6 +48,8 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
+    const double&        max_rate(void) const;
+    void                 max_rate(const double& max_rate);
 };
@@ -494,6 +494,7 @@ void ctobssim::init_members(void)
     m_prefix.clear();
     m_seed        = 1;
     m_apply_edisp = false;
+    m_max_rate    = 1.0e6;
     // Initialise protected members
     m_rans.clear();
@@ -524,6 +525,7 @@ void ctobssim::copy_members(const ctobssim&amp; app)
     m_prefix      = app.m_prefix;
     m_seed        = app.m_seed;
     m_apply_edisp = app.m_apply_edisp;
+    m_max_rate    = app.m_max_rate;
     // Copy protected members
     m_area             = app.m_area;
@@ -577,6 +579,7 @@ void ctobssim::get_parameters(void)
     // Get other parameters
     m_seed        = (*this)["seed"].integer();
     m_apply_edisp = (*this)["edisp"].boolean();
+    m_max_rate    = (*this)["maxrate"].real();
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
@@ -768,6 +771,36 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
                         }
                         GTime tslice(duration, "sec");
+                        // If photon rate exceeds the maximum photon rate
+                        // then throw an exception
+                        if (rate > m_max_rate) {
+                            std::string modnam = (model->name().length() > 0) ?
+                                                 model->name() : "Unknown";
+                            std::string msg    = "Photon rate "+
+                                                 gammalib::str(rate)+
+                                                 " photons/sec for model \""+
+                                                 modnam+"\" exceeds maximum"
+                                                 " allowed photon rate of "+
+                                                 gammalib::str(m_max_rate)+
+                                                 " photons/sec. Please check"
+                                                 " the model parameters for"
+                                                 " model \""+modnam+"\" or"
+                                                 " increase the value of the"
+                                                 " hidden \"maxrate\""
+                                                 " parameter.";
+                            throw GException::invalid_value(G_SIMULATE_SOURCE, msg);
+                        }
+
+                        // Dump length of time slice and rate
+                        if (logExplicit()) {
+                            *wrklog << gammalib::parformat("Photon rate", indent);
+                            *wrklog << rate << " photons/sec";
+                            if (model->name().length() > 0) {
+                                *wrklog << " [" << model->name() << "]";
+                            }
+                            *wrklog << std::endl;
+                        }
+
                         // To reduce memory requirements we split long time
                         // intervals into several slices.
                         GTime tstart = tmin;
@@ -795,6 +828,10 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
                                 *wrklog << " - ";
                                 *wrklog << tstop.convert(m_cta_ref);
                                 *wrklog << " s" << std::endl;
+                                if (model->name().length() > 0) {
+                                    *wrklog << " [" << model->name() << "]";
+                                }
+                                *wrklog << std::endl;
                             }
                             // Get photons
@@ -66,6 +66,8 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
+    const double&        max_rate(void) const;
+    void                 max_rate(const double& max_rate);
 protected:
     // Protected methods
@@ -95,6 +97,7 @@ protected:
     mutable bool          m_save_and_dispose; //!< Save and dispose immediately
     double                m_area;             //!< Surface area for simulation (cm2)
     int                   m_max_photons;      //!< Maximum number of photons/slice
+    double                m_max_rate;         //!< Maximum photon rate
     std::vector<GRan>     m_rans;             //!< Random number generators
     int                   m_event_id;         //!< Event identifier
 };
@@ -108,7 +111,32 @@ protected:
 inline
 const GObservations& ctobssim::obs(void) const
 {
-    return m_obs;
+    return (m_obs);
+}
+
+
+/***********************************************************************//**
+ * @brief Return maximum photon rate
+ *
+ * @return Reference to maximum photon rate.
+ ***************************************************************************/
+inline
+const double& ctobssim::max_rate(void) const
+{
+    return (m_max_rate);
+}
+
+
+/***********************************************************************//**
+ * @brief Set maximum photon rate
+ *
+ * @param[in] max_rate Maximum photon rate.
+ ***************************************************************************/
+inline
+void ctobssim::max_rate(const double& max_rate)
+{
+    m_max_rate = max_rate;
+    return;
 }
 #endif /* CTOBSSIM_HPP */
@@ -29,16 +29,17 @@ irf,   s, a, &quot;cta_dummy_irf&quot;,,, &quot;Instrument response function&quot;
 #
 # Simulation parameters
 #======================
-seed,   i, h, 1,,, "Random number generator seed"
-ra,     r, a, 83.63,0,360, "RA of pointing (degrees)"
-dec,    r, a, 22.01,-90,90, "Dec of pointing (degrees)"
-rad,    r, a, 5.0,0,180, "Radius of FOV (degrees)"
-tmin,   r, a, 0.0,,, "Start time (MET in s)"
-tmax,   r, a, 1800.0,,, "End time (MET in s)"
-emin,   r, a, 0.1,,, "Lower energy limit (TeV)"
-emax,   r, a, 100.0,,, "Upper energy limit (TeV)"
-edisp,  b, h, no,,, "Apply energy dispersion?"
-deadc,  r, h, 0.95,0,1, "Average deadtime correction factor"
+seed,    i, h, 1,,, "Random number generator seed"
+ra,      r, a, 83.63,0,360, "RA of pointing (degrees)"
+dec,     r, a, 22.01,-90,90, "Dec of pointing (degrees)"
+rad,     r, a, 5.0,0,180, "Radius of FOV (degrees)"
+tmin,    r, a, 0.0,,, "Start time (MET in s)"
+tmax,    r, a, 1800.0,,, "End time (MET in s)"
+emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
+emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
+edisp,   b, h, no,,, "Apply energy dispersion?"
+deadc,   r, h, 0.95,0,1, "Average deadtime correction factor"
+maxrate, r, h, 1.0e6,,, "Maximum photon rate"
 #
 # Standard parameters
 #====================