Exclude consistently binned and non CTA observations from the processing.

Jürgen Knödlseder
1 parent 51c98949
Showing 5 changed files with 188 additions and 198 deletions
scripts/cslightcrv.py
src/ctbin/ctbin.cpp
src/ctobssim/ctobssim.cpp
src/ctselect/ctselect.cpp
src/ctskymap/ctskymap.cpp
@@ -112,7 +112,6 @@ class cslightcrv(ctools.cscript):
             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"))
             pars.append(gammalib.GApplicationPar("tbinfile","f","a","tbins.fits","","", "File defining the time binning"))
-            #pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each time bin"))
             pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins per light curve bin (0=unbinned)"))
             pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit of events (TeV)"))
             pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit of events (TeV)"))
@@ -244,37 +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) {
-
-            // Write header for observation
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
             }
+            continue;
+        }
-            // 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;
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
+            if (logTerse()) {
+                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
@@ -327,93 +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) {
-
-                // Write header for observation
+            // Skip observation if it's not CTA
+            if (obs == NULL) {
                 if (logTerse()) {
-                    if (obs->name().length() > 1) {
-                        wrklog.header3("Observation "+obs->name());
-                    }
-                    else {
-                        wrklog.header3("Observation");
-                    }
+                    log << "Warning: Skipping "+m_obs[i]->instrument();
+                    log << " observation \"";
+                    log << m_obs[i]->name() << "\"" << std::endl;
                 }
+                continue;
+            }
-                // Skip observation if we don't have an unbinned observation
-                if (obs->eventtype() != "EventList") {
+            // Skip observation if we have a binned observation
+            if (obs->eventtype() == "CountsCube") {
+                if (logTerse()) {
+                    log << "Warning: Skipping binned observation \"";
+                    log << obs->name()+"\"" << std::endl;
+                }
+                continue;
+            }
-                    // Log that we skip the this observation
-                    if (logTerse()) {
-                        log << "Warning: Skipping binned observation \""+obs->name()+"\"";
-                        log << std::endl;
-                    }
-                    continue;
+            // Write header for observation
+            if (logTerse()) {
+                if (obs->name().length() > 1) {
+                    wrklog.header3("Observation "+obs->name());
+                }
+                else {
+                    wrklog.header3("Observation");
                 }
+            }
-                // 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;
+            // 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;
-                // Save number of events before entering simulation
-                int events_before = obs_clone.events()->size();
+            // Save number of events before entering simulation
+            int events_before = obs_clone.events()->size();
-                // Simulate source events
-                simulate_source(&obs_clone, models, m_rans[i], &wrklog);
+            // Simulate source events
+            simulate_source(&obs_clone, models, m_rans[i], &wrklog);
-                // Simulate source events
-                simulate_background(&obs_clone, models, m_rans[i], &wrklog);
+            // Simulate source events
+            simulate_background(&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;
-                }
-
-                // 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();
-                    }
-
-                    // Store output file name in original observation
-                    obs->eventfile(outfile);
+            // If requested, event lists are saved immediately
+            if (m_save_and_dispose) {
-                    // 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
@@ -249,85 +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) {
-
-            // Write header for observation
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
             }
+            continue;
+        }
-            // Increment counter
-            n_observations++;
-
-            // 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;
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
+            if (logTerse()) {
+                log << "Warning: Skipping binned observation \"";
+                log << obs->name()+"\"" << std::endl;
             }
+            continue;
+        }
-            // 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;
+        // Write header for observation
+        if (logTerse()) {
+            if (obs->name().length() > 1) {
+                log.header3("Observation "+obs->name());
+            }
+            else {
+                log.header3("Observation");
             }
+        }
+
+        // Increment counter
+        n_observations++;
+
+        // 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
@@ -250,54 +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) {
-
-            // Write header for current observation
+        // Skip observation if it's not CTA
+        if (obs == NULL) {
             if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation");
-                }
+                log << "Warning: Skipping "+m_obs[i]->instrument();
+                log << " observation \"";
+                log << m_obs[i]->name() << "\"" << std::endl;
             }
+            continue;
+        }
-            // Skip observation if we don't have an unbinned observation
-            if (obs->eventtype() != "EventList") {
+        // Skip observation if we have a binned observation
+        if (obs->eventtype() == "CountsCube") {
+            if (logTerse()) {
+                log << "Warning: Skipping binned observation \"";
+                log << obs->name()+"\"" << std::endl;
+            }
+            continue;
+        }
-                // Log that we skip the this observation
-                if (logTerse()) {
-                    log << "Warning: Skipping binned observation \""+obs->name()+"\"";
-                    log << std::endl;
-                }
-                continue;
+        // Write header for current observation
+        if (logTerse()) {
+            if (obs->name().length() > 1) {
+                log.header3("Observation "+obs->name());
+            }
+            else {
+                log.header3("Observation");
             }
+        }
-            // Map events into sky map
-            map_events(obs);
+        // Map events into sky map
+        map_events(obs);
-            // Increment observation counter
-            num_obs++;
-
-        } // endif: CTA observation found
+        // 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;
 }