Adapted ctools to handle observation type consistently

Michael Mayer
1 parent 23dd86c9
Showing 10 changed files with 91 additions and 149 deletions
scripts/cslightcrv.py
src/ctbin/ctbin.cpp
src/ctbkgcube/ctbkgcube.cpp
src/ctexpcube/ctexpcube.cpp
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
src/ctobssim/ctobssim.cpp
src/ctpsfcube/ctpsfcube.cpp
src/ctselect/ctselect.cpp
src/ctskymap/ctskymap.cpp
@@ -107,7 +107,7 @@ class cslightcrv(ctools.cscript):
             pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
             pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
-            pars.append(gammalib.GApplicationPar("tbinalg","s","h","GTI","FILE|LIN|GTI","", "Algorithm for defining time bins"))
+            pars.append(gammalib.GApplicationPar("tbinalg","s","a","GTI","FILE|LIN|GTI","", "Algorithm for defining time bins"))
             pars.append(gammalib.GApplicationPar("tmin","r","a","51544.5","","days","Lightcurve start time [MJD]"))
             pars.append(gammalib.GApplicationPar("tmax","r","a","51544.6","","days","Lightcurve stop time [MJD]"))
             pars.append(gammalib.GApplicationPar("tbins","i","a","5","","","Number of time bins"))
@@ -257,6 +257,17 @@ void ctbin::run(void)
                 }
             }
+            // Skip observation if we don't have an unbinned observation
+            if (obs->eventtype() != "EventList") {
+
+                // Log that we skip the this observation
+                if (logTerse()) {
+                    log << "Warning: Skipping binned observation \""+obs->name()+"\"";
+                    log << std::endl;
+                }
+                continue;
+            }
+
             // Fill the cube
             fill_cube(obs);
@@ -561,11 +572,25 @@ void ctbin::obs_cube(void)
         // Attach event cube to CTA observation
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
         if (obs != NULL) {
-            obs->events(this->cube());
-            obs->eventfile("");
-        }
-    }
+            // Only change the event type if we had an unbinned observation
+            if (obs->eventtype() == "EventList") {
+
+                obs->events(this->cube());
+                obs->eventfile("");
+             } // endif: observation was unbinned
+            else { // input observation was binned and skipped
+
+                // Create new and empty cube
+                GCTAEventCube cube = GCTAEventCube(m_cube, m_ebounds, obs->gti());
+
+                // Assign empty cube to have a new binned observation
+                obs->events(cube);
+            }
+
+        } // endif: obervation was valid
+
+    } // endif: we only had one observation in the container
     // ... otherwise put a single CTA observation in container
     else {
@@ -308,7 +308,7 @@ void ctbkgcube::run(void)
     m_obs.models(m_bkgmdl);
     // Fill background cube from observations
-    m_background.fill(m_obs);
+    m_background.fill(m_obs, &log);
     // Create a background model for the output background cube and append
     // that model to the input model in place of the original
@@ -248,7 +248,7 @@ void ctexpcube::run(void)
     }
     // Fill exposure
-    m_expcube.fill(m_obs);
+    m_expcube.fill(m_obs, &log);
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -427,7 +427,7 @@ void ctmodel::get_parameters(void)
     // If there are no observations in container then load them via user
     // parameters
     if (m_obs.size() == 0) {
-        get_obs();
+        m_obs = get_observations();
     }
     // Read model definition file if required
