From 4104319bea27dc3c08b07e46672e23a1b74018a6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 5 Oct 2018 16:30:34 +0200 Subject: [PATCH] Complete implementation of general CTA background model (#2694) --- inst/cta/include/GCTAModelRadial.hpp | 38 ++++++++------------------------------ inst/cta/pyext/GCTAModelRadial.i | 19 ++++++++++++++----- inst/cta/src/GCTAModelBackground.cpp | 89 ++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------- inst/cta/src/GCTAModelRadial.cpp | 139 +++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------ 4 files changed, 101 insertions(+), 184 deletions(-) diff --git a/inst/cta/include/GCTAModelRadial.hpp b/inst/cta/include/GCTAModelRadial.hpp index 4871fc2..efd5be9 100644 --- a/inst/cta/include/GCTAModelRadial.hpp +++ b/inst/cta/include/GCTAModelRadial.hpp @@ -56,19 +56,15 @@ public: // Operators virtual GCTAModelRadial& operator=(const GCTAModelRadial& model); - virtual GModelPar& operator[](const int& index); - virtual const GModelPar& operator[](const int& index) const; - virtual GModelPar& operator[](const std::string& name); - virtual const GModelPar& operator[](const std::string& name) const; // Implemented virtual methods - virtual double eval(const GCTAInstDir& dir, - const GEnergy& energy, - const GTime& time, - const bool& gradients = false) const; - virtual GCTAInstDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const; + virtual double eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GCTAInstDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; // Pure virtual methods virtual void clear(void) = 0; @@ -83,31 +79,13 @@ public: // Derived class pure virtual methods virtual double eval(const double& offset, const bool& gradients = false) const = 0; - virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0; - - // Methods - int size(void) const; + virtual GCTAInstDir mc(GRan& ran) const = 0; protected: // Protected methods void init_members(void); void copy_members(const GCTAModelRadial& model); void free_members(void); - - // Proteced members - std::vector m_pars; //!< Parameter pointers }; - -/***********************************************************************//** - * @brief Return number of model parameters - * - * @return Number of model parameters. - ***************************************************************************/ -inline -int GCTAModelRadial::size(void) const -{ - return ((int)m_pars.size()); -} - #endif /* GCTAMODELRADIAL_HPP */ diff --git a/inst/cta/pyext/GCTAModelRadial.i b/inst/cta/pyext/GCTAModelRadial.i index e83b139..ea6773c 100644 --- a/inst/cta/pyext/GCTAModelRadial.i +++ b/inst/cta/pyext/GCTAModelRadial.i @@ -42,20 +42,29 @@ public: GCTAModelRadial(const GCTAModelRadial& model); virtual ~GCTAModelRadial(void); + // Implemented virtual methods + virtual double eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GCTAInstDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + + // Pure virtual methods virtual void clear(void) = 0; virtual GCTAModelRadial* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const double& offset, - const bool& gradients = false) const = 0; - virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0; virtual double omega(void) const = 0; virtual void read(const GXmlElement& xml) = 0; virtual void write(GXmlElement& xml) const = 0; - // Methods - int size(void) const; + // Derived class pure virtual methods + virtual double eval(const double& offset, + const bool& gradients = false) const = 0; + virtual GCTAInstDir mc(GRan& ran) const = 0; }; diff --git a/inst/cta/src/GCTAModelBackground.cpp b/inst/cta/src/GCTAModelBackground.cpp index 38c7557..75f94ec 100644 --- a/inst/cta/src/GCTAModelBackground.cpp +++ b/inst/cta/src/GCTAModelBackground.cpp @@ -551,18 +551,26 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation& obs, // TODO: Handle special case of cube spatial model. This will // replace the spectral model by a spectral nodes model. + // Get solid angle of spatial model. This only works for an energy + // independent spatial model!!! + double solidangle = m_spatial->npred(GEnergy(), GTime(), obs); + // Loop over all energy boundaries for (int ieng = 0; ieng < ebounds.size(); ++ieng) { // Compute the background rate in model within the energy - // boundaries from spectral component (units: cts/s). - // Note that the time here is ontime. Deadtime correction will - // be done later. - double rate = spectral->flux(ebounds.emin(ieng), ebounds.emax(ieng)); + // boundaries from spectral component (units: cts/s/sr). + double flux = spectral->flux(ebounds.emin(ieng), ebounds.emax(ieng)); + + // Derive expecting rate (units: cts/s). Note that the time here + // is good time. Deadtime correction will be done later. + double rate = flux * solidangle; // Debug option: dump rate #if defined(G_DUMP_MC) std::cout << "GCTAModelBackground::mc(\"" << name() << "\": "; + std::cout << "flux=" << flux << " cts/s/sr, "; + std::cout << "solidangle=" << solidangle << " sr, "; std::cout << "rate=" << rate << " cts/s)" << std::endl; #endif @@ -593,8 +601,8 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation& obs, #if defined(G_DUMP_MC) std::cout << " Interval " << itime; std::cout << " events=" << n_events << std::endl; - int n_killed_by_deadtime = 0; - int n_killed_by_roi = 0; + int n_removed_by_deadtime = 0; + int n_trial_outside_roi = 0; #endif // Loop over events @@ -605,7 +613,7 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation& obs, if (deadc < 1.0) { if (ran.uniform() > deadc) { #if defined(G_DUMP_MC) - n_killed_by_deadtime++; + n_removed_by_deadtime++; #endif continue; } @@ -617,42 +625,53 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation& obs, times[i], ran); - // Get Monte Carlo event direction from spatial model. - // This only will set the DETX and DETY coordinates. - GCTAInstDir instdir = m_spatial->mc(energy, times[i], ran); + // Get an instrument direction within the RoI. This is + // potentially tried 100 times so that if we really can't + // a valid instrument direction the code is not locked up + for (int k = 0; k < 100; ++k) { - // Derive sky direction from instrument coordinates - GSkyDir skydir = pnt.skydir(instdir); + // Get Monte Carlo event direction from spatial model. + // This only will set the DETX and DETY coordinates. + GCTAInstDir instdir = m_spatial->mc(energy, times[i], ran); - // Set sky direction in GCTAInstDir object - instdir.dir(skydir); + // Derive sky direction from instrument coordinates + GSkyDir skydir = pnt.skydir(instdir); - // Allocate event - GCTAEventAtom event; + // Set sky direction in GCTAInstDir object + instdir.dir(skydir); - // Set event attributes - event.dir(instdir); - event.energy(energy); - event.time(times[i]); + // Allocate event + GCTAEventAtom event; - // Append event to list if it falls in ROI - if (events.roi().contains(event)) { - list->append(event); - } - #if defined(G_DUMP_MC) - else { - n_killed_by_roi++; - } - #endif + // Set event attributes + event.dir(instdir); + event.energy(energy); + event.time(times[i]); + + // Append event to list if it falls in RoI and break + // the look + if (events.roi().contains(event)) { + list->append(event); + break; + } + + // ... otherwise optionally bookkeep the trial + #if defined(G_DUMP_MC) + else { + n_trial_outside_roi++; + } + #endif + + } // endfor: trial look for instrument direction within RoI } // endfor: looped over all events - // Debug option: provide statisics + // Debug option: provide statisics #if defined(G_DUMP_MC) - std::cout << " Killed by deadtime="; - std::cout << n_killed_by_deadtime << std::endl; - std::cout << " Killed by ROI="; - std::cout << n_killed_by_roi << std::endl; + std::cout << " Events removed due to deadtime="; + std::cout << n_removed_by_deadtime << std::endl; + std::cout << " Trials outside RoI="; + std::cout << n_trial_outside_roi << std::endl; #endif } // endfor: looped over all GTIs @@ -987,7 +1006,7 @@ void GCTAModelBackground::set_pointers(void) * @param[in] spatial XML element. * @return Pointer to spatial model. * - * Returns pointer to a radial model that is defined in an XML element. + * Returns pointer to a spatial model that is defined in an XML element. ***************************************************************************/ GCTAModelSpatial* GCTAModelBackground::xml_spatial(const GXmlElement& spatial) const { diff --git a/inst/cta/src/GCTAModelRadial.cpp b/inst/cta/src/GCTAModelRadial.cpp index f428690..ed028c5 100644 --- a/inst/cta/src/GCTAModelRadial.cpp +++ b/inst/cta/src/GCTAModelRadial.cpp @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GCTAModelRadial.cpp - * @brief Abstract radial acceptance model class implementation + * @brief Abstract radial background model class implementation * @author Juergen Knoedlseder */ @@ -66,7 +66,7 @@ GCTAModelRadial::GCTAModelRadial(void) /***********************************************************************//** * @brief Copy constructor * - * @param[in] model Radial acceptance model. + * @param[in] model Radial background model. ***************************************************************************/ GCTAModelRadial::GCTAModelRadial(const GCTAModelRadial& model) { @@ -103,9 +103,9 @@ GCTAModelRadial::~GCTAModelRadial(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] model Radial acceptance model. + * @param[in] model Radial background model. ***************************************************************************/ -GCTAModelRadial& GCTAModelRadial::operator= (const GCTAModelRadial& model) +GCTAModelRadial& GCTAModelRadial::operator=(const GCTAModelRadial& model) { // Execute only if object is not identical if (this != &model) { @@ -126,100 +126,6 @@ GCTAModelRadial& GCTAModelRadial::operator= (const GCTAModelRadial& model) } -/***********************************************************************//** - * @brief Returns model parameter - * - * @param[in] index Parameter index [0,...,size()-1]. - * - * @exception GException::out_of_range - * Parameter index is out of range. - ***************************************************************************/ -GModelPar& GCTAModelRadial::operator[](const int& index) -{ - // Compile option: raise exception if index is out of range - #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) - throw GException::out_of_range(G_ACCESS1, index, 0, size()-1); - #endif - - // Return reference - return *(m_pars[index]); -} - - -/***********************************************************************//** - * @brief Returns model parameter (const version) - * - * @param[in] index Parameter index [0,...,size()-1]. - * - * @exception GException::out_of_range - * Parameter index is out of range. - ***************************************************************************/ -const GModelPar& GCTAModelRadial::operator[](const int& index) const -{ - // Compile option: raise exception if index is out of range - #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) - throw GException::out_of_range(G_ACCESS1, index, 0, size()-1); - #endif - - // Return reference - return *(m_pars[index]); -} - - -/***********************************************************************//** - * @brief Returns reference to model parameter - * - * @param[in] name Parameter name. - * - * @exception GException::par_not_found - * Parameter with specified name not found in container. - ***************************************************************************/ -GModelPar& GCTAModelRadial::operator[](const std::string& name) -{ - // Get parameter index - int index = 0; - for (; index < size(); ++index) { - if (m_pars[index]->name() == name) - break; - } - - // Throw exception if parameter name was not found - if (index >= size()) - throw GException::par_not_found(G_ACCESS2, name); - - // Return reference - return *(m_pars[index]); -} - - -/***********************************************************************//** - * @brief Returns reference to model parameter (const version) - * - * @param[in] name Parameter name. - * - * @exception GException::par_not_found - * Parameter with specified name not found in container. - ***************************************************************************/ -const GModelPar& GCTAModelRadial::operator[](const std::string& name) const -{ - // Get parameter index - int index = 0; - for (; index < size(); ++index) { - if (m_pars[index]->name() == name) - break; - } - - // Throw exception if parameter name was not found - if (index >= size()) - throw GException::par_not_found(G_ACCESS2, name); - - // Return reference - return *(m_pars[index]); -} - - /*========================================================================== = = = Public methods = @@ -230,39 +136,50 @@ const GModelPar& GCTAModelRadial::operator[](const std::string& name) const * @brief Evaluate function * * @param[in] dir Event direction. - * @param[in] energy Event energy. - * @param[in] time Event time. + * @param[in] energy Event energy (not used). + * @param[in] time Event time (not used). * @param[in] gradients Compute gradients? * @return Function value * - * @todo Implement method + * Evaluate radial model for a given event direction. The energy and time of + * the event are not used. ***************************************************************************/ double GCTAModelRadial::eval(const GCTAInstDir& dir, const GEnergy& energy, const GTime& time, const bool& gradients) const { - // Return - return 0.0; + // Compute offset angle in degrees + double offset = dir.theta() * gammalib::rad2deg; + + // Evaluate function + double value = eval(offset, gradients); + + // Return value + return value; } /***********************************************************************//** * @brief Returns MC instrument direction * - * @param[in] energy Event energy. - * @param[in] time Event time. + * @param[in] energy Event energy (not used). + * @param[in] time Event time (not used). * @param[in,out] ran Random number generator. * @return Instrument direction * - * @todo Implement method + * Return random instrument direction. The energy and time of the event are + * not used. ***************************************************************************/ GCTAInstDir GCTAModelRadial::mc(const GEnergy& energy, const GTime& time, GRan& ran) const { - // Return - return GCTAInstDir(); + // Get random instrument direction + GCTAInstDir dir = mc(ran); + + // Return instrument direction + return dir; } @@ -277,9 +194,6 @@ GCTAInstDir GCTAModelRadial::mc(const GEnergy& energy, ***************************************************************************/ void GCTAModelRadial::init_members(void) { - // Initialise members - m_pars.clear(); - // Return return; } @@ -292,9 +206,6 @@ void GCTAModelRadial::init_members(void) ***************************************************************************/ void GCTAModelRadial::copy_members(const GCTAModelRadial& model) { - // Copy members - m_pars = model.m_pars; - // Return return; } -- libgit2 0.21.4