diff --git a/ChangeLog b/ChangeLog index aad047b..565e106 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,7 @@ Clean-up public interface Adapt to GammaLib 0.9.0 interface ctbin now always merges event lists into a single counts cube + ctmodel now always creates a single model cube Add ctexpcube tool (generates exposure cube) Add ctpsfcube tool (generates PSF cube) Add ctbkgcube tool (generates background cube) diff --git a/NEWS b/NEWS index 6601715..9297c6e 100644 --- a/NEWS +++ b/NEWS @@ -98,8 +98,26 @@ extension in the FITS file. None -10. ctmodel - CTA counts map modelling --------------------------------------- +10. ctmodel - CTA counts cube modelling +--------------------------------------- +WARNING: The philosophy behind the ctmodel tool changed fundamentally! + +ctmodel now always produces a single model cube on output. There are +two important input parameters: "infile" and "obsfile". "infile" +specifies a counts or model cube that is used to define the binning +of the model cube. The functionnality for this parameter is as before. +If "infile=NONE", the cube will be defined according to task parameters +that specify the spatial and the energy binning. The "obsfile" +parameter is new and allows specifying the observation definition +for which the model cube should be computed. If an event list or a +counts cube is presented, the pointing and Good Time Intervals will +be extracted from the FITS file header. If an observation definition +XML file is presented, ctmodel will now loop over all CTA observations +in the file and add the model values to the model cube. If no +response information is associated to the "obsfile" (which is typically +the case for event lists or counts cubes), the "caldb" and "irf" +parameters will be used to identify the response. + ctmodel now has full access to the calibration database, allowing the use for 2D response functions in FITS files. diff --git a/pyext/ctmodel.i b/pyext/ctmodel.i index c6c33b7..5bbe12f 100644 --- a/pyext/ctmodel.i +++ b/pyext/ctmodel.i @@ -1,5 +1,5 @@ /*************************************************************************** - * ctmodel - CTA counts model tool * + * ctmodel - CTA model cube generation tool * * ----------------------------------------------------------------------- * * copyright (C) 2012-2014 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctmodel.i - * @brief CTA counts model tool definition + * @brief CTA model cube tool definition * @author Juergen Knoedlseder */ %{ @@ -51,6 +51,7 @@ public: void run(void); void save(void); const GObservations& obs(void) const; + const GCTAEventCube& cube(void) const; }; diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 3c0a33e..cb823db 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctmodel - CTA counts model tool * + * ctmodel - CTA model cube generation tool * * ----------------------------------------------------------------------- * * copyright (C) 2012-2014 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctmodel.cpp - * @brief CTA counts model tool implementation + * @brief CTA model cube tool implementation * @author Juergen Knoedlseder */ @@ -33,8 +33,7 @@ #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_SETUP_OBS "ctmodel::setup_obs()" -#define G_MODEL_MAP "ctmodel::model_map(GCTAObservation*)" +#define G_GET_EBOUNDS "ctmodel::get_ebounds()" /* __ Debug definitions __________________________________________________ */ @@ -212,10 +211,10 @@ void ctmodel::execute(void) // Signal that some parameters should be read ahead m_read_ahead = true; - // Create the model map(s) + // Create model cube run(); - // Save the model map(s) into FITS file + // Save model cube into FITS file save(); // Return @@ -276,23 +275,15 @@ void ctmodel::run(void) // Write header if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Generate model maps"); - } - else { - log.header1("Generate model map"); - } + log.header1("Generate model cube"); } - // Initialise observation counter - int n_observations = 0; + // Initialise counts cube + init_cube(); // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { - // Initialise event input and output filenames - m_infiles.push_back(""); - // Get CTA observation GCTAObservation* obs = dynamic_cast(m_obs[i]); @@ -309,35 +300,18 @@ void ctmodel::run(void) } } - // Increment number of observations - n_observations++; - - // Save event file name (for possible saving) - m_infiles[i] = obs->eventfile(); - - // Generate model map - model_map(obs, m_obs.models()); + // Fill the cube + fill_cube(obs); } // endif: CTA observation found } // endfor: looped over observations - // If more than a single observation has been handled then make sure - // that an XML file will be used for storage - if (n_observations > 1) { - m_use_xml = true; - } - - // Write observation(s) into logger + // Log cube if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations after model map generation"); - } - else { - log.header1("Observation after model map generation"); - } - log << m_obs << std::endl; + log.header1("Model cube"); + log << m_cube << std::endl; } // Restore energy dispersion flag for all CTA observations @@ -354,43 +328,135 @@ void ctmodel::run(void) /***********************************************************************//** - * @brief Save model map(s) - * - * This method saves the model map(s) into FITS file(s). There are two - * modes, depending on the m_use_xml flag. + * @brief Save model cube * - * If m_use_xml is true, all model map(s) will be saved into FITS files, - * where the output filenames are constructued from the input filenames by - * prepending the m_prefix string to name. Any path information will be - * stripped form the input name, hence event files will be written into the - * local working directory (unless some path information is present - * in the prefix). In addition, an XML file will be created that gathers - * the filename information for the model map(s). If an XML file was present - * on input, all metadata information will be copied from this input file. - * - * If m_use_xml is false, the model map will be saved into a FITS file. + * Saves the model cube into a FITS file specified using the "outfile" + * task parameter. ***************************************************************************/ void ctmodel::save(void) { // Write header if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Save observations"); - } - else { - log.header1("Save observation"); - } + log.header1("Save cube"); } - // Case A: Save model map(s) and XML metadata information - if (m_use_xml) { - save_xml(); - } + // Make sure we have the FITS filename + m_outfile = (*this)["outfile"].filename(); - // Case B: Save model map as FITS file - else { - save_fits(); + // Save model cube into FITS file + m_cube.save(m_outfile, clobber()); + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctmodel::init_members(void) +{ + // Initialise members + m_infile.clear(); + m_obsfile.clear(); + m_outfile.clear(); + m_caldb.clear(); + m_irf.clear(); + m_srcmdl.clear(); + m_proj.clear(); + m_coordsys.clear(); + m_ebinalg.clear(); + m_ebinfile.clear(); + m_ra = 0.0; + m_dec = 0.0; + m_deadc = 1.0; + m_tmin = 0.0; + m_tmax = 0.0; + m_emin = 0.0; + m_emax = 0.0; + m_enumbins = 0; + m_xref = 0.0; + m_yref = 0.0; + m_binsz = 0.0; + m_nxpix = 0; + m_nypix = 0; + m_apply_edisp = false; + + // Initialise protected members + m_obs.clear(); + m_cube.clear(); + m_ebounds.clear(); + m_gti.clear(); + m_read_ahead = false; + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctmodel::copy_members(const ctmodel& app) +{ + // Copy attributes + m_infile = app.m_infile; + m_obsfile = app.m_obsfile; + m_outfile = app.m_outfile; + m_caldb = app.m_caldb; + m_irf = app.m_irf; + m_srcmdl = app.m_srcmdl; + m_proj = app.m_proj; + m_coordsys = app.m_coordsys; + m_ebinalg = app.m_ebinalg; + m_ebinfile = app.m_ebinfile; + m_ra = app.m_ra; + m_dec = app.m_dec; + m_deadc = app.m_deadc; + m_tmin = app.m_tmin; + m_tmax = app.m_tmax; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_enumbins = app.m_enumbins; + m_xref = app.m_xref; + m_yref = app.m_yref; + m_binsz = app.m_binsz; + m_nxpix = app.m_nxpix; + m_nypix = app.m_nypix; + m_apply_edisp = app.m_apply_edisp; + + // Copy protected members + m_obs = app.m_obs; + m_cube = app.m_cube; + m_ebounds = app.m_ebounds; + m_gti = app.m_gti; + m_read_ahead = app.m_read_ahead; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctmodel::free_members(void) +{ + // Write separator into logger + if (logTerse()) { + log << std::endl; } // Return @@ -406,62 +472,77 @@ void ctmodel::save(void) ***************************************************************************/ void ctmodel::get_parameters(void) { - // If we do not have any observations in the container then get all - // information to build (at least) one - if (m_obs.size() == 0) { + // Read input and optionally output cube filenames + m_infile = (*this)["infile"].filename(); + if (m_read_ahead) { + m_outfile = (*this)["outfile"].filename(); + } - // Read general parameters. Only read the output parameters if - // required, and also only read the model XML filename if no - // model is yet in the container. - m_infile = (*this)["infile"].filename(); - if (m_read_ahead) { - m_outfile = (*this)["outfile"].filename(); - m_prefix = (*this)["prefix"].string(); - } - m_caldb = (*this)["caldb"].string(); - m_irf = (*this)["irf"].string(); - if (m_obs.models().size() == 0) { - m_srcmdl = (*this)["srcmdl"].filename(); + // Read model definiton file if required + if (m_obs.models().size() == 0) { + m_srcmdl = (*this)["srcmdl"].filename(); + } + + // Get energy dispersion flag parameters + m_apply_edisp = (*this)["edisp"].boolean(); + + // If the observation container is empty then read the parameters that + // are necessary to fill it + if (m_obs.size() == 0) { + + // Read observation definition filename and response parameters + m_obsfile = (*this)["obsfile"].filename(); + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + + // If no observation definition file has been specified then read all + // parameters that are necessary to create an observation from scratch + // (see method setup_obs) + if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) { + m_ra = (*this)["ra"].real(); + m_dec = (*this)["dec"].real(); + m_deadc = (*this)["deadc"].real(); + m_tmin = (*this)["tmin"].real(); + m_tmax = (*this)["tmax"].real(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); } - // If there is no input filename then read all parameters that - // are required to build a model map from scratch - if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); - m_deadc = (*this)["deadc"].real(); - m_tmin = (*this)["tmin"].real(); - m_tmax = (*this)["tmax"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_enumbins = (*this)["enumbins"].integer(); - m_proj = (*this)["proj"].string(); - m_coordsys = (*this)["coordsys"].string(); - m_xref = (*this)["xref"].real(); - m_yref = (*this)["yref"].real(); - m_binsz = (*this)["binsz"].real(); - m_nxpix = (*this)["nxpix"].integer(); - m_nypix = (*this)["nypix"].integer(); + } + + // Check if we need response information + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(m_obs[i]); + if (obs != NULL) { + if (!obs->has_response()) { + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + break; + } } } - // ... otherwise, read only the parameters that are required for the - // model generation. - else { - if (m_read_ahead) { - m_outfile = (*this)["outfile"].filename(); - m_prefix = (*this)["prefix"].string(); + // If no cube file has been specified then read all parameters that + // are necessary to create the cube from scratch (see method init_cube) + if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { + m_ebinalg = (*this)["ebinalg"].string(); + if (m_ebinalg == "FILE") { + m_ebinfile = (*this)["ebinfile"].filename(); } - m_caldb = (*this)["caldb"].string(); - m_irf = (*this)["irf"].string(); - if (m_obs.models().size() == 0) { - m_srcmdl = (*this)["srcmdl"].filename(); + else { + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_enumbins = (*this)["enumbins"].integer(); } + m_proj = (*this)["proj"].string(); + m_coordsys = (*this)["coordsys"].string(); + m_xref = (*this)["xref"].real(); + m_yref = (*this)["yref"].real(); + m_binsz = (*this)["binsz"].real(); + m_nxpix = (*this)["nxpix"].integer(); + m_nypix = (*this)["nypix"].integer(); } - // Get other parameters - m_apply_edisp = (*this)["edisp"].boolean(); - // Return return; } @@ -470,37 +551,33 @@ void ctmodel::get_parameters(void) /***********************************************************************//** * @brief Setup observation container * - * @exception GException::no_cube - * No event cube found in CTA observation. - * - * This method sets up the observation container for processing. There are - * two cases: + * This method makes sure that we have an observation container with an + * associated model and that all CTA observations have a valid response + * functions. * - * If there are no observations in the actual observation container, the - * method will check in "infile" parameter. If this parameter is "NONE" or - * empty, the task parameters will be used to construct a model map. - * Otherwise, the method first tries to interpret the "infile" parameter as - * a counts map, and attemps loading of the file in an event cube. If this - * fails, the method tries to interpret the "infile" parameter as an - * observation definition XML file. If this also fails, an exception will - * be thrown. + * If on input the observation container is empty, the method sets up an + * observation container by either loading the information from an + * observation definition XML file, by loading an event list or a counts + * cube as a single observation, or by constructing an observation from + * the scratch from the task parameters. * - * If observations exist already in the observation container, the method - * will simply keep them. + * If no model exists yet in the observation container, the method will load + * the model from the file given by the task parameter "srcmdl". * - * Test if all CTA observations contain counts maps. - * - * Finally, if no models exist so far in the observation container, the - * models will be loaded from the model XML file. + * The method then loops over all observations in the container and makes + * sure that all CTA observations have a valid response function. If no + * response function exists so far, the response function will be set + * using the "caldb" and "irf" task parameters. ***************************************************************************/ void ctmodel::setup_obs(void) { - // If there are no observations in the container then try to build some + // If there are no observations in the container then get them either + // from the observation definition file or from the task parameters if (m_obs.size() == 0) { - // If no input filename has been specified, then create a model map - // from the task parameters - if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { + // If no observation definition file has been specified, then create + // an observation container from the task parameters + if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) { // Set pointing direction GCTAPointing pnt; @@ -508,27 +585,27 @@ void ctmodel::setup_obs(void) skydir.radec_deg(m_ra, m_dec); pnt.dir(skydir); - // Setup energy range covered by model - GEnergy emin(m_emin, "TeV"); - GEnergy emax(m_emax, "TeV"); - GEbounds ebds(m_enumbins, emin, emax); - // Setup time interval covered by model GGti gti; GTime tmin(m_tmin); GTime tmax(m_tmax); gti.append(tmin, tmax); - // Setup skymap - GSkymap map = GSkymap(m_proj, m_coordsys, - m_xref, m_yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, m_enumbins); - - // Create model cube from sky map - GCTAEventCube cube(map, ebds, gti); + // Setup energy range covered by model + GEnergy emin(m_emin, "TeV"); + GEnergy emax(m_emax, "TeV"); + GEbounds ebounds(1, emin, emax); - // Allocate CTA observation + // Allocate CTA observation and empty event list GCTAObservation obs; + GCTAEventList list; + + // Set event list GTI and energy boundaries and append event list + // to observation (we need in fact only the GTI for the + // computations) + list.gti(gti); + list.ebounds(ebounds); + obs.events(list); // Set CTA observation attributes obs.pointing(pnt); @@ -536,34 +613,25 @@ void ctmodel::setup_obs(void) obs.livetime(gti.ontime()*m_deadc); obs.deadc(m_deadc); - // Set event cube in observation - obs.events(cube); - // Append CTA observation to container m_obs.append(obs); - // Signal that no XML file should be used for storage - m_use_xml = false; - - } // endif: created model map from task parameters + } // endif: created single CTA observation from task parameters // ... otherwise try to load information from the file else { - // First try to open the file as a counts map + // First try to open the file as an event list or counts map try { // Allocate CTA observation GCTAObservation obs; - // Load counts map in CTA observation - obs.load(m_infile); + // Load CTA observation from FITS file + obs.load(m_obsfile); // Append CTA observation to container m_obs.append(obs); - - // Signal that no XML file should be used for storage - m_use_xml = false; } @@ -572,10 +640,7 @@ void ctmodel::setup_obs(void) // Load observations from XML file. This will throw // an exception if it fails. - m_obs.load(m_infile); - - // Signal that XML file should be used for storage - m_use_xml = true; + m_obs.load(m_obsfile); } @@ -589,8 +654,7 @@ void ctmodel::setup_obs(void) m_obs.models(GModels(m_srcmdl)); } - // Check if all CTA observations contain an event cube and setup response - // for all observations + // Setup response for all observations for (int i = 0; i < m_obs.size(); ++i) { // Is this observation a CTA observation? @@ -599,12 +663,6 @@ void ctmodel::setup_obs(void) // Yes ... if (obs != NULL) { - // Throw an exception if this observation does not contain - // an event cube - if (dynamic_cast(obs->events()) == NULL) { - throw GException::no_cube(G_SETUP_OBS); - } - // Set response if we don't have one if (!obs->has_response()) { @@ -635,211 +693,133 @@ void ctmodel::setup_obs(void) /***********************************************************************//** - * @brief Generate model map - * - * @param[in] obs CTA observation pointer. - * @param[in] models Model container. + * @brief Initialise model cube * - * @exception GException::no_cube - * No event cube found in CTA observation. + * Initialises model cube by either reading the cube definition from an + * existing cube file specified by the "infile" parameter, or by constructing + * a cube from the energy and spatial binning information provided by the + * task parameters. The GTIs of the model cube will be undefined on + * exit. ***************************************************************************/ -void ctmodel::model_map(GCTAObservation* obs, const GModels& models) +void ctmodel::init_cube(void) { - // Continue only if observation pointer is valid - if (obs != NULL) { + // Clear model cube and GTIs + m_cube.clear(); + m_gti.clear(); - // Get event cube pointer - GCTAEventCube* cube = - const_cast(dynamic_cast(obs->events())); + // If no cube file has been specified then create a cube from the task + // parameters + if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { - // Throw an exception if the observation does not hold and event - // cube - if (cube == NULL) { - throw GException::no_cube(G_MODEL_MAP); - } + // Set dummy GTI that is needed for event cube creation + GGti gti; + gti.append(GTime(0.0), GTime(0.1234)); + + // Set energy boundaries + get_ebounds(); - // Initialise statistics - double sum = 0.0; + // Setup skymap + GSkymap map = GSkymap(m_proj, m_coordsys, + m_xref, m_yref, -m_binsz, m_binsz, + m_nxpix, m_nypix, m_enumbins); - // Loop over all events in counts map - for (int i = 0; i < cube->size(); ++i) { + // Create model cube + m_cube = GCTAEventCube(map, m_ebounds, gti); - // Get event bin - GCTAEventBin* bin = (*cube)[i]; - - // Compute model value for event bin - double model = - models.eval(*(const_cast(bin)), *obs) * - bin->size(); + } - // Store value - bin->counts(model); + // ... otherwise load cube from file + else { - // Sum all events - sum += model; - } + // Load cube + m_cube.load(m_infile); - // Log results - if (logTerse()) { - log << gammalib::parformat("Model events in cube"); - log << sum << std::endl; + // Set all cube bins to zero + for (int i = 0; i < m_cube.size(); ++i) { + m_cube[i]->counts(0.0); } - } // endif: observation pointer was not valid + } // Return return; } - -/*========================================================================== - = = - = Private methods = - = = - ==========================================================================*/ /***********************************************************************//** - * @brief Initialise class members + * @brief Fill model into model cube + * + * @param[in] obs CTA observation. + * + * Adds the expected number of events for a given observation to the events + * that are already found in the model cube. The method also updates the + * GTI of the model cube so that cube GTI is a list of the GTIs of all + * observations that were used to generate the model cube. ***************************************************************************/ -void ctmodel::init_members(void) +void ctmodel::fill_cube(const GCTAObservation* obs) { - // Initialise members - m_infile.clear(); - m_outfile.clear(); - m_prefix.clear(); - m_caldb.clear(); - m_irf.clear(); - m_srcmdl.clear(); - m_proj.clear(); - m_coordsys.clear(); - m_ra = 0.0; - m_dec = 0.0; - m_deadc = 1.0; - m_tmin = 0.0; - m_tmax = 0.0; - m_emin = 0.0; - m_emax = 0.0; - m_enumbins = 0; - m_xref = 0.0; - m_yref = 0.0; - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; - m_apply_edisp = false; - - // Initialise protected members - m_obs.clear(); - m_infiles.clear(); - m_use_xml = false; - m_read_ahead = false; - - // Set logger properties - log.date(true); - - // Return - return; -} + // Continue only if observation pointer is valid + if (obs != NULL) { + // Get energy boundaries and GTI references for observation + const GEbounds& ebounds = obs->events()->ebounds(); + const GGti& gti = obs->events()->gti(); -/***********************************************************************//** - * @brief Copy class members - * - * @param[in] app Application. - ***************************************************************************/ -void ctmodel::copy_members(const ctmodel& app) -{ - // Copy attributes - m_infile = app.m_infile; - m_outfile = app.m_outfile; - m_prefix = app.m_prefix; - m_caldb = app.m_caldb; - m_irf = app.m_irf; - m_srcmdl = app.m_srcmdl; - m_proj = app.m_proj; - m_coordsys = app.m_coordsys; - m_ra = app.m_ra; - m_dec = app.m_dec; - m_deadc = app.m_deadc; - m_tmin = app.m_tmin; - m_tmax = app.m_tmax; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_enumbins = app.m_enumbins; - m_xref = app.m_xref; - m_yref = app.m_yref; - m_binsz = app.m_binsz; - m_nxpix = app.m_nxpix; - m_nypix = app.m_nypix; - m_apply_edisp = app.m_apply_edisp; + // Initialise statistics + double sum = 0.0; + int num_outside_ebds = 0; - // Copy protected members - m_obs = app.m_obs; - m_infiles = app.m_infiles; - m_use_xml = app.m_use_xml; - m_read_ahead = app.m_read_ahead; + // Setup cube GTIs for this observation + m_cube.gti(obs->events()->gti()); - // Return - return; -} + // Loop over all cube bins + for (int i = 0; i < m_cube.size(); ++i) { + // Get cube bin + GCTAEventBin* bin = m_cube[i]; -/***********************************************************************//** - * @brief Delete class members - ***************************************************************************/ -void ctmodel::free_members(void) -{ - // Write separator into logger - if (logTerse()) { - log << std::endl; - } + // Skip bin if it is outside the energy range of the observation + if (!ebounds.contains(bin->energy())) { + num_outside_ebds++; + continue; + } - // Return - return; -} + // Get actual bin value + double value = bin->counts(); + + // Compute model value for cube bin + double model = m_obs.models().eval(*bin, *obs) * bin->size(); + // Add model to actual value + value += model; + sum += model; -/***********************************************************************//** - * @brief Set output file name. - * - * @param[in] filename Input file name. - * - * Converts an input filename into an output filename by prepending the - * prefix stored in the member m_prefix to the input filename. Any path will - * be stripped from the input filename. Also a trailing ".gz" will be - * stripped. - ***************************************************************************/ -std::string ctmodel::set_outfile_name(const std::string& filename) const -{ - // Split input filename into path elements - std::vector elements = gammalib::split(filename, "/"); + // Store value + bin->counts(value); - // The last path element is the filename - std::string outname = m_prefix + elements[elements.size()-1]; + } // endfor: looped over all cube bins - // Strip any ".gz" - outname = gammalib::strip_chars(outname, ".gz"); - - // Return output filename - return outname; -} + // Append GTIs of observation to list of GTIs + m_gti.extend(gti); + // Update GTIs + m_cube.gti(m_gti); -/***********************************************************************//** - * @brief Save model map in FITS format. - * - * Save the model map as a FITS file. The filename of the FITS file is - * specified by the m_outfile member. - ***************************************************************************/ -void ctmodel::save_fits(void) -{ - // Get output filename - m_outfile = (*this)["outfile"].filename(); + // Log results + if (logTerse()) { + log << gammalib::parformat("Model events in cube"); + log << sum << std::endl; + log << gammalib::parformat("Cube bins outside energy range"); + log << num_outside_ebds << std::endl; + } - // Get CTA observation from observation container - GCTAObservation* obs = dynamic_cast(m_obs[0]); + // Log cube + if (logExplicit()) { + log.header2("Model cube"); + log << m_cube << std::endl; + } - // Save model map - save_model_map(obs, m_outfile); + } // endif: observation was valid // Return return; @@ -847,83 +827,76 @@ void ctmodel::save_fits(void) /***********************************************************************//** - * @brief Save model map(s) in XML format. + * @brief Get the energy boundaries + * + * @exception GException::invalid_value + * No valid energy boundary extension found. * - * Save the model map(s) into FITS files and write the file path information - * into a XML file. The filename of the XML file is specified by the - * m_outfile member, the filename(s) of the model map(s) are built by - * prepending the prefix given by the m_prefix member to the input model - * map(s) filenames. Any path present in the input filename will be stripped, - * i.e. the model map(s) will be written in the local working directory - * (unless a path is specified in the m_prefix member). + * Get the energy boundaries according to the user parameters. The method + * supports loading of energy boundary information from the EBOUNDS or + * ENERGYBINS extension, or setting energy boundaries using a linear + * or logarithmical spacing. + * + * @todo Ultimately this method should go in a support library as it is + * used by several ctools. ***************************************************************************/ -void ctmodel::save_xml(void) +void ctmodel::get_ebounds(void) { - // Get output filename and prefix - m_outfile = (*this)["outfile"].filename(); - m_prefix = (*this)["prefix"].string(); - - // Issue warning if output filename has no .xml suffix - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); - if (suffix != ".xml") { - log << "*** WARNING: Name of observation definition output file \""+ - m_outfile+"\"" << std::endl; - log << "*** WARNING: does not terminate with \".xml\"." << std::endl; - log << "*** WARNING: This is not an error, but might be misleading." - " It is recommended" << std::endl; - log << "*** WARNING: to use the suffix \".xml\" for observation" - " definition files." << std::endl; - } - - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Handle only CTA observations - if (obs != NULL) { - - // Set event output file name - std::string outfile = set_outfile_name(m_infiles[i]); - - // Store output file name in observation - obs->eventfile(outfile); - - // Save event list - save_model_map(obs, outfile); - - } // endif: observation was a CTA observations + // Initialse energy boundary information + m_ebounds.clear(); + + // Check whether energy binning information should be read from a FITS + // file ... + if (m_ebinalg == "FILE") { + + // Open energy boundary file using the EBOUNDS or ENERGYBINS + // extension. Throw an exception if opening fails. + GFits file(m_ebinfile); + if (file.contains("EBOUNDS")) { + file.close(); + m_ebounds.load(m_ebinfile,"EBOUNDS"); + } + else if (file.contains("ENERGYBINS")) { + file.close(); + m_ebounds.load(m_ebinfile,"ENERGYBINS"); + } + else { + file.close(); + std::string msg = "No extension with name \"EBOUNDS\" or" + " \"ENERGYBINS\" found in FITS file" + " \""+m_ebinfile+"\".\n" + "An \"EBOUNDS\" or \"ENERGYBINS\" extension" + " is required if the parameter \"ebinalg\"" + " is set to \"FILE\"."; + throw GException::invalid_value(G_GET_EBOUNDS, msg); + } + + // Set enumbins parameter to number of ebounds + m_enumbins = m_ebounds.size(); - } // endfor: looped over observations + } // endif: ebinalg was "FILE" - // Save observations in XML file - m_obs.save(m_outfile); + // ... otherwise use a linear or a logarithmically-spaced energy binning + else { - // Return - return; -} + // Initialise log mode for ebinning + bool log = true; + // check if algorithm is linear + if (m_ebinalg == "LIN") { + log = false; + } -/***********************************************************************//** - * @brief Save a single model map into a FITS file - * - * @param[in] obs Pointer to CTA observation. - * @param[in] outfile Output file name. - * - * This method saves a single model map into a FITS file. The method does - * nothing if the observation pointer is not valid. - ***************************************************************************/ -void ctmodel::save_model_map(const GCTAObservation* obs, - const std::string& outfile) const -{ - // Save only if observation is valid - if (obs != NULL) { + // todo: should we also check if m_ebinalg is "LOG" + // and throw an exception if neither LIN/LOG/FILE + // is given? - // Save observation into FITS file - obs->save(outfile, clobber()); + // Setup energy range covered by data + GEnergy emin(m_emin, "TeV"); + GEnergy emax(m_emax, "TeV"); + m_ebounds = GEbounds(m_enumbins, emin, emax, log); - } // endif: observation was valid + } //endif: ebinalg was not "FILE" // Return return; diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 88b2cea..4695d36 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctmodel - CTA counts model tool * + * ctmodel - CTA model cube generation tool * * ----------------------------------------------------------------------- * * copyright (C) 2012-2014 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctmodel.hpp - * @brief CTA counts model tool definition + * @brief CTA model cube tool definition * @author Juergen Knoedlseder */ @@ -69,25 +69,23 @@ public: void run(void); void save(void); const GObservations& obs(void) const; + const GCTAEventCube& cube(void) const; protected: // Protected methods - void init_members(void); - void copy_members(const ctmodel& app); - void free_members(void); - void get_parameters(void); - void setup_obs(void); - void model_map(GCTAObservation* obs, const GModels& models); - std::string set_outfile_name(const std::string& filename) const; - void save_fits(void); - void save_xml(void); - void save_model_map(const GCTAObservation* obs, - const std::string& outfile) const; - + void init_members(void); + void copy_members(const ctmodel& app); + void free_members(void); + void get_parameters(void); + void setup_obs(void); + void init_cube(void); + void fill_cube(const GCTAObservation* obs); + void get_ebounds(void); + // User parameters - std::string m_infile; //!< Input counts map or XML file - std::string m_outfile; //!< Output model map or XML file - std::string m_prefix; //!< Prefix for multiple model maps + std::string m_infile; //!< Input counts cube + std::string m_obsfile; //!< Event list, counts map or observation definition XML file + std::string m_outfile; //!< Output model cube std::string m_caldb; //!< Calibration database std::string m_irf; //!< Instrument response functions std::string m_srcmdl; //!< Source model @@ -96,6 +94,8 @@ protected: double m_deadc; //!< Deadtime correction factor double m_tmin; //!< Start time double m_tmax; //!< Stop time + std::string m_ebinalg; //!< Algorithm for energy binning + std::string m_ebinfile; //!< FITS-file containing energy binning double m_emin; //!< Lower energy double m_emax; //!< Upper energy int m_enumbins; //!< Number of energy bins @@ -109,10 +109,11 @@ protected: bool m_apply_edisp; //!< Apply energy dispersion? // Protected members - GObservations m_obs; //!< Observation container - std::vector m_infiles; //!< Input event filenames - bool m_use_xml; //!< Use XML file instead of FITS file - bool m_read_ahead; //!< Read ahead parameters + GObservations m_obs; //!< Observation container + GCTAEventCube m_cube; //!< Model cube + GEbounds m_ebounds; //!< Energy boundaries + GGti m_gti; //!< Model cube GTIs + bool m_read_ahead; //!< Read ahead parameters }; @@ -127,4 +128,16 @@ const GObservations& ctmodel::obs(void) const return m_obs; } + +/***********************************************************************//** + * @brief Return model cube + * + * @return Reference to model cube + ***************************************************************************/ +inline +const GCTAEventCube& ctmodel::cube(void) const +{ + return m_cube; +} + #endif /* CTMODEL_HPP */ diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index 11a9cd8..77e510d 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctmodel.par # -# Version: 00-02-00 # +# Version: 00-03-00 # # # # Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -19,22 +19,19 @@ #=================== infile, f, a, "NONE",,, "Input counts map or observation definition file" outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file" -prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file" -caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" -irf, s, a, "cta_dummy_irf",,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +edisp, b, h, no,,, "Apply energy dispersion?" # -# Instrument parameters -#====================== -ra, r, a, 83.63, 0,360, "RA of pointing (degrees)" -dec, r, a, 22.01, -90, 90, "Dec of pointing (degrees)" -deadc, r, h, 0.95, 0, 1, "Average deadtime correction factor" -# -# Time range parameters -#====================== -tmin, r, a, 0.0, 0, , "Start time (MET in s)" -tmax, r, a, 1800.0, 0, , "End time (MET in s)" +# Observation parameters +#======================= +obsfile, f, a, "NONE",,, "Observation definition file" +caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" +irf, s, a, "cta_dummy_irf",,, "Instrument response function" +ra, r, a, 83.63,0,360, "RA of pointing (degrees)" +dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)" +tmin, r, a, 0.0,0,, "Start time (MET in s)" +tmax, r, a, 1800.0,0,, "End time (MET in s)" +deadc, r, h, 0.95,0,1, "Average deadtime correction factor" # # Energy binning parameters #========================== diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 996226c..2b7ba3d 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -68,9 +68,10 @@ class Test(gammalib.GPythonTestSuite): model = ctools.ctmodel() model["infile"].filename("NONE") model["outfile"].filename("modmap.fits") + model["srcmdl"].filename(self.model_name) + model["obsfile"].filename("NONE") model["caldb"].string(self.caldb) model["irf"].string(self.irf) - model["srcmdl"].filename(self.model_name) model["ra"].real(83.63) model["dec"].real(22.01) model["tmin"].real(0.0) diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 682d7de..8853d0c 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -197,6 +197,7 @@ $ECHO -n "Test ctmodel: " # # Run 1 $ctmodel infile="cntmap.fits" \ + obsfile="cntmap.fits" \ outfile="modmap1.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ @@ -212,6 +213,7 @@ fi # # Run 3 $ctmodel infile="NONE" \ + obsfile="NONE" \ outfile="modmap3.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ @@ -241,6 +243,7 @@ fi # # Run 2 $ctmodel infile="cntmap2.fits" \ + obsfile="cntmap2.fits" \ outfile="modmap2.fits" \ caldb="irf" \ irf="cta_dummy_irf" \