diff --git a/scripts/cslightcrv.py b/scripts/cslightcrv.py index 05d1db6..6252230 100755 --- a/scripts/cslightcrv.py +++ b/scripts/cslightcrv.py @@ -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")) diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 7070f2f..5b875c0 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -244,26 +244,40 @@ void ctbin::run(void) // Get CTA observation GCTAObservation* obs = dynamic_cast(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(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 { diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index 731910f..a6efa81 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -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 diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp index 8f8cb15..90f610e 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -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) { diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index b1a9f08..862c0fa 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -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. diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index e3d559e..9c37587 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -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; diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 17aed86..075fbdd 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -327,82 +327,89 @@ void ctobssim::run(void) // Get pointer on CTA observation GCTAObservation* obs = dynamic_cast(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 diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp index c11bdf8..bbb083a 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -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) { diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 200cc20..6e310dc 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -249,74 +249,88 @@ void ctselect::run(void) // Get CTA observation GCTAObservation* obs = dynamic_cast(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 diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 69bd7d2..c5b9b6d 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -250,43 +250,43 @@ void ctskymap::run(void) // Get CTA observation GCTAObservation* obs = dynamic_cast(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; }