@@ -505,126 +505,6 @@ void ctmodel::get_parameters(void)
 /***********************************************************************//**
- * @brief Get observation container
- *
- * Get an observation container according to the user parameters. The method
- * supports loading of a individual FITS file or an observation definition
- * file in XML format. If the input filename is empty, parameters are read
- * to build a CTA observation from scratch.
- ***************************************************************************/
-void ctmodel::get_obs(void)
-{
-    // Get the filename from the input parameters
-    std::string filename = (*this)["inobs"].filename();
-
-    // If no observation definition file has been specified then read all
-    // parameters that are necessary to create an observation from scratch
-    if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) {
-
-        // Get response cube filenames
-        std::string expcube = (*this)["expcube"].filename();
-        std::string psfcube = (*this)["psfcube"].filename();
-        std::string bkgcube = (*this)["bkgcube"].filename();
-
-        // If the filenames are valid then build an observation from cube
-        // response information
-        if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") &&
-            (gammalib::strip_whitespace(expcube) != "") &&
-            (gammalib::strip_whitespace(psfcube) != "") &&
-            (gammalib::strip_whitespace(bkgcube) != "")) {
-
-            // Get exposure, PSF and background cubes
-            GCTACubeExposure   exposure(expcube);
-            GCTACubePsf        psf(psfcube);
-            GCTACubeBackground background(bkgcube);
-
-            // Create energy boundaries
-            GEbounds ebounds = create_ebounds();
-
-            // Create dummy sky map cube
-            GSkymap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size());
-            m_append_cube = true;
-
-            // Create event cube
-            GCTAEventCube cube(map, ebounds, exposure.gti());
-
-            // Create CTA observation
-            GCTAObservation cta;
-            cta.events(cube);
-            cta.response(exposure, psf, background);
-
-            // Append observation to container
-            m_obs.append(cta);
-
-        } // endif: cube response information was available
-
-        // ... otherwise build an observation from IRF response information
-        else {
-
-            // Create CTA observation
-            GCTAObservation cta = create_cta_obs();
-
-            // Set response
-            set_obs_response(&cta);
-
-            // Append observation to container
-            m_obs.append(cta);
-            
-        }
-
-    } // endif: filename was "NONE" or ""
-
-    // ... otherwise we have a file name
-    else {
-
-        // If file is a FITS file then create an empty CTA observation
-        // and load file into observation
-        if (gammalib::is_fits(filename)) {
-
-            // Allocate empty CTA observation
-            GCTAObservation cta;
-
-            // Load data
-            cta.load(filename);
-
-            // Set response
-            set_obs_response(&cta);
-
-            // Append observation to container
-            m_obs.append(cta);
-
-            // Signal that no XML file should be used for storage
-            m_use_xml = false;
-
-        }
-
-        // ... otherwise load file into observation container
-        else {
-
-            // Load observations from XML file
-            m_obs.load(filename);
-
-            // For all observations that have no response, set the response
-            // from the task parameters
-            set_response(m_obs);
-
-            // Set observation boundary parameters (emin, emax, rad)
-            set_obs_bounds(m_obs);
-
-            // Signal that XML file should be used for storage
-            m_use_xml = true;
-
-        } // endelse: file was an XML file
-
-    }
-
-    // Return
-    return;
-
-}
-
-
-/***********************************************************************//**
  * @brief Fill model into model cube
  *
  * @param[in] obs CTA observation.
@@ -79,7 +79,6 @@ protected:
     void copy_members(const ctmodel& app);
     void free_members(void);
     void get_parameters(void);
-    void get_obs(void);
     void fill_cube(const GCTAObservation* obs);
     bool has_cube(void) const;
@@ -340,6 +340,17 @@ void ctobssim::run(void)
                     }
                 }
+                // Skip observation if we don't have an unbinned observation
+                if (obs->eventtype() != "EventList") {
+
+                    // Log that we skip the this observation
+                    if (logTerse()) {
+                        log << "Warning: Skipping binned observation \""+obs->name()+"\"";
+                        log << std::endl;
+                    }
+                    continue;
+                }
+
                 // Work on a clone of the CTA observation. This makes sure that
                 // any memory allocated for computing (for example a response
                 // cache) is properly de-allocated on exit of this run
@@ -247,7 +247,7 @@ void ctpsfcube::run(void)
     }
     // Fill PSF 
-    m_psfcube.fill(m_obs);
+    m_psfcube.fill(m_obs, &log);
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -268,6 +268,17 @@ void ctselect::run(void)
             // Save event file name (for possible saving)
             m_infiles[i] = obs->eventfile();
+            // Skip observation if we don't have an unbinned observation
+            if (obs->eventtype() != "EventList") {
+
+                // Log that we skip the this observation
+                if (logTerse()) {
+                    log << "Warning: Skipping binned observation \""+obs->name()+"\"";
+                    log << std::endl;
+                }
+                continue;
+            }
+
             // Fall through in case that the event file is empty
             if (obs->events()->size() == 0) {
                 if (logTerse()) {
@@ -1218,28 +1229,33 @@ void ctselect::save_event_list(const GCTAObservation* obs,
     // Save only if observation is valid
     if (obs != NULL) {
-        // Save observation into FITS file
-        obs->save(outfile, clobber());
-
-        // Copy all extensions other than EVENTS and GTI from the input to
-        // the output event list. The EVENTS and GTI extensions are written
-        // by the save method, all others that may eventually be present
-        // have to be copied by hand.
-        GFits infits(infile);
-        GFits outfits(outfile);
-        for (int extno = 1; extno < infits.size(); ++extno) {
-            GFitsHDU* hdu = infits.at(extno);
-            if (hdu->extname() != "EVENTS" && hdu->extname() != "GTI") {
-                outfits.append(*hdu);
+        // Save only if we have an event list
+        if (obs->eventtype() == "EventList") {
+
+            // Save observation into FITS file
+            obs->save(outfile, clobber());
+
+            // Copy all extensions other than EVENTS and GTI from the input to
+            // the output event list. The EVENTS and GTI extensions are written
+            // by the save method, all others that may eventually be present
+            // have to be copied by hand.
+            GFits infits(infile);
+            GFits outfits(outfile);
+            for (int extno = 1; extno < infits.size(); ++extno) {
+                GFitsHDU* hdu = infits.at(extno);
+                if (hdu->extname() != "EVENTS" && hdu->extname() != "GTI") {
+                    outfits.append(*hdu);
+                }
             }
-        }
-        // Close input file
-        infits.close();
+            // Close input file
+            infits.close();
+
+            // Save file to disk and close it (we need both operations)
+            outfits.save(true);
+            outfits.close();
-        // Save file to disk and close it (we need both operations)
-        outfits.save(true);
-        outfits.close();
+        } // endif: observation was unbinned
     } // endif: observation was valid
@@ -263,6 +263,17 @@ void ctskymap::run(void)
                 }
             }
+            // Skip observation if we don't have an unbinned observation
+            if (obs->eventtype() != "EventList") {
+
+                // Log that we skip the this observation
+                if (logTerse()) {
+                    log << "Warning: Skipping binned observation \""+obs->name()+"\"";
+                    log << std::endl;
+                }
+                continue;
+            }
+
             // Map events into sky map
             map_events(obs);