Andreas Specovius / ctools

Commit d124c4ae60b018064d108934e0baa6989f5d5034

Authored by Michael Mayer
1 parent ae29cd19

Changed ctlike to calculate errors explicitly and added first version of cttsmap

Inline Side-by-side
Showing 10 changed files with 1091 additions and 2 deletions
configure.ac
  View file @d124c4a
... ... @@ -886,6 +886,7 @@ AC_CONFIG_FILES([Makefile
886 886 src/ctpsfcube/Makefile
887 887 src/ctbkgcube/Makefile
888 888 src/ctcubemask/Makefile
  889 + src/cttsmap/Makefile
889 890 src/ctools-setup
890 891 test/Makefile
891 892 test/test_cscripts.sh
... ...
pyext/ctools/ctools.i
  View file @d124c4a
... ... @@ -55,3 +55,4 @@
55 55 %include "ctpsfcube.i"
56 56 %include "ctbkgcube.i"
57 57 %include "ctcubemask.i"
  58 +%include "cttsmap.i"
... ...
pyext/cttsmap.i 0 → 100644
  View file @d124c4a
  1 +/***************************************************************************
  2 + * cttsmap - CTA data binning tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2014 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 cttsmap.i
  23 + * @brief CTA data binning tool Python interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "cttsmap.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************
  34 +//**
  35 +* @class cttsmap
  36 +*
  37 +* @brief CTA data binning tool Python interface
  38 +***************************************************************************/
  39 +class cttsmap : public GApplication {
  40 +public:
  41 + // Constructors and destructors
  42 + cttsmap(void);
  43 + explicit cttsmap(const GObservations& obs);
  44 + cttsmap(int argc, char *argv[]);
  45 + cttsmap(const cttsmap& app);
  46 + virtual ~cttsmap(void);
  47 +
  48 + // Methods
  49 + void clear(void);
  50 + void execute(void);
  51 + void run(void);
  52 + void save(void);
  53 + const GObservations& obs(void) const;
  54 + const GSkymap& tsmap(void) const;
  55 + const GSkymap& fluxmap(void) const;
  56 + const GSkymap& indexmap(void) const;
  57 + const GSkymap& statusmap(void) const;
  58 +
  59 +
  60 +};
  61 +
  62 +
  63 +/***********************************************************************
  64 + //**
  65 +* @brief CTA data binning tool Python extension
  66 +***************************************************************************/
  67 +%extend cttsmap {
  68 + cttsmap copy() {
  69 + return (*self);
  70 + }
  71 +}
... ...
src/Makefile.am
  View file @d124c4a
... ... @@ -14,7 +14,8 @@ SUBDIRS = ctbin \
14 14 ctexpcube \
15 15 ctpsfcube \
16 16 ctbkgcube \
17   - ctcubemask
  17 + ctcubemask \
  18 + cttsmap
18 19  
19 20 # Scripts to be distributed
20 21 scripts = ctools-setup \
... ... @@ -35,7 +36,8 @@ libctools_la_LIBADD = ctbin/libctbin.la \
35 36 ctexpcube/libctexpcube.la \
36 37 ctpsfcube/libctpsfcube.la \
37 38 ctbkgcube/libctbkgcube.la \
38   - ctcubemask/libctcubemask.la
  39 + ctcubemask/libctcubemask.la \
  40 + cttsmap/libcttsmap.la
39 41  
40 42 # Distribute scripts
41 43 EXTRA_DIST = $(scripts)
... ...
src/ctlike/ctlike.cpp
  View file @d124c4a
... ... @@ -543,6 +543,9 @@ void ctlike::optimize_lm(void)
543 543 m_obs.optimize(*opt);
544 544 }
545 545  
  546 + // Compute errors
  547 + m_obs.errors(*opt);
  548 +
546 549 // Store maximum log likelihood value
547 550 m_logL = -(opt->value());
548 551  
... ...
src/cttsmap/Makefile.am 0 → 100644
  View file @d124c4a
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Program name
  4 +bin_PROGRAMS = cttsmap
  5 +
  6 +# Sources
  7 +sources = cttsmap.cpp
  8 +
  9 +# Libtool library
  10 +noinst_LTLIBRARIES = libcttsmap.la
  11 +libcttsmap_la_SOURCES = $(sources)
  12 +
  13 +# Program definition
  14 +cttsmap_SOURCES = main.cpp
  15 +cttsmap_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  16 +cttsmap_CPPFLAGS = @CPPFLAGS@
  17 +cttsmap_LDADD = libcttsmap.la
  18 +
  19 +# Parameter files
  20 +pfilesdir = $(prefix)/syspfiles
  21 +pfiles_DATA = cttsmap.par
... ...
src/cttsmap/cttsmap.cpp 0 → 100644
  View file @d124c4a
  1 +/***************************************************************************
  2 + * cttsmap - TS map computation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2014 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 cttsmap.cpp
  23 + * @brief TS map computation tool
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "cttsmap.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +
  37 +/* __ Debug definitions __________________________________________________ */
  38 +
  39 +/* __ Coding definitions _________________________________________________ */
  40 +
  41 +
  42 +/*==========================================================================
  43 + = =
  44 + = Constructors/destructors =
  45 + = =
  46 + ==========================================================================*/
  47 +
  48 +/***********************************************************************//**
  49 + * @brief Void constructor
  50 + ***************************************************************************/
  51 +cttsmap::cttsmap(void) : GApplication(CTTSMAP_NAME, CTTSMAP_VERSION)
  52 +{
  53 + // Initialise members
  54 + init_members();
  55 +
  56 + // Write header into logger
  57 + log_header();
  58 +
  59 + // Return
  60 + return;
  61 +}
  62 +
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Observations constructor
  66 + *
  67 + * param[in] obs Observation container.
  68 + *
  69 + * This method creates an instance of the class by copying an existing
  70 + * observations container.
  71 + ***************************************************************************/
  72 +cttsmap::cttsmap(const GObservations& obs) :
  73 + GApplication(CTTSMAP_NAME, CTTSMAP_VERSION)
  74 +{
  75 + // Initialise members
  76 + init_members();
  77 +
  78 + // Set observations
  79 + m_obs = obs;
  80 +
  81 + // Write header into logger
  82 + log_header();
  83 +
  84 + // Return
  85 + return;
  86 +}
  87 +
  88 +
  89 +
  90 +/***********************************************************************//**
  91 + * @brief Command line constructor
  92 + *
  93 + * @param[in] argc Number of arguments in command line.
  94 + * @param[in] argv Array of command line arguments.
  95 + ***************************************************************************/
  96 +cttsmap::cttsmap(int argc, char *argv[]) :
  97 + GApplication(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
  98 +{
  99 + // Initialise members
  100 + init_members();
  101 +
  102 + // Write header into logger
  103 + log_header();
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Copy constructor
  112 + *
  113 + * @param[in] app Application.
  114 + ***************************************************************************/
  115 +cttsmap::cttsmap(const cttsmap& app) : GApplication(app)
  116 +{
  117 + // Initialise members
  118 + init_members();
  119 +
  120 + // Copy members
  121 + copy_members(app);
  122 +
  123 + // Return
  124 + return;
  125 +}
  126 +
  127 +
  128 +/***********************************************************************//**
  129 + * @brief Destructor
  130 + ***************************************************************************/
  131 +cttsmap::~cttsmap(void)
  132 +{
  133 + // Free members
  134 + free_members();
  135 +
  136 + // Return
  137 + return;
  138 +}
  139 +
  140 +
  141 +/*==========================================================================
  142 + = =
  143 + = Operators =
  144 + = =
  145 + ==========================================================================*/
  146 +
  147 +/***********************************************************************//**
  148 + * @brief Assignment operator
  149 + *
  150 + * @param[in] app Application.
  151 + ***************************************************************************/
  152 +cttsmap& cttsmap::operator=(const cttsmap& 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 cttsmap::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 cttsmap class. It is invoked when
  204 + * the executable is called from command line.
  205 + *
  206 + * The method reads the task parameters, computed the TS values
  207 + * and stores it into map(s), and writes the results into FITS files on disk.
  208 + ***************************************************************************/
  209 +void cttsmap::execute(void)
  210 +{
  211 + // Signal that some parameters should be read ahead
  212 + m_read_ahead = true;
  213 +
  214 + // Compute the maps
  215 + run();
  216 +
  217 + // Save the maps into FITS file
  218 + save();
  219 +
  220 + // Return
  221 + return;
  222 +}
  223 +
  224 +
  225 +/***********************************************************************//**
  226 + * @brief Computes the TS maps
  227 + *
  228 + * This method moves a point-like source along a grid of coordinates
  229 + * and refits the given models including this additional source.
  230 + * The fit parameters of the source and its TS values are filled into
  231 + * respective position in the map.
  232 + ***************************************************************************/
  233 +void cttsmap::run(void)
  234 +{
  235 + // If we're in debug mode then all output is also dumped on the screen
  236 + if (logDebug()) {
  237 + log.cout(true);
  238 + }
  239 +
  240 + // Get task parameters
  241 + get_parameters();
  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("Bin observations");
  266 + }
  267 + else {
  268 + log.header1("Bin observation");
  269 + }
  270 + }
  271 +
  272 + // Initialise maps
  273 + init_maps();
  274 +
  275 + // Initialise point source model
  276 + init_testsource();
  277 +
  278 + // Get bins to be computed
  279 + // To split the computation on several jobs,
  280 + // the user is able to provide a subset of bins
  281 + // which should be computed
  282 + int binmin = (m_binmin == -1) ? 0 : m_binmin;
  283 + int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
  284 +
  285 + GOptimizerLM* opt = (logTerse()) ? new GOptimizerLM(log)
  286 + : new GOptimizerLM();
  287 +
  288 + // Get likelihood of null hypothesis if not given before
  289 + if (m_logL0 == 0.0) {
  290 + m_obs.optimize(*opt);
  291 + m_logL0 = -(opt->value());
  292 + }
  293 +
  294 + GModels models_orig = m_obs.models();
  295 + GModels models = m_obs.models();
  296 +
  297 + // Loop over all observations in the container
  298 + for (int i = binmin; i < binmax; ++i) {
  299 +
  300 + // Get the coordinate of current bin
  301 + GSkyDir bincentre = m_tsmap.inx2dir(i);
  302 +
  303 + // Verbosity option
  304 + if ((*this)["chatter"].integer()>2) {
  305 + std::cout<<"Computing TS for bin number "<<i<<" at "<<bincentre.print()<<std::endl;
  306 + }
  307 +
  308 + // Set the point source at current bin position
  309 + m_testsource["RA"].value(bincentre.ra_deg());
  310 + m_testsource["DEC"].value(bincentre.dec_deg());
  311 +
  312 + // append test source to model container
  313 + models.append(m_testsource);
  314 +
  315 + // Assign models to observation
  316 + m_obs.models(models);
  317 +
  318 + // Optimize observation container
  319 + m_obs.optimize(*opt);
  320 +
  321 + // Retrieve the Likelihood value
  322 + double logL1 = -(opt->value());
  323 +
  324 + // Compute TS value
  325 + double ts = 2.0 * (logL1-m_logL0);
  326 +
  327 + // Get test source model instance
  328 + GModels best_fit_model = m_obs.models();
  329 + GModel* testsource = best_fit_model["TestSource"];
  330 +
  331 + // Assign values to the maps
  332 + m_tsmap(i) = ts;
  333 + m_fluxmap(i) = (*testsource)["Integral"].value();
  334 + m_indexmap(i) = (*testsource)["Index"].value();
  335 + m_statusmap(i) = 1.0;
  336 +
  337 + // Remove test source from container
  338 + models.remove("TestSource");
  339 +
  340 + // Bring models to initial state
  341 + m_obs.models(models_orig);
  342 +
  343 + } // endfor: looped over bins
  344 +
  345 + // Return
  346 + return;
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief Save maps
  352 + *
  353 + * This method saves the TS map and additional maps into a FITS file.
  354 + ***************************************************************************/
  355 +void cttsmap::save(void)
  356 +{
  357 + // Write header
  358 + if (logTerse()) {
  359 + log << std::endl;
  360 + if (m_obs.size() > 1) {
  361 + log.header1("Save observations");
  362 + }
  363 + else {
  364 + log.header1("Save observation");
  365 + }
  366 + }
  367 +
  368 + // Get output filename
  369 + m_outfile = (*this)["outfile"].filename();
  370 +
  371 + // Create fits file instance
  372 + GFits fitsfile = GFits();
  373 +
  374 + // Write the sky maps to the FITS file
  375 + m_tsmap.write(fitsfile);
  376 + m_fluxmap.write(fitsfile);
  377 + m_indexmap.write(fitsfile);
  378 +
  379 + // Change names of extensions
  380 + fitsfile[0]->extname("TS MAP");
  381 + fitsfile[1]->extname("FLUX MAP");
  382 + fitsfile[2]->extname("SPECTRAL INDEX MAP");
  383 +
  384 + // Add computation log if not all bins are computed
  385 + if (m_binmin != -1 || m_binmax != -1) {
  386 + m_statusmap.write(fitsfile);
  387 + fitsfile[3]->extname("STATUS MAP");
  388 + }
  389 +
  390 + // Save FITS file
  391 + fitsfile.saveto(m_outfile,true);
  392 +
  393 + // Return
  394 + return;
  395 +}
  396 +
  397 +
  398 +/*==========================================================================
  399 + = =
  400 + = Private methods =
  401 + = =
  402 + ==========================================================================*/
  403 +
  404 +/***********************************************************************//**
  405 + * @brief Initialise class members
  406 + ***************************************************************************/
  407 +void cttsmap::init_members(void)
  408 +{
  409 + // Initialise members
  410 + m_evfile.clear();
  411 + m_outfile.clear();
  412 + m_proj.clear();
  413 + m_coordsys.clear();
  414 + m_xref = 0.0;
  415 + m_yref = 0.0;
  416 + m_binsz = 0.0;
  417 + m_nxpix = 0;
  418 + m_nypix = 0;
  419 +
  420 + // Initialise protected members
  421 + m_obs.clear();
  422 + m_read_ahead = false;
  423 + m_binmin = -1;
  424 + m_binmax = -1;
  425 + m_logL0 = 0.0;
  426 + m_tsmap.clear();
  427 + m_fluxmap.clear();
  428 + m_indexmap.clear();
  429 + m_statusmap.clear();
  430 + m_testsource.clear();
  431 + m_source_integral = 1e-12;
  432 + m_source_integral_min = 1e-16;
  433 + m_source_integral_max = 1e-8;
  434 + m_source_index = -2.5;
  435 + m_source_index_min = - 5.0;
  436 + m_source_index_max = -1.0;
  437 + m_source_emin = GEnergy(1,"TeV");
  438 + m_source_emax = GEnergy(100, "TeV");
  439 + m_free_index = false;
  440 +
  441 + // Set logger properties
  442 + log.date(true);
  443 +
  444 + // Return
  445 + return;
  446 +}
  447 +
  448 +
  449 +/***********************************************************************//**
  450 + * @brief Copy class members
  451 + *
  452 + * @param[in] app Application.
  453 + ***************************************************************************/
  454 +void cttsmap::copy_members(const cttsmap& app)
  455 +{
  456 + // Copy attributes
  457 + m_evfile = app.m_evfile;
  458 + m_outfile = app.m_outfile;
  459 + m_proj = app.m_proj;
  460 + m_coordsys = app.m_coordsys;
  461 + m_xref = app.m_xref;
  462 + m_yref = app.m_yref;
  463 + m_binsz = app.m_binsz;
  464 + m_nxpix = app.m_nxpix;
  465 + m_nypix = app.m_nypix;
  466 +
  467 + // Copy protected members
  468 + m_binmin = app.m_binmin;
  469 + m_binmax = app.m_binmax;
  470 + m_logL0 = app.m_logL0;
  471 + m_obs = app.m_obs;
  472 + m_read_ahead = app.m_read_ahead;
  473 + m_tsmap = app.m_tsmap;
  474 + m_fluxmap = app.m_fluxmap;
  475 + m_indexmap = app.m_indexmap;
  476 + m_statusmap = app.m_statusmap;
  477 + m_testsource = app.m_testsource;
  478 + m_source_integral = app.m_source_integral;
  479 + m_source_integral_min = app.m_source_integral_min;
  480 + m_source_integral_max = app.m_source_integral_max;
  481 + m_source_index = app.m_source_integral_max;
  482 + m_source_index_min = app.m_source_index_min;
  483 + m_source_index_max = app.m_source_index_max;
  484 + m_source_emin = app.m_source_emin;
  485 + m_source_emax = app.m_source_emax;
  486 + m_free_index = app.m_free_index;
  487 +
  488 + // Return
  489 + return;
  490 +}
  491 +
  492 +
  493 +/***********************************************************************//**
  494 + * @brief Delete class members
  495 + ***************************************************************************/
  496 +void cttsmap::free_members(void)
  497 +{
  498 + // Write separator into logger
  499 + if (logTerse()) {
  500 + log << std::endl;
  501 + }
  502 +
  503 + // Return
  504 + return;
  505 +}
  506 +
  507 +
  508 +/***********************************************************************//**
  509 + * @brief Get application parameters
  510 + *
  511 + * Get all task parameters from parameter file or (if required) by querying
  512 + * the user. Most parameters are only required if no observation exists so
  513 + * far in the observation container. In this case, a single CTA observation
  514 + * will be added to the container, using the definition provided in the
  515 + * parameter file.
  516 + ***************************************************************************/
  517 +void cttsmap::get_parameters(void)
  518 +{
  519 + // If there are no observations in container then add a single CTA
  520 + // observation using the parameters from the parameter file
  521 + if (m_obs.size() == 0) {
  522 +
  523 + // Get name of CTA events file
  524 + m_evfile = (*this)["evfile"].filename();
  525 +
  526 + // Allocate CTA observation
  527 + GCTAObservation obs;
  528 +
  529 + // Try first to open as FITS file
  530 + try {
  531 +
  532 + // Load event list in CTA observation
  533 + obs.load(m_evfile);
  534 +
  535 + // Append CTA observation to container
  536 + m_obs.append(obs);
  537 +
  538 + }
  539 +
  540 + // ... otherwise try to open as XML file
  541 + catch (GException::fits_open_error &e) {
  542 +
  543 + // Load observations from XML file
  544 + m_obs.load(m_evfile);
  545 +
  546 + }
  547 +
  548 + m_obs.models((*this)["srcmdl"].filename());
  549 +
  550 + } // endif: there was no observation in the container
  551 +
  552 + // Get map parameters
  553 + m_xref = (*this)["xref"].real();
  554 + m_yref = (*this)["yref"].real();
  555 + m_proj = (*this)["proj"].string();
  556 + m_coordsys = (*this)["coordsys"].string();
  557 + m_binsz = (*this)["binsz"].real();
  558 + m_nxpix = (*this)["nxpix"].integer();
  559 + m_nypix = (*this)["nypix"].integer();
  560 +
  561 + // Get test source parameters
  562 + m_source_integral = (*this)["source_integral"].real();
  563 + m_source_integral_min = (*this)["source_integral_min"].real();
  564 + m_source_integral_max = (*this)["source_integral_max"].real();
  565 + m_source_index = (*this)["source_index"].real();
  566 + m_source_index_min = (*this)["source_index_min"].real();
  567 + m_source_index_max = (*this)["source_index_max"].real();
  568 + m_source_emin = GEnergy((*this)["source_emin"].real(),"TeV");
  569 + m_source_emax = GEnergy((*this)["source_emax"].real(),"TeV");
  570 + m_free_index = (*this)["free_index"].boolean();
  571 +
  572 + // Get optional splitting parameters
  573 + m_binmin = (*this)["binmin"].integer();
  574 + m_binmax = (*this)["binmax"].integer();
  575 + m_logL0 = (*this)["logL0"].real();
  576 +
  577 + // Optionally read ahead parameters so that they get correctly
  578 + // dumped into the log file
  579 + if (m_read_ahead) {
  580 + m_outfile = (*this)["outfile"].filename();
  581 + }
  582 +
  583 + // Return
  584 + return;
  585 +}
  586 +
  587 +
  588 +/***********************************************************************//**
  589 + * @brief Initialise map information
  590 + *
  591 + * Initialises the skymaps.
  592 + ***************************************************************************/
  593 +void cttsmap::init_maps(void)
  594 +{
  595 + // Initialise map information
  596 + m_tsmap.clear();
  597 +
  598 + // Create skymap
  599 + m_tsmap = GSkymap(m_proj, m_coordsys,
  600 + m_xref, m_yref, -m_binsz, m_binsz,
  601 + m_nxpix, m_nypix, 1);
  602 +
  603 + // Initialise map information
  604 + m_fluxmap.clear();
  605 +
  606 + // Create skymap
  607 + m_fluxmap = GSkymap(m_proj, m_coordsys,
  608 + m_xref, m_yref, -m_binsz, m_binsz,
  609 + m_nxpix, m_nypix, 1);
  610 +
  611 + // Initialise map information
  612 + m_indexmap.clear();
  613 +
  614 + // Create skymap
  615 + m_indexmap = GSkymap(m_proj, m_coordsys,
  616 + m_xref, m_yref, -m_binsz, m_binsz,
  617 + m_nxpix, m_nypix, 1);
  618 +
  619 + // Initialise map information
  620 + m_statusmap.clear();
  621 +
  622 + // Create status map
  623 + m_statusmap = GSkymap(m_proj, m_coordsys,
  624 + m_xref, m_yref, -m_binsz, m_binsz,
  625 + m_nxpix, m_nypix, 1);
  626 +
  627 + // Return
  628 + return;
  629 +}
  630 +
  631 +
  632 +/***********************************************************************//**
  633 + * @brief Initialise the test source
  634 + *
  635 + * Initialises the test point source model.
  636 + ***************************************************************************/
  637 +void cttsmap::init_testsource(void)
  638 +{
  639 + // Create spatial model
  640 + GModelSpatialPointSource spat = GModelSpatialPointSource(m_tsmap.inx2dir(0));
  641 + spat["RA"].fix();
  642 + spat["DEC"].fix();
  643 +
  644 + // Create spectral model
  645 + GModelSpectralPlaw2 plaw2 = GModelSpectralPlaw2(m_source_integral,m_source_index,m_source_emin,m_source_emax);
  646 +
  647 + // Create sky model
  648 + m_testsource = GModelSky(spat,plaw2);
  649 + m_testsource.name("TestSource");
  650 +
  651 + // Set source parameter boundaries
  652 + m_testsource["Integral"].min(m_source_integral_min);
  653 + m_testsource["Integral"].max(m_source_integral_max);
  654 + m_testsource["Index"].min(m_source_index_min);
  655 + m_testsource["Index"].max(m_source_index_max);
  656 +
  657 + // Free spectral index if required
  658 + if (m_free_index) {
  659 + m_testsource["Index"].free();
  660 + }
  661 + else {
  662 + m_testsource["Index"].fix();
  663 + }
  664 +
  665 + // Return
  666 + return;
  667 +}
  668 +
... ...
src/cttsmap/cttsmap.hpp 0 → 100644
  View file @d124c4a
  1 +/***************************************************************************
  2 + * cttsmap - TS map tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2014 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 cttsmap.hpp
  23 + * @brief TS map calculation tool
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +#ifndef CTTSMAP_HPP
  28 +#define CTTSMAP_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +
  34 +/* __Definitions _________________________________________________________ */
  35 +#define CTTSMAP_NAME "cttsmap"
  36 +#define CTTSMAP_VERSION "00-01-00"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class cttsmap
  41 + *
  42 + * @brief TS map calculation interface
  43 + *
  44 + * This class computes a set of maps from any kind of observations:
  45 + * - TS map
  46 + * - flux map
  47 + * - map of spectral index
  48 + * The class operates on predefined observation containers, an individual
  49 + * event list or an observation definition XML file.
  50 + *
  51 + * During the computation a putative point-like source is moved along a
  52 + * grid of coordinates. The best fit results (TS, flux, index) are stored in maps
  53 + * which are saved in the output FITS files.
  54 + *
  55 + ***************************************************************************/
  56 +class cttsmap : public GApplication {
  57 +public:
  58 + // Constructors and destructors
  59 + cttsmap(void);
  60 + explicit cttsmap(const GObservations& obs);
  61 + cttsmap(int argc, char *argv[]);
  62 + cttsmap(const cttsmap& app);
  63 + virtual ~cttsmap(void);
  64 +
  65 + // Operators
  66 + cttsmap& operator=(const cttsmap& app);
  67 +
  68 + // Methods
  69 + void clear(void);
  70 + void execute(void);
  71 + void run(void);
  72 + void save(void);
  73 + const GObservations& obs(void) const;
  74 + const GSkymap& tsmap(void) const;
  75 + const GSkymap& fluxmap(void) const;
  76 + const GSkymap& indexmap(void) const;
  77 + const GSkymap& statusmap(void) const;
  78 +
  79 +protected:
  80 + // Protected methods
  81 + void init_members(void);
  82 + void copy_members(const cttsmap& app);
  83 + void free_members(void);
  84 + void get_parameters(void);
  85 + void init_maps(void);
  86 + void init_testsource(void);
  87 +
  88 + // User parameters
  89 + std::string m_evfile; //!< Input event list or XML file
  90 + std::string m_outfile; //!< Output counts map or XML file
  91 + std::string m_proj; //!< WCS projection
  92 + std::string m_coordsys; //!< Coordinate system
  93 + double m_xref; //!< Longitude reference coordinate
  94 + double m_yref; //!< Latitude reference coordinate
  95 + double m_binsz; //!< Pixel size
  96 + int m_nxpix; //!< Number of pixels in longitude
  97 + int m_nypix; //!< Number of pixels in latitude
  98 +
  99 + // Parameters to control speed and job splitting
  100 + int m_binmin; //!< Map bin number from which computation should start
  101 + int m_binmax; //!< Map bin number where map computation should end
  102 + double m_logL0; //!< Likelihood value of null hypothesis
  103 +
  104 + // Protected members
  105 + GObservations m_obs; //!< Observation container
  106 + bool m_read_ahead; //!< Read ahead parameters
  107 + GSkymap m_tsmap; //!< ts map
  108 + GSkymap m_fluxmap; //!< flux map
  109 + GSkymap m_indexmap; //!< spectral index map
  110 + GSkymap m_statusmap; //!< map of computed bins
  111 + GModelSky m_testsource; //!< test source for TS computation
  112 + double m_source_integral; //!< Integral flux of test source
  113 + double m_source_integral_min; //!< minimum boundary of integral flux
  114 + double m_source_integral_max; //!< maximum boundary of integral flux
  115 + double m_source_index; //!< assumed spectral index of test source
  116 + double m_source_index_min; //!< minimum boundary of spectral index of test source
  117 + double m_source_index_max; //!< maximum boundary of spectral index of test source
  118 + GEnergy m_source_emin; //!< Minimum energy for source integral
  119 + GEnergy m_source_emax; //!< Maximum energy for source integral
  120 + bool m_free_index; //!< flag if spectral index of test source is left free or fixed
  121 +
  122 +};
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Return observation container
  126 + *
  127 + * @return Reference to observation container
  128 + ***************************************************************************/
  129 +inline
  130 +const GObservations& cttsmap::obs(void) const
  131 +{
  132 + return m_obs;
  133 +}
  134 +
  135 +/***********************************************************************//**
  136 + * @brief Return TS skymap
  137 + *
  138 + * @return Reference to TS skymap
  139 + ***************************************************************************/
  140 +inline
  141 +const GSkymap& cttsmap::tsmap(void) const
  142 +{
  143 + return m_tsmap;
  144 +}
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Return flux skymap
  148 + *
  149 + * @return Reference to flux skymap
  150 + ***************************************************************************/
  151 +inline
  152 +const GSkymap& cttsmap::fluxmap(void) const
  153 +{
  154 + return m_fluxmap;
  155 +}
  156 +
  157 +/***********************************************************************//**
  158 + * @brief Return spectral index skymap
  159 + *
  160 + * @return Reference to spectral index skymap
  161 + ***************************************************************************/
  162 +inline
  163 +const GSkymap& cttsmap::indexmap(void) const
  164 +{
  165 + return m_indexmap;
  166 +}
  167 +
  168 +/***********************************************************************//**
  169 + * @brief Return status skymap
  170 + *
  171 + * @return Reference to status skymap
  172 + ***************************************************************************/
  173 +inline
  174 +const GSkymap& cttsmap::statusmap(void) const
  175 +{
  176 + return m_statusmap;
  177 +}
  178 +
  179 +
  180 +#endif /* CTTSMAP_HPP */
... ...
src/cttsmap/cttsmap.par 0 → 100644
  View file @d124c4a
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: cttsmap.par #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@desy.de #
  12 +# DESY #
  13 +# #
  14 +# Purpose: Parameter file for the cttsmap application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +evfile, f, a, "events.fits",,,"Input event list or observation definition file"
  21 +srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +outfile, f, a, "tsmap.fits",,, "Output counts cube"
  23 +
  24 +#
  25 +# Spatial binning parameters
  26 +#===========================
  27 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  28 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  29 +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
  30 +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  31 +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)"
  32 +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)"
  33 +proj, s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
  34 +
  35 +#
  36 +# Test source parameters
  37 +#===========================
  38 +source_integral, r, h, 1e-10,,, "Integral start parameter [cm-2 s-1]"
  39 +source_integral_min, r, h, 1e-16,,, "Integral minimum boundary [cm-2 s-1]"
  40 +source_integral_max, r, h, 1e-3,,, "Integral maximum boundary [cm-2 s-1]"
  41 +source_index, r, h, -2.5,,, "Index start parameter"
  42 +source_index_min, r, h, -5.0,,, "Index minimum boundary"
  43 +source_index_max, r, h, -1.0,,, "Index maximum boundary"
  44 +source_emin, r, h, 1.0,,, "Minimum energy for flux computation"
  45 +source_emax, r, h, 100.0,,, "Maximum energy for flux computation"
  46 +free_index, b, h, no,,, "leave source index free?"
  47 +
  48 +#
  49 +# Parameters for splitting and speed purpose
  50 +#===========================
  51 +binmin, i, h, -1,,, "First bin to compute"
  52 +binmax, i, h, -1,,, "Last bin to compute"
  53 +logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis"
  54 +
  55 +
  56 +#
  57 +# Standard parameters
  58 +#====================
  59 +chatter, i, h, 2,0,4, "Chattiness of output"
  60 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  61 +debug, b, h, no,,, "Debugging mode activated"
  62 +mode, s, h, "ql",,, "Mode of automatic parameters"
  63 +logfile, f, h, "ctbin.log",,, "Log filename"
... ...
src/cttsmap/main.cpp 0 → 100644
  View file @d124c4a
  1 +/***************************************************************************
  2 + * ctbin - CTA data binning main code *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
  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 ctbin/main.cpp
  23 + * @brief CTA data binning tool main code
  24 + * @author J. Knodlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "cttsmap.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 cttsmap application. It allocates a
  42 + * cttsmap 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 + cttsmap 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 + " cttsmap. 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 +}
... ...