Konstancja Satalecka / ctools

Commit 9fb10f84beebea498ad71a1290a6ce6c87cf379f

Authored by knodlseder
1 parent a3b2546a

initial version

Inline Side-by-side
Showing 7 changed files with 618 additions and 0 deletions
pyext/ctobssim.i 0 → 100644
  View file @9fb10f8
  1 +/***************************************************************************
  2 + * ctobssim - CTA observation simulation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file ctobssim.i
  15 + * @brief CTA observation simulation tool python interface definition
  16 + * @author J. Knodlseder
  17 + */
  18 +%{
  19 +/* Put headers and other declarations here that are needed for compilation */
  20 +#include "ctobssim.hpp"
  21 +%}
  22 +%include gammalib.i
  23 +
  24 +
  25 +/***********************************************************************//**
  26 + * @class ctobssim
  27 + *
  28 + * @brief CTA data selection tool python interface defintion.
  29 + ***************************************************************************/
  30 +class ctobssim : public GApplication {
  31 +public:
  32 + // Constructors and destructors
  33 + ctobssim(void);
  34 + ctobssim(int argc, char *argv[]);
  35 + ~ctobssim(void);
  36 +
  37 + // Methods
  38 + void run(void);
  39 + void get_parameters(void);
  40 + GPhotons simulate_photons(void);
  41 + GCTAEventList* simulate_events(const GPhotons& photons);
  42 + void save_events(const GCTAEventList* events);
  43 +};
  44 +
  45 +
  46 +/***********************************************************************//**
  47 + * @brief CTA observation simulation tool extension
  48 + ***************************************************************************/
  49 +%extend ctobssim {
  50 + ctobssim copy() {
  51 + return (*self);
  52 + }
  53 +}
