diff --git a/doc/source/reference_manual/ctobssim.rst b/doc/source/reference_manual/ctobssim.rst index 8ed50fb..27754f5 100644 --- a/doc/source/reference_manual/ctobssim.rst +++ b/doc/source/reference_manual/ctobssim.rst @@ -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 ------------------- diff --git a/pyext/ctobssim.i b/pyext/ctobssim.i index eb3149d..d44fb9f 100644 --- a/pyext/ctobssim.i +++ b/pyext/ctobssim.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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); }; diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 5cf4e33..a1c9390 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -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& 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& 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& 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 diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp index 0abae1a..8476ca7 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -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 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 */ diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index 6b8e59f..548e595 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -29,16 +29,17 @@ irf, s, a, "cta_dummy_irf",,, "Instrument response function" # # 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 #====================