From 9fb10f84beebea498ad71a1290a6ce6c87cf379f Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Fri, 21 Jan 2011 08:02:40 +0000 Subject: [PATCH] initial version --- pyext/ctobssim.i | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctobssim/ctobssim.cpp | 383 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctobssim/ctobssim.hpp | 73 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctobssim/ctobssim.par | 44 ++++++++++++++++++++++++++++++++++++++++++++ src/ctobssim/main.cpp | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/data/crab_cntmap.fits.gz | Bin 0 -> 8803 bytes test/data/crab_events.fits.gz | Bin 0 -> 23492 bytes 7 files changed, 618 insertions(+), 0 deletions(-) create mode 100644 pyext/ctobssim.i create mode 100644 src/ctobssim/ctobssim.cpp create mode 100644 src/ctobssim/ctobssim.hpp create mode 100644 src/ctobssim/ctobssim.par create mode 100644 src/ctobssim/main.cpp create mode 100644 test/data/crab_cntmap.fits.gz create mode 100644 test/data/crab_events.fits.gz diff --git a/pyext/ctobssim.i b/pyext/ctobssim.i new file mode 100644 index 0000000..f3f4cf4 --- /dev/null +++ b/pyext/ctobssim.i @@ -0,0 +1,53 @@ +/*************************************************************************** + * ctobssim - CTA observation simulation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file ctobssim.i + * @brief CTA observation simulation tool python interface definition + * @author J. Knodlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctobssim.hpp" +%} +%include gammalib.i + + +/***********************************************************************//** + * @class ctobssim + * + * @brief CTA data selection tool python interface defintion. + ***************************************************************************/ +class ctobssim : public GApplication { +public: + // Constructors and destructors + ctobssim(void); + ctobssim(int argc, char *argv[]); + ~ctobssim(void); + + // Methods + void run(void); + void get_parameters(void); + GPhotons simulate_photons(void); + GCTAEventList* simulate_events(const GPhotons& photons); + void save_events(const GCTAEventList* events); +}; + + +/***********************************************************************//** + * @brief CTA observation simulation tool extension + ***************************************************************************/ +%extend ctobssim { + ctobssim copy() { + return (*self); + } +} diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp new file mode 100644 index 0000000..edbb5d3 --- /dev/null +++ b/src/ctobssim/ctobssim.cpp @@ -0,0 +1,383 @@ +/*************************************************************************** + * ctobssim - CTA observation simulator tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file ctobssim.cpp + * @brief CTA observation simulator tool implementation + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctobssim.hpp" +#include "GTools.hpp" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctobssim::ctobssim(void) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctobssim::ctobssim(int argc, char *argv[]) : + GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctobssim::~ctobssim(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Simulate event data + ***************************************************************************/ +void ctobssim::run(void) +{ + // Get parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Simulate photons"); + } + + // Simulate photons + GPhotons photons = simulate_photons(); + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Simulate events"); + } + + // Simulate events + GCTAEventList* events = simulate_events(photons); + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Save events in FITS file"); + } + + // Save events + save_events(events); + + // Release events + delete events; + + // Write separator into logger + if (logTerse()) + log << std::endl; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all task parameters from parameter file or (if required) by querying + * the user. + ***************************************************************************/ +void ctobssim::get_parameters(void) +{ + // Get parameters + m_infile = par("infile")->value(); + m_outfile = par("outfile")->value(); + m_caldb = par("caldb")->value(); + m_irf = par("irf")->value(); + m_seed = par("seed")->integer(); + m_ra = par("ra")->real(); + m_dec = par("dec")->real(); + m_rad = par("rad")->real(); + m_tmin = par("tmin")->real(); + m_tmax = par("tmax")->real(); + m_emin = par("emin")->real(); + m_emax = par("emax")->real(); + + // Initialise random number generator + m_ran.seed(m_seed); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Simulate photons + ***************************************************************************/ +GPhotons ctobssim::simulate_photons(void) +{ + // Allocate MC parameters + GSkyDir dir; + GEnergy emin; + GEnergy emax; + GTime tmin; + GTime tmax; + + // Define MC parameters + dir.radec_deg(m_ra, m_dec); + emin.TeV(m_emin); + emax.TeV(m_emax); + tmin.met(m_tmin); + tmax.met(m_tmax); + + // Load models + GModels models(m_infile); + + // Simulate photons + GPhotons photons = models.mc(m_area, dir, m_rad, emin, emax, tmin, tmax, m_ran); + + // Dump simulation results + if (logNormal()) { + log << parformat("Simulated photons") << str(photons.size()) << std::endl; + log << parformat("Simulation area"); + log << str(m_area) << " cm2" << std::endl; + log << parformat("Simulation cone"); + log << "RA=" << dir.ra_deg() << " deg"; + log << ", Dec=" << dir.dec_deg() << "deg"; + log << ", r=" << m_rad << " deg" << std::endl; + log << parformat("Time interval"); + log << tmin << " - " << tmax << std::endl; + log << parformat("Energy range"); + log << emin << " - " << emax << std::endl; + log << parformat("Random number seed") << str(m_seed) << std::endl; + } + + // Return photons + return photons; +} + + +/***********************************************************************//** + * @brief Simulate events + * + * @param[in] photons Photon list. + * + * Simulate events from a photon list using a given response function and + * a given pointing. + ***************************************************************************/ +GCTAEventList* ctobssim::simulate_events(const GPhotons& photons) +{ + // Setup CTA response + GCTAResponse rsp(m_irf, m_caldb); + + // Setup CTA pointing + GCTAPointing pnt; + + // Allocate CTA event list + GCTAEventList* events = new GCTAEventList; + + // Reserves space for events + events->reserve(photons.size()); + + // Simulate events + for (int i = 0; i < photons.size(); ++i) { + GCTAEventAtom* event = rsp.mc(m_area, photons[i], pnt, m_ran); + if (event != NULL) + events->append(*event); + delete event; + } + + // Dump simulation results + if (logNormal()) { + log << parformat("Simulated events") << str(events->size()) << std::endl; + } + + // Return event list + return events; +} + + +/***********************************************************************//** + * @brief Save events + * + * @todo Set header keywords + ***************************************************************************/ +void ctobssim::save_events(const GCTAEventList* events) +{ + // Allocate events file + GFits file; + + // Write events + events->write(&file); + + // Set keywords in EVENTS header + set_events_keywords(&file); + + // Allocate Gti + GGti gti; + + // Setup single GTI covering the selected time range + GTime tstart; + GTime tstop; + tstart.met(m_tmin); + tstop.met(m_tmax); + gti.add(tstart, tstop); + + // Write GTI + gti.write(&file); + + // Dump FITS file + if (logExplicit()) { + log << file << std::endl; + } + + // Save events file. We need to use to saveto() method to make sure + // that any existing file is indeed overwritten + file.saveto(m_outfile, clobber()); + + // Close FITS file + file.close(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set keywords of EVENTS header + * + * @param[in] file FITS file. + ***************************************************************************/ +void ctobssim::set_events_keywords(GFits* file) +{ + // Get EVENTS HDU + GFitsHDU* hdu = file->hdu("EVENTS"); + + // Set observation time information + hdu->card("TSTART", m_tmin, "[s] Mission time of start of observation"); + hdu->card("TSTOP", m_tmax, "[s] Mission time of end of observation"); + hdu->card("TELAPSE", m_tmax-m_tmin, "[s] Mission elapsed time"); + hdu->card("ONTIME", m_tmax-m_tmin, "[s] Total good time including deadtime"); + hdu->card("LIVETIME", m_tmax-m_tmin, "[s] Total livetime"); + + // Set pointing information + hdu->card("RA_OBJ", m_ra, "[deg] Target Right Ascension"); + hdu->card("DEC_OBJ", m_dec, "[deg] Target Declination"); + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctobssim::init_members(void) +{ + // Initialise members + m_infile.clear(); + m_outfile.clear(); + m_caldb.clear(); + m_irf.clear(); + m_seed = 1; + m_ra = 0.0; + m_dec = 0.0; + m_rad = 0.0; + m_tmin = 0.0; + m_tmax = 0.0; + m_emin = 0.0; + m_emax = 0.0; + + // Set fixed parameters + m_area = 3200000.0 * 1.0e4; //!< 3200000.0 m^2 + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctobssim::free_members(void) +{ + // Return + return; +} diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp new file mode 100644 index 0000000..a70de6c --- /dev/null +++ b/src/ctobssim/ctobssim.hpp @@ -0,0 +1,73 @@ +/*************************************************************************** + * ctobssim - CTA observation simulator tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file ctobssim.hpp + * @brief CTA observation simulator tool definition + * @author J. Knodlseder + */ + +#ifndef CTOBSSIM_HPP +#define CTOBSSIM_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTOBSSIM_NAME "ctobssim" +#define CTOBSSIM_VERSION "v1r0p0" + + +/***********************************************************************//** + * @class ctobssim + * + * @brief CTA observation simulator tool interface defintion. + ***************************************************************************/ +class ctobssim : public GApplication { +public: + // Constructors and destructors + ctobssim(void); + ctobssim(int argc, char *argv[]); + ~ctobssim(void); + + // Methods + void run(void); + void get_parameters(void); + GPhotons simulate_photons(void); + GCTAEventList* simulate_events(const GPhotons& photons); + void save_events(const GCTAEventList* events); + void set_events_keywords(GFits* file); + +protected: + // Protected methods + void init_members(void); + void free_members(void); + + // User parameters + std::string m_infile; //!< Input model + std::string m_outfile; //!< Output events file + std::string m_caldb; //!< Calibration database repository + std::string m_irf; //!< Instrument response function + int m_seed; //!< Random number generator seed + double m_ra; //!< RA of ROI centre + double m_dec; //!< DEC of ROI centre + double m_rad; //!< ROI radius + double m_tmin; //!< Start time (MET) + double m_tmax; //!< Stop time (MET) + double m_emin; //!< Lower energy (TeV) + double m_emax; //!< Upper energy (TeV) + double m_area; //!< Surface area for simulation (cm2) + GRan m_ran; //!< Random number generator +}; + +#endif /* CTOBSSIM_HPP */ diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par new file mode 100644 index 0000000..ade8454 --- /dev/null +++ b/src/ctobssim/ctobssim.par @@ -0,0 +1,44 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctobssim.par # +# Version: v1r0p0 # +# # +# Author: Jurgen Knodlseder # +# knodlseder@cesr.fr # +# IRAP # +# # +# Purpose: Parameter file for the ctobssim application # +# # +############################################################################# +# +# File information +#================= +infile, f, a, "test/data/crab.xml",,, "Model" +outfile,f, a, "events.fits",,, "Output event data file" +# +# Observation parameters +#======================= +caldb, s, a, "test/irf",,, "Calibration database" +irf, s, a, "kb_E_50h_v3",,, "Instrument response function" +# +# Simulation parameters +#====================== +seed, i, h, 1, , , "Random number generator seed" +ra, r, a, 83.63, 0,360, "RA for ROI centre (degrees)" +dec, r, a, 22.01, -90, 90, "Dec for ROI centre (degrees)" +rad, r, a, 10.0, 0,180, "Radius of ROI (degrees)" +tmin, r, a, 0.0, 0, , "Start time (MET in s)" +tmax, r, a, 1800.0, 0, , "End time (MET in s)" +emin, r, a, 0.1, 0, , "Lower energy limit (TeV)" +emax, r, a, 100.0, 0, , "Upper energy limit (TeV)" +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, "ql",,, "Mode of automatic parameters" diff --git a/src/ctobssim/main.cpp b/src/ctobssim/main.cpp new file mode 100644 index 0000000..4b6675c --- /dev/null +++ b/src/ctobssim/main.cpp @@ -0,0 +1,65 @@ +/*************************************************************************** + * ctobssim - CTA observation simulator tool main code * + * ----------------------------------------------------------------------- * + * copyright (C) 2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file ctobssim.cpp + * @brief CTA observation simulator tool main code + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctobssim.hpp" + + +/***********************************************************************//** + * @brief Main entry point + * + * @param[in] argc Number of arguments + * @param[in] argv Arguments + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctobssim application(argc, argv); + + // Run application + try { + application.run(); + rc = 0; + } + catch (std::exception &e) { + + // Extract error message + std::string message = e.what(); + std::string signal = "*** ERROR encounterted in the execution of" + " ctobssim. Run aborted ..."; + + // Write error in logger + application.log << signal << std::endl; + application.log << message << std::endl; + + // Write error on standard output + std::cout << signal << std::endl; + std::cout << message << std::endl; + + } // endcatch: catched any application error + + // Return + return rc; +} diff --git a/test/data/crab_cntmap.fits.gz b/test/data/crab_cntmap.fits.gz new file mode 100644 index 0000000..ba86920 Binary files /dev/null and b/test/data/crab_cntmap.fits.gz differ diff --git a/test/data/crab_events.fits.gz b/test/data/crab_events.fits.gz new file mode 100644 index 0000000..54196f0 Binary files /dev/null and b/test/data/crab_events.fits.gz differ -- libgit2 0.21.4