Merge branch '1439-adapt-ctools-to-consistently-handle-observation-types' into integration

Jürgen Knödlseder
2 parents 40d6cfa9 5c3a8508
Showing 10 changed files with 239 additions and 306 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"))
@@ -244,26 +244,40 @@ void ctbin::run(void)
         // Get CTA observation
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
-        // Continue only if observation is a CTA observation
-        if (obs != NULL) {
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
+            if (logTerse()) {
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
+            }
+            continue;
+        }
  
-            // Write header for observation
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping binned observation \"";
+                log << obs->name()+"\"" << std::endl;
             }
+            continue;
+        }
  
-            // Fill the cube
-            fill_cube(obs);
+        // Write header for observation
+        if (logTerse()) {
+            if (obs->name().length() > 1) {
+                log.header3("Observation "+obs->name());
+            }
+            else {
+                log.header3("Observation");
+            }
+        }
  
-            // Dispose events to free memory
-            obs->dispose_events();
+        // Fill the cube
+        fill_cube(obs);
  
-        } // endif: CTA observation found
+        // Dispose events to free memory
+        obs->dispose_events();
  
     } // endfor: looped over observations
  
@@ -561,11 +575,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;
  
@@ -327,82 +327,89 @@ void ctobssim::run(void)
             // Get pointer on CTA observation
             GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
-            // Continue only if observation is a CTA observation
-            if (obs != NULL) {
+            // Skip observation if it's not CTA
+            if (obs == NULL) {
+                if (logTerse()) {
+                    log << "Warning: Skipping "+m_obs[i]->instrument();
+                    log << " observation \"";
+                    log << m_obs[i]->name() << "\"" << std::endl;
+                }
+                continue;
+            }
  
-                // Write header for observation
+            // Skip observation if we have a binned observation
+            if (obs->eventtype() == "CountsCube") {
                 if (logTerse()) {
-                    if (obs->name().length() > 1) {
-                        wrklog.header3("Observation "+obs->name());
-                    }
-                    else {
-                        wrklog.header3("Observation");
-                    }
+                    log << "Warning: Skipping binned observation \"";
+                    log << obs->name()+"\"" << 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
-                GCTAObservation obs_clone = *obs;
+            // Write header for observation
+            if (logTerse()) {
+                if (obs->name().length() > 1) {
+                    wrklog.header3("Observation "+obs->name());
+                }
+                else {
+                    wrklog.header3("Observation");
+                }
+            }
  
-                // Save number of events before entering simulation
-                int events_before = obs_clone.events()->size();
+            // 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
+            GCTAObservation obs_clone = *obs;
  
-                // Simulate source events
-                simulate_source(&obs_clone, models, m_rans[i], &wrklog);
+            // Save number of events before entering simulation
+            int events_before = obs_clone.events()->size();
  
-                // Simulate source events
-                simulate_background(&obs_clone, models, m_rans[i], &wrklog);
+            // Simulate source events
+            simulate_source(&obs_clone, models, m_rans[i], &wrklog);
  
-                // Dump simulation results
-                if (logNormal()) {
-                    wrklog << gammalib::parformat("MC events");
-                    wrklog << obs_clone.events()->size() - events_before;
-                    wrklog << " (all models)";
-                    wrklog << std::endl;
-                }
+            // Simulate source events
+            simulate_background(&obs_clone, models, m_rans[i], &wrklog);
  
-                // Append the event list to the original observation
-                obs->events(*(obs_clone.events()));
+            // Dump simulation results
+            if (logNormal()) {
+                wrklog << gammalib::parformat("MC events");
+                wrklog << obs_clone.events()->size() - events_before;
+                wrklog << " (all models)";
+                wrklog << std::endl;
+            }
  
-                // If requested, event lists are saved immediately
-                if (m_save_and_dispose) {
+            // Append the event list to the original observation
+            obs->events(*(obs_clone.events()));
  
-                    // Set event output file name. If multiple observations are
-                    // handled, build the filename from prefix and observation
-                    // index. Otherwise use the outfile parameter.
-                    std::string outfile;
-                    if (m_use_xml) {
-                        m_prefix = (*this)["prefix"].string();
-                        outfile  = m_prefix + gammalib::str(i) + ".fits";
-                    }
-                    else {
-                        outfile  = (*this)["outevents"].filename();
-                    }
+            // If requested, event lists are saved immediately
+            if (m_save_and_dispose) {
  
-                    // Store output file name in original observation
-                    obs->eventfile(outfile);
-
-                    // Save observation into FITS file. This is a critical zone
-                    // to avoid multiple threads writing simultaneously
-                    #pragma omp critical
-                    {
-                        obs_clone.save(outfile, clobber());
-                    }
+                // Set event output file name. If multiple observations are
+                // handled, build the filename from prefix and observation
+                // index. Otherwise use the outfile parameter.
+                std::string outfile;
+                if (m_use_xml) {
+                    m_prefix = (*this)["prefix"].string();
+                    outfile  = m_prefix + gammalib::str(i) + ".fits";
+                }
+                else {
+                    outfile  = (*this)["outevents"].filename();
+                }
  
-                    // Dispose events
-                    obs->dispose_events();
+                // Store output file name in original observation
+                obs->eventfile(outfile);
  
+                // Save observation into FITS file. This is a critical zone
+                // to avoid multiple threads writing simultaneously
+                #pragma omp critical
+                {
+                    obs_clone.save(outfile, clobber());
                 }
  
-                // ... otherwise append the event list to the original observation
-                /*
-                else {
-                    obs->events(*(obs_clone.events()));
-                }
-                */
+                // Dispose events
+                obs->dispose_events();
  
-            } // endif: CTA observation found
+            }
  
         } // endfor: looped over observations
  
@@ -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) {
@@ -249,74 +249,88 @@ void ctselect::run(void)
         // Get CTA observation
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
-        // Continue only if observation is a CTA observation
-        if (obs != NULL) {
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
+            if (logTerse()) {
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
+            }
+            continue;
+        }
  
-            // Write header for observation
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping binned observation \"";
+                log << obs->name()+"\"" << std::endl;
             }
+            continue;
+        }
  
