Luigi Tibaldo / ctools

Commit f048b0ed5a5ea50c751f03e9af93fa7db7660422

Authored by Jürgen Knödlseder
1 parent 0604461b

Create ctmodel tool. 

The ctmodel tool generates a model map from a source model description. The definition of the model map may either be taken from a counts map, from an XMK observation definition file, or from the task parameters.
Inline Side-by-side
Showing 5 changed files with 1126 additions and 0 deletions
src/ctmodel/Makefile.am 0 → 100644
  View file @f048b0e
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Program name
  4 +bin_PROGRAMS = ctmodel
  5 +
  6 +# Sources
  7 +sources = ctmodel.cpp
  8 +
  9 +# Libtool library
  10 +noinst_LTLIBRARIES = libctmodel.la
  11 +libctmodel_la_SOURCES = $(sources)
  12 +
  13 +# Program definition
  14 +ctmodel_SOURCES = main.cpp
  15 +ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  16 +ctmodel_CPPFLAGS = @CPPFLAGS@
  17 +ctmodel_LDADD = libctmodel.la
  18 +
  19 +# Parameter files
  20 +pfilesdir = $(prefix)/syspfiles
  21 +pfiles_DATA = ctmodel.par
... ...
src/ctmodel/ctmodel.cpp 0 → 100644
  View file @f048b0e
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctmodel.cpp
  23 + * @brief CTA counts model tool implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctmodel.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_SETUP_OBS "ctmodel::setup_obs()"
  37 +#define G_MODEL_MAP "ctmodel::model_map(GCTAObservation*)"
  38 +
  39 +/* __ Debug definitions __________________________________________________ */
  40 +
  41 +/* __ Coding definitions _________________________________________________ */
  42 +
  43 +
  44 +/*==========================================================================
  45 + = =
  46 + = Constructors/destructors =
  47 + = =
  48 + ==========================================================================*/
  49 +
  50 +/***********************************************************************//**
  51 + * @brief Void constructor
  52 + ***************************************************************************/
  53 +ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
  54 +{
  55 + // Initialise members
  56 + init_members();
  57 +
  58 + // Write header into logger
  59 + log_header();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Observations constructor
  68 + *
  69 + * This method creates an instance of the class by copying an existing
  70 + * observations container.
  71 + ***************************************************************************/
  72 +ctmodel::ctmodel(GObservations obs) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Set observations
  78 + m_obs = obs;
  79 +
  80 + // Write header into logger
  81 + log_header();
  82 +
  83 + // Return
  84 + return;
  85 +}
  86 +
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Command line constructor
  91 + *
  92 + * @param[in] argc Number of arguments in command line.
  93 + * @param[in] argv Array of command line arguments.
  94 + ***************************************************************************/
  95 +ctmodel::ctmodel(int argc, char *argv[]) :
  96 + GApplication(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
  97 +{
  98 + // Initialise members
  99 + init_members();
  100 +
  101 + // Write header into logger
  102 + log_header();
  103 +
  104 + // Return
  105 + return;
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Copy constructor
  111 + *
  112 + * @param[in] app Application.
  113 + ***************************************************************************/
  114 +ctmodel::ctmodel(const ctmodel& app) : GApplication(app)
  115 +{
  116 + // Initialise members
  117 + init_members();
  118 +
  119 + // Copy members
  120 + copy_members(app);
  121 +
  122 + // Return
  123 + return;
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Destructor
  129 + ***************************************************************************/
  130 +ctmodel::~ctmodel(void)
  131 +{
  132 + // Free members
  133 + free_members();
  134 +
  135 + // Return
  136 + return;
  137 +}
  138 +
  139 +
  140 +/*==========================================================================
  141 + = =
  142 + = Operators =
  143 + = =
  144 + ==========================================================================*/
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Assignment operator
  148 + *
  149 + * @param[in] app ctmodel application.
  150 + * @return Returns ctmodel application.
  151 + ***************************************************************************/
  152 +ctmodel& ctmodel::operator= (const ctmodel& app)
  153 +{
  154 + // Execute only if object is not identical
  155 + if (this != &app) {
  156 +
  157 + // Copy base class members
  158 + this->GApplication::operator=(app);
  159 +
  160 + // Free members
  161 + free_members();
  162 +
  163 + // Initialise members
  164 + init_members();
  165 +
  166 + // Copy members
  167 + copy_members(app);
  168 +
  169 + } // endif: object was not identical
  170 +
  171 + // Return this object
  172 + return *this;
  173 +}
  174 +
  175 +
  176 +/*==========================================================================
  177 + = =
  178 + = Public methods =
  179 + = =
  180 + ==========================================================================*/
  181 +
  182 +/***********************************************************************//**
  183 + * @brief Clear instance
  184 + ***************************************************************************/
  185 +void ctmodel::clear(void)
  186 +{
  187 + // Free members
  188 + free_members();
  189 + this->GApplication::free_members();
  190 +
  191 + // Initialise members
  192 + this->GApplication::init_members();
  193 + init_members();
  194 +
  195 + // Return
  196 + return;
  197 +}
  198 +
  199 +
  200 +/***********************************************************************//**
  201 + * @brief Execute application
  202 + *
  203 + * This is the main execution method of the ctmodel class. It is invoked
  204 + * when the executable is called from command line. The method generates
  205 + * the model maps and saves the results.
  206 + ***************************************************************************/
  207 +void ctmodel::execute(void)
  208 +{
  209 + // Signal that some parameters should be read ahead
  210 + m_read_ahead = true;
  211 +
  212 + // Create the model map(s)
  213 + run();
  214 +
  215 + // Save the model map(s) into FITS file
  216 + save();
  217 +
  218 + // Return
  219 + return;
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
  224 + * @brief Generate the model map(s)
  225 + *
  226 + * This method reads the task parameters from the parfile, sets up the
  227 + * observation container, loops over all CTA observations in the container
  228 + * and generates a model map for each CTA observation.
  229 + ***************************************************************************/
  230 +void ctmodel::run(void)
  231 +{
  232 + // If we're in debug mode then all output is also dumped on the screen
  233 + if (logDebug()) {
  234 + log.cout(true);
  235 + }
  236 +
  237 + // Get task parameters
  238 + get_parameters();
  239 +
  240 + // Setup observation container
  241 + setup_obs();
  242 +
  243 + // Write parameters into logger
  244 + if (logTerse()) {
  245 + log_parameters();
  246 + log << std::endl;
  247 + }
  248 +
  249 + // Write observation(s) into logger
  250 + if (logTerse()) {
  251 + log << std::endl;
  252 + if (m_obs.size() > 1) {
  253 + log.header1("Observations");
  254 + }
  255 + else {
  256 + log.header1("Observation");
  257 + }
  258 + log << m_obs << std::endl;
  259 + }
  260 +
  261 + // Write header
  262 + if (logTerse()) {
  263 + log << std::endl;
  264 + if (m_obs.size() > 1) {
  265 + log.header1("Generate model maps");
  266 + }
  267 + else {
  268 + log.header1("Generate model map");
  269 + }
  270 + }
  271 +
  272 + // Initialise observation counter
  273 + int n_observations = 0;
  274 +
  275 + // Loop over all observations in the container
  276 + for (int i = 0; i < m_obs.size(); ++i) {
  277 +
  278 + // Initialise event input and output filenames
  279 + m_infiles.push_back("");
  280 +
  281 + // Get CTA observation
  282 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  283 +
  284 + // Continue only if observation is a CTA observation
  285 + if (obs != NULL) {
  286 +
  287 + // Write header for observation
  288 + if (logTerse()) {
  289 + if (obs->name().length() > 1) {
  290 + log.header3("Observation "+obs->name());
  291 + }
  292 + else {
  293 + log.header3("Observation");
  294 + }
  295 + }
  296 +
  297 + // Increment number of observations
  298 + n_observations++;
  299 +
  300 + // Save event file name (for possible saving)
  301 + m_infiles[i] = obs->eventfile();
  302 +
  303 + // Generate model map
  304 + model_map(obs, m_obs.models());
  305 +
  306 + } // endif: CTA observation found
  307 +
  308 + } // endfor: looped over observations
  309 +
  310 + // If more than a single observation has been handled then make sure
  311 + // that an XML file will be used for storage
  312 + if (n_observations > 1) {
  313 + m_use_xml = true;
  314 + }
  315 +
  316 + // Write observation(s) into logger
  317 + if (logTerse()) {
  318 + log << std::endl;
  319 + if (m_obs.size() > 1) {
  320 + log.header1("Observations after model map generation");
  321 + }
  322 + else {
  323 + log.header1("Observation after model map generation");
  324 + }
  325 + log << m_obs << std::endl;
  326 + }
  327 +
  328 + // Return
  329 + return;
  330 +}
  331 +
  332 +
  333 +/***********************************************************************//**
  334 + * @brief Save model map(s)
  335 + *
  336 + * This method saves the model map(s) into FITS file(s). There are two
  337 + * modes, depending on the m_use_xml flag.
  338 + *
  339 + * If m_use_xml is true, all model map(s) will be saved into FITS files,
  340 + * where the output filenames are constructued from the input filenames by
  341 + * prepending the m_prefix string to name. Any path information will be
  342 + * stripped form the input name, hence event files will be written into the
  343 + * local working directory (unless some path information is present
  344 + * in the prefix). In addition, an XML file will be created that gathers
  345 + * the filename information for the model map(s). If an XML file was present
  346 + * on input, all metadata information will be copied from this input file.
  347 + *
  348 + * If m_use_xml is false, the model map will be saved into a FITS file.
  349 + ***************************************************************************/
  350 +void ctmodel::save(void)
  351 +{
  352 + // Write header
  353 + if (logTerse()) {
  354 + log << std::endl;
  355 + if (m_obs.size() > 1) {
  356 + log.header1("Save observations");
  357 + }
  358 + else {
  359 + log.header1("Save observation");
  360 + }
  361 + }
  362 +
  363 + // Case A: Save model map(s) and XML metadata information
  364 + if (m_use_xml) {
  365 + save_xml();
  366 + }
  367 +
  368 + // Case B: Save model map as FITS file
  369 + else {
  370 + save_fits();
  371 + }
  372 +
  373 + // Return
  374 + return;
  375 +}
  376 +
  377 +
  378 +/***********************************************************************//**
  379 + * @brief Get application parameters
  380 + *
  381 + * Get all task parameters from parameter file or (if required) by querying
  382 + * the user. The parameters are read in the correct order.
  383 + ***************************************************************************/
  384 +void ctmodel::get_parameters(void)
  385 +{
  386 + // If we do not have any observations in the container then get all
  387 + // information to build (at least) one
  388 + if (m_obs.size() == 0) {
  389 +
  390 + // Read general parameters. Only read the output parameters if
  391 + // required, and also only read the model XML filename if no
  392 + // model is yet in the container.
  393 + m_infile = (*this)["infile"].filename();
  394 + if (m_read_ahead) {
  395 + m_outfile = (*this)["outfile"].filename();
  396 + m_prefix = (*this)["prefix"].string();
  397 + }
  398 + m_caldb = (*this)["caldb"].string();
  399 + m_irf = (*this)["irf"].string();
  400 + if (m_obs.models().size() == 0) {
  401 + m_srcmdl = (*this)["srcmdl"].filename();
  402 + }
  403 +
  404 + // If there is no input filename then read all parameters that
  405 + // are required to build a model map from scratch
  406 + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) {
  407 + m_tmin = (*this)["tmin"].real();
  408 + m_tmax = (*this)["tmax"].real();
  409 + m_emin = (*this)["emin"].real();
  410 + m_emax = (*this)["emax"].real();
  411 + m_enumbins = (*this)["enumbins"].integer();
  412 + m_proj = (*this)["proj"].string();
  413 + m_coordsys = (*this)["coordsys"].string();
  414 + m_xref = (*this)["xref"].real();
  415 + m_yref = (*this)["yref"].real();
  416 + m_binsz = (*this)["binsz"].real();
  417 + m_nxpix = (*this)["nxpix"].integer();
  418 + m_nypix = (*this)["nypix"].integer();
  419 + }
  420 + }
  421 +
  422 + // ... otherwise, read only the parameters that are required for the
  423 + // model generation.
  424 + else {
  425 + if (m_read_ahead) {
  426 + m_outfile = (*this)["outfile"].filename();
  427 + m_prefix = (*this)["prefix"].string();
  428 + }
  429 + m_caldb = (*this)["caldb"].string();
  430 + m_irf = (*this)["irf"].string();
  431 + if (m_obs.models().size() == 0) {
  432 + m_srcmdl = (*this)["srcmdl"].filename();
  433 + }
  434 + }
  435 +
  436 + // Return
  437 + return;
  438 +}
  439 +
  440 +
  441 +/***********************************************************************//**
  442 + * @brief Setup observation container
  443 + *
  444 + * @exception GException::no_cube
  445 + * No event cube found in CTA observation.
  446 + *
  447 + * This method sets up the observation container for processing. There are
  448 + * two cases:
  449 + *
  450 + * If there are no observations in the actual observation container, the
  451 + * method will check in "infile" parameter. If this parameter is "NONE" or
  452 + * empty, the task parameters will be used to construct a model map.
  453 + * Otherwise, the method first tries to interpret the "infile" parameter as
  454 + * a counts map, and attemps loading of the file in an event cube. If this
  455 + * fails, the method tries to interpret the "infile" parameter as an
  456 + * observation definition XML file. If this also fails, an exception will
  457 + * be thrown.
  458 + *
  459 + * If observations exist already in the observation container, the method
  460 + * will simply keep them.
  461 + *
  462 + * Test if all CTA observations contain counts maps.
  463 + *
  464 + * Finally, if no models exist so far in the observation container, the
  465 + * models will be loaded from the model XML file.
  466 + ***************************************************************************/
  467 +void ctmodel::setup_obs(void)
  468 +{
  469 + // If there are no observations in the container then try to build some
  470 + if (m_obs.size() == 0) {
  471 +
  472 + // If no input filename has been specified, then create a model map
  473 + // from the task parameters
  474 + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) {
  475 +
  476 + // Setup energy range covered by model
  477 + GEnergy emin;
  478 + GEnergy emax;
  479 + GEbounds ebds;
  480 + emin.TeV(m_emin);
  481 + emax.TeV(m_emax);
  482 + ebds.setlog(emin, emax, m_enumbins);
  483 +
  484 + // Setup time interval covered by model
  485 + GGti gti;
  486 + GTime tmin;
  487 + GTime tmax;
  488 + tmin.met(m_tmin);
  489 + tmax.met(m_tmax);
  490 + gti.append(tmin, tmax);
  491 +
  492 + // Setup skymap
  493 + GSkymap map = GSkymap(m_proj, m_coordsys,
  494 + m_xref, m_yref, m_binsz, m_binsz,
  495 + m_nxpix, m_nypix, m_enumbins);
  496 +
  497 + // Create model cube from sky map
  498 + GCTAEventCube cube(map, ebds, gti);
  499 +
  500 + // Allocate CTA observation
  501 + GCTAObservation obs;
  502 +
  503 + // Set event cube in observation
  504 + obs.events(&cube);
  505 +
  506 + // Append CTA observation to container
  507 + m_obs.append(obs);
  508 +
  509 + // Signal that no XML file should be used for storage
  510 + m_use_xml = false;
  511 +
  512 + } // endif: created model map from task parameters
  513 +
  514 + // ... otherwise try to load information from the file
  515 + else {
  516 +
  517 + // First try to open the file as a counts map
  518 + try {
  519 +
  520 + // Allocate CTA observation
  521 + GCTAObservation obs;
  522 +
  523 + // Load counts map in CTA observation
  524 + obs.load_binned(m_infile);
  525 +
  526 + // Append CTA observation to container
  527 + m_obs.append(obs);
  528 +
  529 + // Signal that no XML file should be used for storage
  530 + m_use_xml = false;
  531 +
  532 + }
  533 +
  534 + // ... otherwise try to open as XML file
  535 + catch (GException::fits_open_error &e) {
  536 +
  537 + // Load observations from XML file. This will throw
  538 + // an exception if it fails.
  539 + m_obs.load(m_infile);
  540 +
  541 + // Signal that XML file should be used for storage
  542 + m_use_xml = true;
  543 +
  544 + }
  545 +
  546 + } // endelse: loaded information from input file
  547 +
  548 + } // endif: there was no observation in the container
  549 +
  550 + // If there are no models associated with the observations then
  551 + // load now the model definition from the XML file
  552 + if (m_obs.models().size() == 0) {
  553 + m_obs.models(GModels(m_srcmdl));
  554 + }
  555 +
  556 + // Check if all CTA observations contain an event cube and setup response
  557 + // for all observations
  558 + for (int i = 0; i < m_obs.size(); ++i) {
  559 +
  560 + // Is this observation a CTA observation?
  561 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  562 +
  563 + // Yes ...
  564 + if (obs != NULL) {
  565 +
  566 + // Throw an exception if this observation does not contain
  567 + // an event cube
  568 + if (dynamic_cast<const GCTAEventCube*>(obs->events()) == NULL) {
  569 + throw GException::no_cube(G_SETUP_OBS);
  570 + }
  571 +
  572 + // Set response
  573 + obs->response(m_irf, m_caldb);
  574 +
  575 + } // endif: observation was a CTA observation
  576 + } // endfor: looped over all observations
  577 +
  578 + // Return
  579 + return;
  580 +}
  581 +
  582 +
  583 +/***********************************************************************//**
  584 + * @brief Generate model map
  585 + *
  586 + * @param[in] obs CTA observation pointer.
  587 + * @param[in] models Model container.
  588 + *
  589 + * @exception GException::no_cube
  590 + * No event cube found in CTA observation.
  591 + ***************************************************************************/
  592 +void ctmodel::model_map(GCTAObservation* obs, const GModels& models)
  593 +{
  594 + // Continue only if observation pointer is valid
  595 + if (obs != NULL) {
  596 +
  597 + // Get event cube pointer
  598 + GCTAEventCube* cube =
  599 + const_cast<GCTAEventCube*>(dynamic_cast<const GCTAEventCube*>(obs->events()));
  600 +
  601 + // Throw an exception if the observation does not hold and event
  602 + // cube
  603 + if (cube == NULL) {
  604 + throw GException::no_cube(G_MODEL_MAP);
  605 + }
  606 +
  607 + // Initialise statistics
  608 + double sum = 0.0;
  609 +
  610 + // Loop over all events in counts map
  611 + for (int i = 0; i < cube->size(); ++i) {
  612 +
  613 + // Get event bin
  614 + GCTAEventBin* bin = (*cube)[i];
  615 +
  616 + // Compute model value for event bin
  617 + double model =
  618 + models.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
  619 + bin->size();
  620 +
  621 + // Store value
  622 + bin->counts(model);
  623 +
  624 + // Sum all events
  625 + sum += model;
  626 + }
  627 +
  628 + // Log results
  629 + if (logTerse()) {
  630 + log << parformat("Model events in cube");
  631 + log << sum << std::endl;
  632 + }
  633 +
  634 + } // endif: observation pointer was not valid
  635 +
  636 + // Return
  637 + return;
  638 +}
  639 +
  640 +
  641 +/*==========================================================================
  642 + = =
  643 + = Private methods =
  644 + = =
  645 + ==========================================================================*/
  646 +
  647 +/***********************************************************************//**
  648 + * @brief Initialise class members
  649 + ***************************************************************************/
  650 +void ctmodel::init_members(void)
  651 +{
  652 + // Initialise members
  653 + m_infile.clear();
  654 + m_outfile.clear();
  655 + m_prefix.clear();
  656 + m_caldb.clear();
  657 + m_irf.clear();
  658 + m_srcmdl.clear();
  659 + m_proj.clear();
  660 + m_coordsys.clear();
  661 + m_tmin = 0.0;
  662 + m_tmax = 0.0;
  663 + m_deadc = 1.0;
  664 + m_emin = 0.0;
  665 + m_emax = 0.0;
  666 + m_enumbins = 0;
  667 + m_xref = 0.0;
  668 + m_yref = 0.0;
  669 + m_binsz = 0.0;
  670 + m_nxpix = 0;
  671 + m_nypix = 0;
  672 +
  673 + // Initialise protected members
  674 + m_obs.clear();
  675 + m_infiles.clear();
  676 + m_use_xml = false;
  677 + m_read_ahead = false;
  678 +
  679 + // Set logger properties
  680 + log.date(true);
  681 +
  682 + // Return
  683 + return;
  684 +}
  685 +
  686 +
  687 +/***********************************************************************//**
  688 + * @brief Copy class members
  689 + *
  690 + * @param[in] app Application.
  691 + ***************************************************************************/
  692 +void ctmodel::copy_members(const ctmodel& app)
  693 +{
  694 + // Copy attributes
  695 + m_infile = app.m_infile;
  696 + m_outfile = app.m_outfile;
  697 + m_prefix = app.m_prefix;
  698 + m_caldb = app.m_caldb;
  699 + m_irf = app.m_irf;
  700 + m_srcmdl = app.m_srcmdl;
  701 + m_proj = app.m_proj;
  702 + m_coordsys = app.m_coordsys;
  703 + m_tmin = app.m_tmin;
  704 + m_tmax = app.m_tmax;
  705 + m_deadc = app.m_deadc;
  706 + m_emin = app.m_emin;
  707 + m_emax = app.m_emax;
  708 + m_enumbins = app.m_enumbins;
  709 + m_xref = app.m_xref;
  710 + m_yref = app.m_yref;
  711 + m_binsz = app.m_binsz;
  712 + m_nxpix = app.m_nxpix;
  713 + m_nypix = app.m_nypix;
  714 +
  715 + // Copy protected members
  716 + m_obs = app.m_obs;
  717 + m_infiles = app.m_infiles;
  718 + m_use_xml = app.m_use_xml;
  719 + m_read_ahead = app.m_read_ahead;
  720 +
  721 + // Return
  722 + return;
  723 +}
  724 +
  725 +
  726 +/***********************************************************************//**
  727 + * @brief Delete class members
  728 + ***************************************************************************/
  729 +void ctmodel::free_members(void)
  730 +{
  731 + // Write separator into logger
  732 + if (logTerse()) {
  733 + log << std::endl;
  734 + }
  735 +
  736 + // Return
  737 + return;
  738 +}
  739 +
  740 +
  741 +/***********************************************************************//**
  742 + * @brief Set output file name.
  743 + *
  744 + * @param[in] filename Input file name.
  745 + *
  746 + * This method converts an input filename into an output filename by
  747 + * prepending the prefix specified by m_prefix to the input filename. Any
  748 + * path will be stripped from the input filename.
  749 + ***************************************************************************/
  750 +std::string ctmodel::set_outfile_name(const std::string& filename) const
  751 +{
  752 + // Split input filename into path elements
  753 + std::vector<std::string> elements = split(filename, "/");
  754 +
  755 + // The last path element is the filename
  756 + std::string outname = m_prefix + elements[elements.size()-1];
  757 +
  758 + // Return output filename
  759 + return outname;
  760 +}
  761 +
  762 +
  763 +/***********************************************************************//**
  764 + * @brief Save model map in FITS format.
  765 + *
  766 + * Save the model map as a FITS file. The filename of the FITS file is
  767 + * specified by the m_outfile member.
  768 + ***************************************************************************/
  769 +void ctmodel::save_fits(void)
  770 +{
  771 + // Get output filename
  772 + m_outfile = (*this)["outfile"].filename();
  773 +
  774 + // Get CTA observation from observation container
  775 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[0]);
  776 +
  777 + // Save model map
  778 + save_model_map(obs, m_outfile);
  779 +
  780 + // Return
  781 + return;
  782 +}
  783 +
  784 +
  785 +/***********************************************************************//**
  786 + * @brief Save model map(s) in XML format.
  787 + *
  788 + * Save the model map(s) into FITS files and write the file path information
  789 + * into a XML file. The filename of the XML file is specified by the
  790 + * m_outfile member, the filename(s) of the model map(s) are built by
  791 + * prepending the prefix given by the m_prefix member to the input model
  792 + * map(s) filenames. Any path present in the input filename will be stripped,
  793 + * i.e. the model map(s) will be written in the local working directory
  794 + * (unless a path is specified in the m_prefix member).
  795 + ***************************************************************************/
  796 +void ctmodel::save_xml(void)
  797 +{
  798 + // Get output filename and prefix
  799 + m_outfile = (*this)["outfile"].filename();
  800 + m_prefix = (*this)["prefix"].string();
  801 +
  802 + // Loop over all observation in the container
  803 + for (int i = 0; i < m_obs.size(); ++i) {
  804 +
  805 + // Get CTA observation
  806 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  807 +
  808 + // Handle only CTA observations
  809 + if (obs != NULL) {
  810 +
  811 + // Set event output file name
  812 + std::string outfile = set_outfile_name(m_infiles[i]);
  813 +
  814 + // Store output file name in observation
  815 + obs->eventfile(outfile);
  816 +
  817 + // Save event list
  818 + save_model_map(obs, outfile);
  819 +
  820 + } // endif: observation was a CTA observations
  821 +
  822 + } // endfor: looped over observations
  823 +
  824 + // Save observations in XML file
  825 + m_obs.save(m_outfile);
  826 +
  827 + // Return
  828 + return;
  829 +}
  830 +
  831 +
  832 +/***********************************************************************//**
  833 + * @brief Save a single model map into a FITS file
  834 + *
  835 + * @param[in] obs Pointer to CTA observation.
  836 + * @param[in] outfile Output file name.
  837 + *
  838 + * This method saves a single model map into a FITS file. The method does
  839 + * nothing if the observation pointer is not valid.
  840 + ***************************************************************************/
  841 +void ctmodel::save_model_map(const GCTAObservation* obs,
  842 + const std::string& outfile) const
  843 +{
  844 + // Save only if observation is valid
  845 + if (obs != NULL) {
  846 +
  847 + // Save observation into FITS file
  848 + obs->save(outfile, clobber());
  849 +
  850 + } // endif: observation was valid
  851 +
  852 + // Return
  853 + return;
  854 +}
... ...
src/ctmodel/ctmodel.hpp 0 → 100644
  View file @f048b0e
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctmodel.hpp
  23 + * @brief CTA counts model tool definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef CTMODEL_HPP
  28 +#define CTMODEL_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +
  34 +/* __Definitions _________________________________________________________ */
  35 +#define CTMODEL_NAME "ctmodel"
  36 +#define CTMODEL_VERSION "00-01-00"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class ctmodel
  41 + *
  42 + * @brief CTA counts model tool interface defintion
  43 + *
  44 + * This class creates counts model map(s). The definition of the counts model
  45 + * can be taken from a predefined observation container, from the counts maps
  46 + * found in an observation definition XML file, from an individual counts
  47 + * map, or from the parameter definition.
  48 + *
  49 + * Results are stored in an observation container that can be written to disk
  50 + * in form of FITS files (model maps) and an updated observation definition
  51 + * XML file.
  52 + ***************************************************************************/
  53 +class ctmodel : public GApplication {
  54 +
  55 +public:
  56 + // Constructors and destructors
  57 + ctmodel(void);
  58 + explicit ctmodel(GObservations obs);
  59 + ctmodel(int argc, char *argv[]);
  60 + ctmodel(const ctmodel& app);
  61 + virtual ~ctmodel(void);
  62 +
  63 + // Operators
  64 + ctmodel& operator= (const ctmodel& app);
  65 +
  66 + // Methods
  67 + void clear(void);
  68 + void execute(void);
  69 + void run(void);
  70 + void save(void);
  71 + GObservations& obs(void) { return m_obs; }
  72 + void get_parameters(void);
  73 + void setup_obs(void);
  74 + void model_map(GCTAObservation* obs, const GModels& models);
  75 +
  76 +protected:
  77 + // Protected methods
  78 + void init_members(void);
  79 + void copy_members(const ctmodel& app);
  80 + void free_members(void);
  81 + std::string set_outfile_name(const std::string& filename) const;
  82 + void save_fits(void);
  83 + void save_xml(void);
  84 + void save_model_map(const GCTAObservation* obs,
  85 + const std::string& outfile) const;
  86 +
  87 + // User parameters
  88 + std::string m_infile; //!< Input counts map or XML file
  89 + std::string m_outfile; //!< Output model map or XML file
  90 + std::string m_prefix; //!< Prefix for multiple model maps
  91 + std::string m_caldb; //!< Calibration database
  92 + std::string m_irf; //!< Instrument response functions
  93 + std::string m_srcmdl; //!< Source model
  94 + double m_tmin; //!< Start time
  95 + double m_tmax; //!< Stop time
  96 + double m_deadc; //!< Deadtime correction factor
  97 + double m_emin; //!< Lower energy
  98 + double m_emax; //!< Upper energy
  99 + int m_enumbins; //!< Number of energy bins
  100 + std::string m_proj; //!< WCS projection
  101 + std::string m_coordsys; //!< Coordinate system
  102 + double m_xref; //!< Longitude reference coordinate
  103 + double m_yref; //!< Latitude reference coordinate
  104 + double m_binsz; //!< Pixel size
  105 + int m_nxpix; //!< Number of pixels in longitude
  106 + int m_nypix; //!< Number of pixels in latitude
  107 +
  108 + // Protected members
  109 + GObservations m_obs; //!< Observation container
  110 + std::vector<std::string> m_infiles; //!< Input event filenames
  111 + bool m_use_xml; //!< Use XML file instead of FITS file
  112 + bool m_read_ahead; //!< Read ahead parameters
  113 +};
  114 +
  115 +#endif /* CTMODEL_HPP */
... ...
src/ctmodel/ctmodel.par 0 → 100644
  View file @f048b0e
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctmodel.par #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Juergen Knoedlseder #
  11 +# Jurgen.Knodlseder@irap.omp.eu #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the ctmodel application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +infile, f, a, "NONE",,, "Input counts map or observation definition file"
  21 +outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file"
  22 +prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file"
  23 +caldb, s, a, "test/irf",,, "Calibration database"
  24 +irf, s, a, "kb_E_50h_v3",,, "Instrument response function"
  25 +srcmdl, f, a, "test/data/crab.xml",,, "Source model"
  26 +#
  27 +# Time range parameters
  28 +#======================
  29 +tmin, r, a, 0.0, 0, , "Start time (MET in s)"
  30 +tmax, r, a, 1800.0, 0, , "End time (MET in s)"
  31 +deadc, r, h, 0.95, 0, 1, "Average deadtime correction factor"
  32 +#
  33 +# Energy binning parameters
  34 +#==========================
  35 +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins"
  36 +emin, r, a, 0.1,,,"Start value for first energy bin in TeV"
  37 +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV"
  38 +enumbins, i, a, 20,,,"Number of energy bins"
  39 +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
  40 +#
  41 +# Spatial binning parameters
  42 +#===========================
  43 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  44 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  45 +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
  46 +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  47 +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)"
  48 +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)"
  49 +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees"
  50 +proj, s, a, "CAR",,, "Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN"
  51 +#
  52 +# Standard parameters
  53 +#====================
  54 +chatter, i, h, 1,0,4, "Chattiness of output"
  55 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  56 +debug, b, h, no,,, "Debugging mode activated"
  57 +mode, s, h, "ql",,, "Mode of automatic parameters"
