Add support for general GModelData simulations.

Background simulations can now be done for any model of type GModelData. This is thanks to the abstract mc() method that has been introduced in GModelData. A test case for the new GCTABackground model has been added to the unit tests.

Add support for general GModelData simulations.
Background simulations can now be done for any model of type GModelData. This is thanks to the abstract mc() method that has been introduced in GModelData. A test case for the new GCTABackground model has been added to the unit tests.
Jürgen Knödlseder
1 parent 4eab7742
Showing 6 changed files with 82 additions and 34 deletions
ChangeLog
NEWS
src/ctobssim/ctobssim.cpp
test/data/model_background_cube.xml
test/data/test_cube.fits
test/test_ctools.sh
-2013-04-09  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2013-11-19  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
         * Version 0.7.0 released
         ========================
  
         Adapt to GammaLib-00-08-00 interface
-        ctobssim now only simulates events within ROI
+        ctobssim now only simulates events within ROI and supports
+          all kinds of data space models
  
 2012-12-21  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-07-00
  
 Author(s) : Juergen Knoedlseder
  
-April 9th, 2013
+19 November 2013
  
  
 1. Introduction
@@ -33,6 +33,10 @@ ra and dec) and region radius (parameter rad). A margin of 0.5 has been
 added internally to the radius so that also sources at the edge but
 slightly outside of the ROI will be simulated correctly.
  
+ctobssim now supports any kind of data space model for background
+simulations, using the abstract mc() method that is now available
+in the GModelData base class.
+
  
 5. ctselect - CTA event selection
 ---------------------------------
@@ -872,49 +872,56 @@ void ctobssim::simulate_background(GCTAObservation* obs,
         // Loop over all models
         for (int i = 0; i < models.size(); ++i) {
  
-            // Get model (NULL if not a radial acceptance model)
-            const GCTAModelRadialAcceptance* model = 
-                  dynamic_cast<const GCTAModelRadialAcceptance*>(models[i]);
+            // Get data model (NULL if not a data model)
+            const GModelData* model =
+                dynamic_cast<const GModelData*>(models[i]);
  
-            // If we have a radial acceptance model then simulate events
-            if (model != NULL) {
+            // If we have a data model that applies to CTA then simulate
+            // events
+            if (model != NULL && model->isvalid("CTA", "")) {
  
-                // Get simulated event list. Note that this method includes
-                // the deadtime correction.
-                GCTAEventList* list = model->mc(*obs, ran);
+                // Get simulated CTA event list. Note that this method
+                // includes the deadtime correction.
+                GCTAEventList* list =
+                    dynamic_cast<GCTAEventList*>(model->mc(*obs, ran));
  
-                // Reserves space for events
-                events->reserve(list->size()+events->size());
+                // Continue only if we got a CTA event list
+                if (list != NULL) {
  
-                // Append events
-                for (int k = 0; k < list->size(); k++) {
+                    // Reserves space for events
+                    events->reserve(list->size()+events->size());
  
-                    // Get event pointer
-                    GCTAEventAtom* event = (*list)[k];
+                    // Append events
+                    for (int k = 0; k < list->size(); k++) {
  
-                    // Use event only if it falls within ROI
-                    if (event->dir().dist_deg(events->roi().centre().dir()) <=
-                        events->roi().radius()) {
+                        // Get event pointer
+                        GCTAEventAtom* event = (*list)[k];
  
-                        // Set event identifier
-                        event->event_id(m_event_id);
-                        m_event_id++;
+                        // Use event only if it falls within ROI
+                        if (event->dir().dist_deg(events->roi().centre().dir()) <=
+                            events->roi().radius()) {
  
-                        // Append event
-                        events->append(*event);
+                            // Set event identifier
+                            event->event_id(m_event_id);
+                            m_event_id++;
  
-                    } // endif: event was within ROI
+                            // Append event
+                            events->append(*event);
  
-                } // endfor: looped over all events
+                        } // endif: event was within ROI
  
-                // Dump simulation results
-                if (logNormal()) {
-                    *wrklog << gammalib::parformat("MC background events");
-                    *wrklog << list->size() << std::endl;
-                }
+                    } // endfor: looped over all events
+
+                    // Dump simulation results
+                    if (logNormal()) {
+                        *wrklog << gammalib::parformat("MC background events");
+                        *wrklog << list->size() << std::endl;
+                    }
+
+                    // Free event list
+                    delete list;
  
-                // Free event list
-                delete list;
+                } // endif: we had a CTA event list
  
             } // endif: model was valid
  
+<?xml version="1.0" standalone="no"?>	
+<source_library title="source library">
+  <source name="CTABackgroundModel" type="CTABackground" instrument="CTA">
+    <spectrum type="PowerLaw">	
+      <parameter name="Prefactor" scale="1.0"  value="1.0"  min="1e-3" max="1e+3"   free="1"/>	
+      <parameter name="Index"     scale="1.0"  value="0.0"  min="-5.0" max="+5.0"   free="1"/>	
+      <parameter name="Scale"     scale="1e12" value="1.0"  min="0.01" max="1000.0" free="0"/>	
+    </spectrum>
+    <spatialModel type="MapCubeFunction" file="data/test_cube.fits">
+      <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
+    </spatialModel>
+  </source>	
+</source_library>
@@ -61,6 +61,8 @@ export PFILES=pfiles
 # Test ctobssim
 # =============
 $ECHO -n "Test ctobssim: "
+#
+# Run 1
 $ctobssim infile="data/crab.xml" \
           outfile="events.fits" \
           caldb="irf" \
@@ -80,6 +82,27 @@ else
   echo " events.fits file is not valid"
   exit 1
 fi
+#
+# Run 2
+$ctobssim infile="data/model_background_cube.xml" \
+          outfile="events2.fits" \
+          caldb="irf" \
+          irf="cta_dummy_irf" \
+          ra=84.17263 \
+          dec=22.01444 \
+          rad=1.0 \
+          tmin=0.0 \
+          tmax=1800.0 \
+          emin=0.1 \
+          emax=100.0
+$ECHO -n "."
+if [ -s "events2.fits" ]
+then
+  $ECHO -n "."
+else
+  echo " events2.fits file is not valid"
+  exit 1
+fi
 $ECHO " ok"