-            // Increment counter
-            n_observations++;
+        // Write header for observation
+        if (logTerse()) {
+            if (obs->name().length() > 1) {
+                log.header3("Observation "+obs->name());
+            }
+            else {
+                log.header3("Observation");
+            }
+        }
  
-            // Save event file name (for possible saving)
-            m_infiles[i] = obs->eventfile();
+        // Increment counter
+        n_observations++;
  
-            // Fall through in case that the event file is empty
-            if (obs->events()->size() == 0) {
-                if (logTerse()) {
-                    log << " Warning: No events in event file \"";
-                    log << m_infiles[i] << "\". Event selection skipped.";
-                    log << std::endl;
-                }
-                continue;
-            }
+        // Save event file name (for possible saving)
+        m_infiles[i] = obs->eventfile();
  
-            // Get temporary file name
-            #if G_USE_MKSTEMP
-            char tpl[]  = "ctselectXXXXXX";
-            int  fileid = mkstemp(tpl);
-            std::string filename(tpl);
-            #else
-            std::string filename = std::tmpnam(NULL);
-            #endif
-
-            // Save observation in temporary file
-            obs->save(filename, true);
-
-            // Log saved FITS file
-            if (logExplicit()) {
-                GFits tmpfile(filename);
+        // Fall through in case that the event file is empty
+        if (obs->events()->size() == 0) {
+            if (logTerse()) {
+                log << " Warning: No events in event file \"";
+                log << m_infiles[i] << "\". Event selection skipped.";
                 log << std::endl;
-                log.header1("FITS file content of temporary file");
-                log << tmpfile << std::endl;
-                tmpfile.close();
             }
+            continue;
+        }
  
-            // Check temporary file
-            std::string message = check_infile(filename);
-            if (message.length() > 0) {
-                throw GException::app_error(G_RUN, message);
-            }
+        // Get temporary file name
+        #if G_USE_MKSTEMP
+        char tpl[]  = "ctselectXXXXXX";
+        int  fileid = mkstemp(tpl);
+        std::string filename(tpl);
+        #else
+        std::string filename = std::tmpnam(NULL);
+        #endif
+
+        // Save observation in temporary file
+        obs->save(filename, true);
+
+        // Log saved FITS file
+        if (logExplicit()) {
+            GFits tmpfile(filename);
+            log << std::endl;
+            log.header1("FITS file content of temporary file");
+            log << tmpfile << std::endl;
+            tmpfile.close();
+        }
+
+        // Check temporary file
+        std::string message = check_infile(filename);
+        if (message.length() > 0) {
+            throw GException::app_error(G_RUN, message);
+        }
  
-            // Load observation from temporary file, including event selection
-            select_events(obs, filename);
+        // Load observation from temporary file, including event selection
+        select_events(obs, filename);
  
-            // Close temporary file
-            #if G_USE_MKSTEMP
-            close(fileid);
-            #endif
+        // Close temporary file
+        #if G_USE_MKSTEMP
+        close(fileid);
+        #endif
  
-            // Remove temporary file
-            std::remove(filename.c_str());
-            
-        } // endif: had a CTA observation
+        // Remove temporary file
+        std::remove(filename.c_str());
  
     } // endfor: looped over all observations
  
@@ -1218,28 +1232,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
  
@@ -250,43 +250,43 @@ void ctskymap::run(void)
         // Get CTA observation
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
-        // Continue only if observation is a CTA observation
-        if (obs != NULL) {
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
+            if (logTerse()) {
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
+            }
+            continue;
+        }
  
-            // Write header for current observation
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping binned observation \"";
+                log << obs->name()+"\"" << std::endl;
             }
+            continue;
+        }
  
-            // Map events into sky map
-            map_events(obs);
-            
-            // Increment observation counter
-            num_obs++;
+        // Write header for current observation
+        if (logTerse()) {
+            if (obs->name().length() > 1) {
+                log.header3("Observation "+obs->name());
+            }
+            else {
+                log.header3("Observation");
+            }
+        }
  
-        } // endif: CTA observation found
+        // Map events into sky map
+        map_events(obs);
+            
+        // Increment observation counter
+        num_obs++;
  
     } // endfor: looped over observations
  
-    // Write observation(s) into logger
-/*
-    if (logTerse()) {
-        log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Map observations");
-        }
-        else {
-            log.header1("Map observation");
-        }
-        log << m_obs << std::endl;
-    }
-*/
-
     // Return
     return;
 }