... ...
src/ctmodel/main.cpp 0 → 100644
  View file @f048b0e
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file main.cpp
  23 + * @brief CTA counts model tool main code
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctmodel.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point of application
  37 + *
  38 + * @param[in] argc Number of command line arguments.
  39 + * @param[in] argv Command line arguments.
  40 + *
  41 + * This is the main entry point of the gtbin application. It allocates a
  42 + * ctmodel object and runs the application.
  43 + ***************************************************************************/
  44 +int main (int argc, char *argv[])
  45 +{
  46 + // Initialise return code
  47 + int rc = 1;
  48 +
  49 + // Create instance of application
  50 + ctmodel application(argc, argv);
  51 +
  52 + // Run application
  53 + try {
  54 + // Execute application
  55 + application.execute();
  56 +
  57 + // Signal success
  58 + rc = 0;
  59 + }
  60 + catch (std::exception &e) {
  61 +
  62 + // Extract error message
  63 + std::string message = e.what();
  64 + std::string signal = "*** ERROR encounterted in the execution of"
  65 + " ctmodel. Run aborted ...";
  66 +
  67 + // Write error in logger
  68 + application.log << signal << std::endl;
  69 + application.log << message << std::endl;
  70 +
  71 + // Write error on standard output
  72 + std::cout << signal << std::endl;
  73 + std::cout << message << std::endl;
  74 +
  75 + } // endcatch: catched any application error
  76 +
  77 + // Return
  78 + return rc;
  79 +}
... ...