... ...
src/ctobssim/ctobssim.cpp 0 → 100644
  View file @9fb10f8
  1 +/***************************************************************************
  2 + * ctobssim - CTA observation simulator tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file ctobssim.cpp
  15 + * @brief CTA observation simulator tool implementation
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +/* __ Includes ___________________________________________________________ */
  20 +#ifdef HAVE_CONFIG_H
  21 +#include <config.h>
  22 +#endif
  23 +#include <stdio.h>
  24 +#include "ctobssim.hpp"
  25 +#include "GTools.hpp"
  26 +
  27 +/* __ Debug definitions __________________________________________________ */
  28 +
  29 +/* __ Coding definitions _________________________________________________ */
  30 +
  31 +
  32 +/*==========================================================================
  33 + = =
  34 + = Constructors/destructors =
  35 + = =
  36 + ==========================================================================*/
  37 +
  38 +/***********************************************************************//**
  39 + * @brief Void constructor
  40 + ***************************************************************************/
  41 +ctobssim::ctobssim(void) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION)
  42 +{
  43 + // Initialise members
  44 + init_members();
  45 +
  46 + // Write header into logger
  47 + log_header();
  48 +
  49 + // Return
  50 + return;
  51 +}
  52 +
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Command line constructor
  56 + *
  57 + * @param[in] argc Number of arguments in command line.
  58 + * @param[in] argv Array of command line arguments.
  59 + ***************************************************************************/
  60 +ctobssim::ctobssim(int argc, char *argv[]) :
  61 + GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
  62 +{
  63 + // Initialise members
  64 + init_members();
  65 +
  66 + // Write header into logger
  67 + log_header();
  68 +
  69 + // Return
  70 + return;
  71 +}
  72 +
  73 +
  74 +/***********************************************************************//**
  75 + * @brief Destructor
  76 + ***************************************************************************/
  77 +ctobssim::~ctobssim(void)
  78 +{
  79 + // Free members
  80 + free_members();
  81 +
  82 + // Return
  83 + return;
  84 +}
  85 +
  86 +
  87 +/*==========================================================================
  88 + = =
  89 + = Operators =
  90 + = =
  91 + ==========================================================================*/
  92 +
  93 +/*==========================================================================
  94 + = =
  95 + = Public methods =
  96 + = =
  97 + ==========================================================================*/
  98 +
  99 +/***********************************************************************//**
  100 + * @brief Simulate event data
  101 + ***************************************************************************/
  102 +void ctobssim::run(void)
  103 +{
  104 + // Get parameters
  105 + get_parameters();
  106 +
  107 + // Write parameters into logger
  108 + if (logTerse()) {
  109 + log_parameters();
  110 + log << std::endl;
  111 + }
  112 +
  113 + // Write header
  114 + if (logTerse()) {
  115 + log << std::endl;
  116 + log.header1("Simulate photons");
  117 + }
  118 +
  119 + // Simulate photons
  120 + GPhotons photons = simulate_photons();
  121 +
  122 + // Write header
  123 + if (logTerse()) {
  124 + log << std::endl;
  125 + log.header1("Simulate events");
  126 + }
  127 +
  128 + // Simulate events
  129 + GCTAEventList* events = simulate_events(photons);
  130 +
  131 + // Write header
  132 + if (logTerse()) {
  133 + log << std::endl;
  134 + log.header1("Save events in FITS file");
  135 + }
  136 +
  137 + // Save events
  138 + save_events(events);
  139 +
  140 + // Release events
  141 + delete events;
  142 +
  143 + // Write separator into logger
  144 + if (logTerse())
  145 + log << std::endl;
  146 +
  147 + // Return
  148 + return;
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Get application parameters
  154 + *
  155 + * Get all task parameters from parameter file or (if required) by querying
  156 + * the user.
  157 + ***************************************************************************/
  158 +void ctobssim::get_parameters(void)
  159 +{
  160 + // Get parameters
  161 + m_infile = par("infile")->value();
  162 + m_outfile = par("outfile")->value();
  163 + m_caldb = par("caldb")->value();
  164 + m_irf = par("irf")->value();
  165 + m_seed = par("seed")->integer();
  166 + m_ra = par("ra")->real();
  167 + m_dec = par("dec")->real();
  168 + m_rad = par("rad")->real();
  169 + m_tmin = par("tmin")->real();
  170 + m_tmax = par("tmax")->real();
  171 + m_emin = par("emin")->real();
  172 + m_emax = par("emax")->real();
  173 +
  174 + // Initialise random number generator
  175 + m_ran.seed(m_seed);
  176 +
  177 + // Return
  178 + return;
  179 +}
  180 +
  181 +
  182 +/***********************************************************************//**
  183 + * @brief Simulate photons
  184 + ***************************************************************************/
  185 +GPhotons ctobssim::simulate_photons(void)
  186 +{
  187 + // Allocate MC parameters
  188 + GSkyDir dir;
  189 + GEnergy emin;
  190 + GEnergy emax;
  191 + GTime tmin;
  192 + GTime tmax;
  193 +
  194 + // Define MC parameters
  195 + dir.radec_deg(m_ra, m_dec);
  196 + emin.TeV(m_emin);
  197 + emax.TeV(m_emax);
  198 + tmin.met(m_tmin);
  199 + tmax.met(m_tmax);
  200 +
  201 + // Load models
  202 + GModels models(m_infile);
  203 +
  204 + // Simulate photons
  205 + GPhotons photons = models.mc(m_area, dir, m_rad, emin, emax, tmin, tmax, m_ran);
  206 +
  207 + // Dump simulation results
  208 + if (logNormal()) {
  209 + log << parformat("Simulated photons") << str(photons.size()) << std::endl;
  210 + log << parformat("Simulation area");
  211 + log << str(m_area) << " cm2" << std::endl;
  212 + log << parformat("Simulation cone");
  213 + log << "RA=" << dir.ra_deg() << " deg";
  214 + log << ", Dec=" << dir.dec_deg() << "deg";
  215 + log << ", r=" << m_rad << " deg" << std::endl;
  216 + log << parformat("Time interval");
  217 + log << tmin << " - " << tmax << std::endl;
  218 + log << parformat("Energy range");
  219 + log << emin << " - " << emax << std::endl;
  220 + log << parformat("Random number seed") << str(m_seed) << std::endl;
  221 + }
  222 +
  223 + // Return photons
  224 + return photons;
  225 +}
  226 +
  227 +
  228 +/***********************************************************************//**
  229 + * @brief Simulate events
  230 + *
  231 + * @param[in] photons Photon list.
  232 + *
  233 + * Simulate events from a photon list using a given response function and
  234 + * a given pointing.
  235 + ***************************************************************************/
  236 +GCTAEventList* ctobssim::simulate_events(const GPhotons& photons)
  237 +{
  238 + // Setup CTA response
  239 + GCTAResponse rsp(m_irf, m_caldb);
  240 +
  241 + // Setup CTA pointing
  242 + GCTAPointing pnt;
  243 +
  244 + // Allocate CTA event list
  245 + GCTAEventList* events = new GCTAEventList;
  246 +
  247 + // Reserves space for events
  248 + events->reserve(photons.size());
  249 +
  250 + // Simulate events
  251 + for (int i = 0; i < photons.size(); ++i) {
  252 + GCTAEventAtom* event = rsp.mc(m_area, photons[i], pnt, m_ran);
  253 + if (event != NULL)
  254 + events->append(*event);
  255 + delete event;
  256 + }
  257 +
  258 + // Dump simulation results
  259 + if (logNormal()) {
  260 + log << parformat("Simulated events") << str(events->size()) << std::endl;
  261 + }
  262 +
  263 + // Return event list
  264 + return events;
  265 +}
  266 +
  267 +
  268 +/***********************************************************************//**
  269 + * @brief Save events
  270 + *
  271 + * @todo Set header keywords
  272 + ***************************************************************************/
  273 +void ctobssim::save_events(const GCTAEventList* events)
  274 +{
  275 + // Allocate events file
  276 + GFits file;
  277 +
  278 + // Write events
  279 + events->write(&file);
  280 +
  281 + // Set keywords in EVENTS header
  282 + set_events_keywords(&file);
  283 +
  284 + // Allocate Gti
  285 + GGti gti;
  286 +
  287 + // Setup single GTI covering the selected time range
  288 + GTime tstart;
  289 + GTime tstop;
  290 + tstart.met(m_tmin);
  291 + tstop.met(m_tmax);
  292 + gti.add(tstart, tstop);
  293 +
  294 + // Write GTI
  295 + gti.write(&file);
  296 +
  297 + // Dump FITS file
  298 + if (logExplicit()) {
  299 + log << file << std::endl;
  300 + }
  301 +
  302 + // Save events file. We need to use to saveto() method to make sure
  303 + // that any existing file is indeed overwritten
  304 + file.saveto(m_outfile, clobber());
  305 +
  306 + // Close FITS file
  307 + file.close();
  308 +
  309 + // Return
  310 + return;
  311 +}
  312 +
  313 +
  314 +/***********************************************************************//**
  315 + * @brief Set keywords of EVENTS header
  316 + *
  317 + * @param[in] file FITS file.
  318 + ***************************************************************************/
  319 +void ctobssim::set_events_keywords(GFits* file)
  320 +{
  321 + // Get EVENTS HDU
  322 + GFitsHDU* hdu = file->hdu("EVENTS");
  323 +
  324 + // Set observation time information
  325 + hdu->card("TSTART", m_tmin, "[s] Mission time of start of observation");
  326 + hdu->card("TSTOP", m_tmax, "[s] Mission time of end of observation");
  327 + hdu->card("TELAPSE", m_tmax-m_tmin, "[s] Mission elapsed time");
  328 + hdu->card("ONTIME", m_tmax-m_tmin, "[s] Total good time including deadtime");
  329 + hdu->card("LIVETIME", m_tmax-m_tmin, "[s] Total livetime");
  330 +
  331 + // Set pointing information
  332 + hdu->card("RA_OBJ", m_ra, "[deg] Target Right Ascension");
  333 + hdu->card("DEC_OBJ", m_dec, "[deg] Target Declination");
  334 +
  335 + // Return
  336 + return;
  337 +}
  338 +
  339 +
  340 +/*==========================================================================
  341 + = =
  342 + = Private methods =
  343 + = =
  344 + ==========================================================================*/
  345 +
  346 +/***********************************************************************//**
  347 + * @brief Initialise class members
  348 + ***************************************************************************/
  349 +void ctobssim::init_members(void)
  350 +{
  351 + // Initialise members
  352 + m_infile.clear();
  353 + m_outfile.clear();
  354 + m_caldb.clear();
  355 + m_irf.clear();
  356 + m_seed = 1;
  357 + m_ra = 0.0;
  358 + m_dec = 0.0;
  359 + m_rad = 0.0;
  360 + m_tmin = 0.0;
  361 + m_tmax = 0.0;
  362 + m_emin = 0.0;
  363 + m_emax = 0.0;
  364 +
  365 + // Set fixed parameters
  366 + m_area = 3200000.0 * 1.0e4; //!< 3200000.0 m^2
  367 +
  368 + // Set logger properties
  369 + log.date(true);
  370 +
  371 + // Return
  372 + return;
  373 +}
  374 +
  375 +
  376 +/***********************************************************************//**
  377 + * @brief Delete class members
  378 + ***************************************************************************/
  379 +void ctobssim::free_members(void)
  380 +{
  381 + // Return
  382 + return;
  383 +}
... ...
src/ctobssim/ctobssim.hpp 0 → 100644
  View file @9fb10f8
  1 +/***************************************************************************
  2 + * ctobssim - CTA observation simulator tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file ctobssim.hpp
  15 + * @brief CTA observation simulator tool definition
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +#ifndef CTOBSSIM_HPP
  20 +#define CTOBSSIM_HPP
  21 +
  22 +/* __ Includes ___________________________________________________________ */
  23 +#include "GammaLib.hpp"
  24 +#include "GCTALib.hpp"
  25 +
  26 +/* __Definitions _________________________________________________________ */
  27 +#define CTOBSSIM_NAME "ctobssim"
  28 +#define CTOBSSIM_VERSION "v1r0p0"
  29 +
  30 +
  31 +/***********************************************************************//**
  32 + * @class ctobssim
  33 + *
  34 + * @brief CTA observation simulator tool interface defintion.
  35 + ***************************************************************************/
  36 +class ctobssim : public GApplication {
  37 +public:
  38 + // Constructors and destructors
  39 + ctobssim(void);
  40 + ctobssim(int argc, char *argv[]);
  41 + ~ctobssim(void);
  42 +
  43 + // Methods
  44 + void run(void);
  45 + void get_parameters(void);
  46 + GPhotons simulate_photons(void);
  47 + GCTAEventList* simulate_events(const GPhotons& photons);
  48 + void save_events(const GCTAEventList* events);
  49 + void set_events_keywords(GFits* file);
  50 +
  51 +protected:
  52 + // Protected methods
  53 + void init_members(void);
  54 + void free_members(void);
  55 +
  56 + // User parameters
  57 + std::string m_infile; //!< Input model
  58 + std::string m_outfile; //!< Output events file
  59 + std::string m_caldb; //!< Calibration database repository
  60 + std::string m_irf; //!< Instrument response function
  61 + int m_seed; //!< Random number generator seed
  62 + double m_ra; //!< RA of ROI centre
  63 + double m_dec; //!< DEC of ROI centre
  64 + double m_rad; //!< ROI radius
  65 + double m_tmin; //!< Start time (MET)
  66 + double m_tmax; //!< Stop time (MET)
  67 + double m_emin; //!< Lower energy (TeV)
  68 + double m_emax; //!< Upper energy (TeV)
  69 + double m_area; //!< Surface area for simulation (cm2)
  70 + GRan m_ran; //!< Random number generator
  71 +};
  72 +
  73 +#endif /* CTOBSSIM_HPP */
... ...
src/ctobssim/ctobssim.par 0 → 100644
  View file @9fb10f8
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctobssim.par #
  8 +# Version: v1r0p0 #
  9 +# #
  10 +# Author: Jurgen Knodlseder #
  11 +# knodlseder@cesr.fr #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the ctobssim application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# File information
  19 +#=================
  20 +infile, f, a, "test/data/crab.xml",,, "Model"
  21 +outfile,f, a, "events.fits",,, "Output event data file"
  22 +#
  23 +# Observation parameters
  24 +#=======================
  25 +caldb, s, a, "test/irf",,, "Calibration database"
  26 +irf, s, a, "kb_E_50h_v3",,, "Instrument response function"
  27 +#
  28 +# Simulation parameters
  29 +#======================
  30 +seed, i, h, 1, , , "Random number generator seed"
  31 +ra, r, a, 83.63, 0,360, "RA for ROI centre (degrees)"
  32 +dec, r, a, 22.01, -90, 90, "Dec for ROI centre (degrees)"
  33 +rad, r, a, 10.0, 0,180, "Radius of ROI (degrees)"
  34 +tmin, r, a, 0.0, 0, , "Start time (MET in s)"
  35 +tmax, r, a, 1800.0, 0, , "End time (MET in s)"
  36 +emin, r, a, 0.1, 0, , "Lower energy limit (TeV)"
  37 +emax, r, a, 100.0, 0, , "Upper energy limit (TeV)"
  38 +#
  39 +# Standard parameters
  40 +#====================
  41 +chatter, i, h, 2,0,4, "Chattiness of output"
  42 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  43 +debug, b, h, no,,, "Debugging mode activated"
  44 +mode, s, h, "ql",,, "Mode of automatic parameters"
... ...
src/ctobssim/main.cpp 0 → 100644
  View file @9fb10f8
  1 +/***************************************************************************
  2 + * ctobssim - CTA observation simulator tool main code *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file ctobssim.cpp
  15 + * @brief CTA observation simulator tool main code
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +/* __ Includes ___________________________________________________________ */
  20 +#ifdef HAVE_CONFIG_H
  21 +#include <config.h>
  22 +#endif
  23 +#include <stdio.h>
  24 +#include "ctobssim.hpp"
  25 +
  26 +
  27 +/***********************************************************************//**
  28 + * @brief Main entry point
  29 + *
  30 + * @param[in] argc Number of arguments
  31 + * @param[in] argv Arguments
  32 + ***************************************************************************/
  33 +int main (int argc, char *argv[])
  34 +{
  35 + // Initialise return code
  36 + int rc = 1;
  37 +
  38 + // Create instance of application
  39 + ctobssim application(argc, argv);
  40 +
  41 + // Run application
  42 + try {
  43 + application.run();
  44 + rc = 0;
  45 + }
  46 + catch (std::exception &e) {
  47 +
  48 + // Extract error message
  49 + std::string message = e.what();
  50 + std::string signal = "*** ERROR encounterted in the execution of"
  51 + " ctobssim. Run aborted ...";
  52 +
  53 + // Write error in logger
  54 + application.log << signal << std::endl;
  55 + application.log << message << std::endl;
  56 +
  57 + // Write error on standard output
  58 + std::cout << signal << std::endl;
  59 + std::cout << message << std::endl;
  60 +
  61 + } // endcatch: catched any application error
  62 +
  63 + // Return
  64 + return rc;
  65 +}
... ...
test/data/crab_cntmap.fits.gz 0 → 100644
View file @9fb10f8
No preview for this file type
test/data/crab_events.fits.gz 0 → 100644
View file @9fb10f8
No preview for this file type