diff --git a/.gitignore b/.gitignore index 83fc7cc..6608ac2 100644 --- a/.gitignore +++ b/.gitignore @@ -59,6 +59,7 @@ src/ctselect/ctselect src/ctskymap/ctskymap src/ctexpcube/ctexpcube src/ctpsfcube/ctpsfcube +src/ctedispcube/ctedispcube src/ctbkgcube/ctbkgcube src/ctcubemask/ctcubemask src/ctbutterfly/ctbutterfly diff --git a/configure.ac b/configure.ac index 78f10ca..80dcdbf 100644 --- a/configure.ac +++ b/configure.ac @@ -960,6 +960,7 @@ AC_CONFIG_FILES([Makefile src/ctselect/Makefile src/ctexpcube/Makefile src/ctpsfcube/Makefile + src/ctedispcube/Makefile src/ctbkgcube/Makefile src/ctcubemask/Makefile src/cttsmap/Makefile diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 6b778b9..3a8aa95 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -40,7 +40,7 @@ class cslightcrv(ctools.cscript): # Set name self.name = "cslightcrv" - self.version = "1.0.0" + self.version = "1.1.0" # Initialise some members self.obs = None @@ -104,7 +104,9 @@ class cslightcrv(ctools.cscript): pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name")) pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) pars.append(gammalib.GApplicationPar("outfile","f","a","lightcurve.fits","","","Output light curve file")) @@ -190,6 +192,9 @@ class cslightcrv(ctools.cscript): # Get source name self.m_srcname = self["srcname"].string() + # Get energy dispersion flag + self.m_edisp = self["edisp"].boolean() + # Get time boundaries self.create_tbounds() @@ -430,7 +435,7 @@ class cslightcrv(ctools.cscript): # Header if self.logTerse(): - self.log.header3("Creating PSF cube") + self.log.header3("Creating point spread function cube") # Create psf cube psfcube = ctools.ctpsfcube(select.obs()) @@ -449,6 +454,30 @@ class cslightcrv(ctools.cscript): psfcube["proj"] = self.m_proj psfcube.run() + # Check if we need to include energy dispersion + if self.m_edisp: + + # Header + if self.logTerse(): + self.log.header3("Creating energy dispersion cube") + + # Create edisp cube + edispcube = ctools.ctedispcube(select.obs()) + edispcube["incube"] = "NONE" + edispcube["usepnt"] = False + edispcube["ebinalg"] = "LOG" + edispcube["xref"] = self.m_xref + edispcube["yref"] = self.m_yref + edispcube["binsz"] = self.m_binsz + edispcube["nxpix"] = self.m_nxpix + edispcube["nypix"] = self.m_nypix + edispcube["enumbins"] = self.m_ebins + edispcube["emin"] = self.m_emin + edispcube["emax"] = self.m_emax + edispcube["coordsys"] = self.m_coordsys + edispcube["proj"] = self.m_proj + edispcube.run() + # Header if self.logTerse(): self.log.header3("Creating background cube") @@ -472,12 +501,17 @@ class cslightcrv(ctools.cscript): # Set new binned observation obs = bin.obs() - - # Set precomputed binned response - obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) - + # Get new models models = bkgcube.models() + + # Set precomputed binned response + if self.m_edisp: + obs[0].response(expcube.expcube(), psfcube.psfcube(), + edispcube.edispcube(), bkgcube.bkgcube()) + else: + obs[0].response(expcube.expcube(), psfcube.psfcube(), + bkgcube.bkgcube()) # Fix background models if required if self.m_fix_bkg: @@ -495,6 +529,7 @@ class cslightcrv(ctools.cscript): # Likelihood like = ctools.ctlike(obs) + like["edisp"] = self.m_edisp like.run() # Skip bin if no event was present diff --git a/cscripts/cspull.py b/cscripts/cspull.py index beaf0a1..32fac3b 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -52,6 +52,7 @@ class cspull(ctools.cscript): self.m_exposure = None self.m_psfcube = None self.m_bckcube = None + self.m_edispcube = None self.m_stackmodels = None self.m_coordsys = "CEL" self.m_proj = "TAN" @@ -250,6 +251,29 @@ class cspull(ctools.cscript): if self.logTerse(): self.log("Computed Psf cube\n") + # Optionally get stacked Edisp + if self.m_edisp: + edispcube = ctools.edispcube(self.obs) + edispcube["incube"] = "NONE" + edispcube["usepnt"] = True + edispcube["ebinalg"] = "LOG" + edispcube["binsz"] = self.m_binsz*10.0 + edispcube["nxpix"] = self.m_npix/10 + edispcube["nypix"] = self.m_npix/10 + edispcube["enumbins"] = self.m_enumbins + edispcube["emin"] = self["emin"].real() + edispcube["emax"] = self["emax"].real() + edispcube["coordsys"] = self.m_coordsys + edispcube["proj"] = self.m_proj + edispcube.run() + + # Store result + self.m_edispcube = edispcube.edispcube().copy() + + # Logging + if self.logTerse(): + self.log("Computed energy dispersion cube\n") + # Get stacked background bkgcube = ctools.ctbkgcube(self.obs) bkgcube["incube"] = "NONE" @@ -325,7 +349,6 @@ class cspull(ctools.cscript): for n in result['colnames']: headers[n] = n writer.writerow(headers) - #writer.writerow(dict((_,_) for _ in result['colnames'])) else: file = open(self.m_outfile.url(), 'a') writer = csv.DictWriter(file, result['colnames']) @@ -391,7 +414,12 @@ class cspull(ctools.cscript): # If stacked, add stacked responses and model if self.obs.size() > 1 and self.m_enumbins > 0: - obs[0].response(self.m_exposure, self.m_psfcube, self.m_bckcube) + if self.m_edisp: + obs[0].response(self.m_exposure, self.m_psfcube, + self.m_edispcube, self.m_bckcube) + else: + obs[0].response(self.m_exposure, self.m_psfcube, + self.m_bckcube) obs.models(self.m_stackmodels) # Determine number of events in simulation diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 4d3968e..e1fab98 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -108,6 +108,7 @@ class csresmap(ctools.cscript): pars.append(gammalib.GApplicationPar("modcube","f","a","NONE","","","Input model cube file (generated with ctmodel)")) pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Input background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) @@ -174,7 +175,6 @@ class csresmap(ctools.cscript): # Skip ctbin step later on self.m_skip_binning = True - # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) @@ -371,7 +371,6 @@ class csresmap(ctools.cscript): # Get model map into GSkyMap object modelmap = model.cube().map().copy() - # Store counts map as residual map. Note that we need a # special construct here to avoid memory leaks. This seems # to be a SWIG feature as SWIG creates a new object when diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 9d18b7b..89a4f49 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -2,7 +2,7 @@ # ========================================================================== # Spectral points generation script. # -# Copyright (C) 2014-2015 Michael Mayer +# Copyright (C) 2014-2016 Michael Mayer # # 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 @@ -40,7 +40,7 @@ class csspec(ctools.cscript): # Set name self.name = "csspec" - self.version = "1.0.0" + self.version = "1.1.0" # Initialise some members self.obs = gammalib.GObservations() @@ -105,6 +105,7 @@ class csspec(ctools.cscript): pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name")) pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Input background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index 662573b..952ff44 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -42,7 +42,7 @@ class cstsdist(ctools.cscript): """ # Set name self.name = "cstsdist" - self.version = "1.0.0" + self.version = "1.1.0" # Initialise some members self.obs = None @@ -99,6 +99,7 @@ class cstsdist(ctools.cscript): pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name")) pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) diff --git a/doc/source/reference_manual/ctedispcube.rst b/doc/source/reference_manual/ctedispcube.rst new file mode 100644 index 0000000..c04a4e2 --- /dev/null +++ b/doc/source/reference_manual/ctedispcube.rst @@ -0,0 +1,138 @@ +.. _ctedispcube: + +ctedispcube +========= + +Generate energy dispersion cube for a counts cube. + + +Synopsis +-------- + +This tool generates an energy dispersion cube for a counts cube. An energy +dispersion cube is a 4-dimensional cube spanned by Right Ascension or +Galactic longitude, Declination or Galactic latitude, energy, and fraction +between recontructed and true photon energy. The energy binning of the cube +may be either linear, logarithmic, or custom defined using an input file. + +ctedispcube requires on input the event list or observation definition file +that has been used in the generation of the counts cube using :doc:`ctbin`. + +It is not recommended to use the counts cube for the energy dispersion +cube definition, although this is formally possible by specifying the counts +cube as ``incube`` parameter. This leads however to a large FITS file on +output since the number of bins in the counts cube will be multiplied by +the number of offset angle bins (typically 200). Since the energy dispersion +varies only little over the field of view of the camera it is +recommended to use a rather coarse spatial binning to keep the file size +manageable (with a typical value of ``binsz=1.0``). + +ctedispcube generates an energy dispersion cube FITS file comprising three +extensions. The primary extension contains a 3-dimensional image that contains +the energy disperison values. The energy and fraction dimensions of the enrgy +dispersion cube are folded into the 3rd dimension of the FITS image. The next +extension named ``EBOUNDS`` contains a binary table that defines the energy +boundaries of the exposure cube. The last extension named ``DELTAS`` contains +a binary table that defines the fraction between recontructed and true photon +energy. + + +General parameters +------------------ + +``inobs [file]`` + Input event list or observation definition XML file. + +``incube [file]`` + Input counts cube file to extract energy dispersion cube definition. + +``caldb [string]`` + Calibration database. + +``irf [string]`` + Instrument response function. + +``(edisp = no) [boolean]`` + Apply energy dispersion for response computation. + +``outcube [file]`` + Output energy dispersion cube file. + +``ebinalg [string]`` + Algorithm for defining energy bins. + +``emin [real]`` + Lower energy value for first energy bin (in TeV). + +``emax [real]`` + Upper energy value for last energy bin (in TeV). + +``enumbins [integer]`` + Number of energy bins. + +``ebinfile [file]`` + Name of the file containing the energy bin definition. + +``(usepnt = no) [boolean]`` + Use CTA pointing direction for cube centre instead of xref/yref parameters? + +``nxpix [integer]`` + Number of cube bins in Right Ascension or Galactic longitude. + +``nypix [integer]`` + Number of cube bins in Declination or Galactic latitude. + +``binsz [real]`` + Cube bin size (in degrees/pixel). + +``coordsys [string]`` + Coordinate system (CEL - celestial, GAL - galactic). + +``xref [real]`` + Right Ascension / Galactic longitude of cube centre (J2000, in degrees). + +``yref [real]`` + Declination / Galactic latitude of cube centre (J2000, in degrees). + +``proj [string]`` + Projection method. + +``(mmax = 10.0) [real]`` + Upper bound of fraction between recontructed and true photon energy. + +``(migrabins = 200) [integer]`` + Number of migration bins. + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output energy dispersion cube file should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = ctedispcube.log) [string]`` + Name of log file. + + +Related tools or scripts +------------------------ + +:doc:`ctbin` diff --git a/doc/source/reference_manual/ctpsfcube.rst b/doc/source/reference_manual/ctpsfcube.rst index 254a82d..8871c0c 100644 --- a/doc/source/reference_manual/ctpsfcube.rst +++ b/doc/source/reference_manual/ctpsfcube.rst @@ -20,7 +20,7 @@ ctpsfcube requires on input the event list or observation definition file that has been used in the generation of the counts cube using :doc:`ctbin`. It is not recommended to use the counts cube for the point spread function -cube definition, although this is formally possibly by specifying the counts +cube definition, although this is formally possible by specifying the counts cube as ``incube`` parameter. This leads however to a large FITS file on output since the number of bins in the counts cube will be multiplied by the number of offset angle bins (typically 200). Since the point spread diff --git a/doc/source/reference_manual/index.rst b/doc/source/reference_manual/index.rst index bc42824..0af6d05 100644 --- a/doc/source/reference_manual/index.rst +++ b/doc/source/reference_manual/index.rst @@ -19,6 +19,7 @@ ctools ctbkgcube --- Generates background cube ctbutterfly --- Compute butterfly ctcubemask --- Filter counts cube + ctedispcube --- Generates energy dispersion cube ctexpcube --- Generates exposure cube ctlike --- Performs maximum likelihood fitting ctmodel --- Computes model counts cube diff --git a/src/Makefile.am b/src/Makefile.am index 25b1a9b..dddc604 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -14,6 +14,7 @@ SUBDIRS = support \ ctskymap \ ctexpcube \ ctpsfcube \ + ctedispcube \ ctbkgcube \ ctcubemask \ cttsmap \ @@ -40,6 +41,7 @@ libctools_la_LIBADD = support/libsupport.la \ ctskymap/libctskymap.la \ ctexpcube/libctexpcube.la \ ctpsfcube/libctpsfcube.la \ + ctedispcube/libctedispcube.la \ ctbkgcube/libctbkgcube.la \ ctcubemask/libctcubemask.la \ cttsmap/libcttsmap.la \ diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index 3e20ea9..d6ecad8 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctbutterfly # -# Version: 1.0.0 # +# Version: 1.1.0 # # # # Author: Michael Mayer # # michael.mayer@physik.hu-berlin.de # @@ -22,6 +22,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" srcname, s, a, "Crab",,,"Source of interest" expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube, f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" diff --git a/src/ctedispcube/Makefile.am b/src/ctedispcube/Makefile.am new file mode 100644 index 0000000..bf905d9 --- /dev/null +++ b/src/ctedispcube/Makefile.am @@ -0,0 +1,25 @@ +## Process this file with automake to produce Makefile.in + +# Set include directories +INCLUDES = -I$(top_srcdir)/src/support + +# Program name +bin_PROGRAMS = ctedispcube + +# Sources +sources = ctedispcube.cpp + +# Libtool library +noinst_LTLIBRARIES = libctedispcube.la +libctedispcube_la_SOURCES = $(sources) +libctedispcube_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@ + +# Program definition +ctedispcube_SOURCES = main.cpp +ctedispcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctedispcube_CPPFLAGS = @CPPFLAGS@ +ctedispcube_LDADD = libctedispcube.la $(top_srcdir)/src/support/libsupport.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctedispcube.par diff --git a/src/ctedispcube/ctedispcube.cpp b/src/ctedispcube/ctedispcube.cpp new file mode 100644 index 0000000..3bc2c8c --- /dev/null +++ b/src/ctedispcube/ctedispcube.cpp @@ -0,0 +1,401 @@ +/*************************************************************************** + * ctedispcube - Energy dispersion cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctedispcube.cpp + * @brief Energy dispersion cube generation tool implementation + * @author Maria Haupt + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctedispcube.hpp" +#include "GTools.hpp" +#include "GWcs.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctedispcube::get_parameters()" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs an empty energy dispersion tool. + ***************************************************************************/ +ctedispcube::ctedispcube(void) : ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * @param[in] obs Observation container. + * + * Constructs an energy dispersion tool from the information that is provided + * in an observation container @p obs. + ***************************************************************************/ +ctedispcube::ctedispcube(const GObservations& obs) : + ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + * + * Constructs an energy dispersion tool by parsing the arguments provided + * on the command line. + ***************************************************************************/ +ctedispcube::ctedispcube(int argc, char *argv[]) : + ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Energy dispersion tool. + * + * Constructs an energy dispersion tool by copying anothere energy + * dispersion tool. + ***************************************************************************/ +ctedispcube::ctedispcube(const ctedispcube& app) : ctool(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + * + * Desctructs an energy dispersion tool. + ***************************************************************************/ +ctedispcube::~ctedispcube(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Energy dispersion tool. + * @return Energy dispersion tool. + * + * Assigns energy dispersion tool. + ***************************************************************************/ +ctedispcube& ctedispcube::operator=(const ctedispcube& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->ctool::operator=(app); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear energy dispersion tool + * + * Set the energy disperison tool to an empty tool. + ***************************************************************************/ +void ctedispcube::clear(void) +{ + // Free members + free_members(); + this->ctool::free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + this->ctool::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate the energy dispersion cube + * + * Generates the energy dispersion cube by looping over all unbinned CTA + * observations in the observation container. + ***************************************************************************/ +void ctedispcube::run(void) +{ + // If we're in debug mode then all output is also dumped on the screen + if (logDebug()) { + log.cout(true); + } + + // Get task parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Observations"); + } + else { + log.header1("Observation"); + } + log << m_obs << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Generate energy dispersion cube"); + } + + // Fill Edisp + m_edispcube.fill(m_obs, &log); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save energy dispersion cube + * + * Saves the energy dispersion cube in a FITS file. The FITS filename is + * provided by the "outcube" parameter. + ***************************************************************************/ +void ctedispcube::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Save Edisp cube"); + } + + // Get output filename + m_outcube = (*this)["outcube"].filename(); + + // Save only if filename is non-empty + if (!m_outcube.is_empty()) { + + // Log filename + if (logTerse()) { + log << "Save energy dispersion cube into file \""+m_outcube+"\"."; + log << std::endl; + } + + // Save energy dispersion cube + m_edispcube.save(m_outcube, clobber()); + + } + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctedispcube::init_members(void) +{ + // Initialise members + m_outcube.clear(); + m_obs.clear(); + m_edispcube.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Energy dispersion cube. + ***************************************************************************/ +void ctedispcube::copy_members(const ctedispcube& app) +{ + // Copy members + m_outcube = app.m_outcube; + m_obs = app.m_obs; + m_edispcube = app.m_edispcube; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctedispcube::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all task parameters from parameter file or (if required) by querying + * the user. The parameters are read in the correct order. + ***************************************************************************/ +void ctedispcube::get_parameters(void) +{ + // If there are no observations in container then load them via user + // parameters + if (m_obs.size() == 0) { + + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); + + // Throw exception if counts cube is given + require_inobs_nocube(G_GET_PARAMETERS); + + // Build observation container + m_obs = get_observations(); + + } // endif: there was no observation in the container + + // Get the incube filename + std::string incube = (*this)["incube"].filename(); + + // Get additional binning parameters + double migramax = (*this)["migramax"].real(); + int migrabins = (*this)["migrabins"].integer(); + + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { + + // Create an event cube based on task parameters + GCTAEventCube cube = create_cube(m_obs); + + // Define edisp cube + m_edispcube = GCTACubeEdisp(cube, migramax, migrabins); + + } + + // ... otherwise setup the exposure cube from the counts map + else { + + // Load event cube from filename + GCTAEventCube cube(incube); + + // Define Edisp cube + m_edispcube = GCTACubeEdisp(cube, migramax, migrabins); + + } // endelse: cube loaded from file + + // Read output filename (if needed) + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); + } + + // Return + return; +} diff --git a/src/ctedispcube/ctedispcube.hpp b/src/ctedispcube/ctedispcube.hpp new file mode 100644 index 0000000..62f206f --- /dev/null +++ b/src/ctedispcube/ctedispcube.hpp @@ -0,0 +1,91 @@ +/*************************************************************************** + * ctedispcube - Energy dispersion cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctedispcube.hpp + * @brief Energy dispersion cube generation tool definition + * @author Chia-Chun Lu + */ + +#ifndef CTEDISPCUBE_HPP +#define CTEDISPCUBE_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" +#include "ctool.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTEDISPCUBE_NAME "ctedispcube" +#define CTEDISPCUBE_VERSION "1.1.0" + + +/***********************************************************************//** + * @class ctedispcube + * + * @brief Energy dispersion cube generation tool + ***************************************************************************/ +class ctedispcube : public ctool { + +public: + // Constructors and destructors + ctedispcube(void); + explicit ctedispcube(const GObservations& obs); + ctedispcube(int argc, char *argv[]); + ctedispcube(const ctedispcube& app); + virtual ~ctedispcube(void); + + // Operators + ctedispcube& operator=(const ctedispcube& app); + + // Methods + void clear(void); + void run(void); + void save(void); + const GCTACubeEdisp& edispcube(void) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctedispcube& app); + void free_members(void); + void get_parameters(void); + + // User parameters + GFilename m_outcube; //!< Output energy dispersion cube file name + + // Protected members + GObservations m_obs; //!< Observation container + GCTACubeEdisp m_edispcube; //!< Energy dispersion cube +}; + + +/***********************************************************************//** + * @brief Return energy dispersion cube + * + * @return Energy dispersion cube + ***************************************************************************/ +inline +const GCTACubeEdisp& ctedispcube::edispcube(void) const +{ + return (m_edispcube); +} + +#endif /* CTEDISPCUBE_HPP */ diff --git a/src/ctedispcube/ctedispcube.par b/src/ctedispcube/ctedispcube.par new file mode 100644 index 0000000..246bb03 --- /dev/null +++ b/src/ctedispcube/ctedispcube.par @@ -0,0 +1,60 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctedispcube.par # +# Version: 1.1.0 # +# # +# Author: Maria Haupt # +# maria.haupt@desy.de # +# DESY # +# # +# Purpose: Parameter file for the ctedispcube application # +# # +############################################################################# +# +# Observational parameters +#=================== +inobs, f, a, "NONE",,, "Input event list or observation definition XML file" +incube, f, a, "NONE",,, "Input counts cube file to extract energy dispersion cube definition" +caldb, s, a, "prod2",,, "Calibration database" +irf, s, a, "South_0.5h",,, "Instrument response function" +outcube, f, a, "edispcube.fits",,, "Output energy dispersion cube file" + +# +# Energy binning parameters +#========================== +ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, a, 20,,, "Number of energy bins" +ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" + +# +# Spatial binning parameters +#=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" +nxpix, i, a, 10,,, "Size of the X axis in pixels" +nypix, i, a, 10,,, "Size of the Y axis in pixels" +binsz, r, a, 1.0,,, "Image scale (in degrees/pixel)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" + +# +# Angular separation binning parameters +#====================================== +migramax, r, h, 2.0,,, "Upper bound of ratio of recontructed to true photon energy" +migrabins, i, h, 100,,, "Number of migration bins" + +# +# 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" +logfile, f, h, "ctedispcube.log",,, "Log filename" diff --git a/src/ctedispcube/main.cpp b/src/ctedispcube/main.cpp new file mode 100644 index 0000000..0a11700 --- /dev/null +++ b/src/ctedispcube/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctedispcube - Energy dispersion cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctedispcube/main.cpp + * @brief Energy dispersion cube tool main code + * @author Maria Haupt + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctedispcube.hpp" + + +/***********************************************************************//** + * @brief Main entry point of application + * + * @param[in] argc Number of command line arguments. + * @param[in] argv Command line arguments. + * + * This is the main entry point of the ctedispcube application. It allocates a + * ctedispcube object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctedispcube application(argc, argv); + + // Run application + try { + // Execute application + application.execute(); + + // Signal success + rc = 0; + } + catch (std::exception &e) { + + // Extract error message + std::string message = e.what(); + std::string signal = "*** ERROR encounterted in the execution of" + " ctedispcube. 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/src/cterror/cterror.par b/src/cterror/cterror.par index 9e7f3c8..86563ae 100644 --- a/src/cterror/cterror.par +++ b/src/cterror/cterror.par @@ -6,7 +6,7 @@ # ------------------------------------------------------------------------- # # # # File: cterror.par # -# Version: 1.0.0 # +# Version: 1.1.0 # # # # Author: Florent Forest # # florent.forest@isae.fr # @@ -23,6 +23,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" srcname, s, a, "Crab",,, "Source of interest" expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp index b0e987f..e74ea09 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -217,18 +217,6 @@ void ctexpcube::run(void) log << std::endl; } - // Set energy dispersion flag for all CTA observations and save old - // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } - // Write observation(s) into logger if (logTerse()) { log << std::endl; @@ -250,14 +238,6 @@ void ctexpcube::run(void) // Fill exposure m_expcube.fill(m_obs, &log); - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } - // Optionally publish exposure cube if (m_publish) { publish(); @@ -348,7 +328,6 @@ void ctexpcube::init_members(void) { // Initialise members m_outcube.clear(); - m_apply_edisp = false; m_publish = false; // Initialise protected members @@ -369,7 +348,6 @@ void ctexpcube::copy_members(const ctexpcube& app) { // Copy attributes m_outcube = app.m_outcube; - m_apply_edisp = app.m_apply_edisp; m_publish = app.m_publish; // Copy protected members @@ -414,9 +392,6 @@ void ctexpcube::get_parameters(void) } // endif: there was no observation in the container - // Read energy dispersion flag - m_apply_edisp = (*this)["edisp"].boolean(); - // Get the incube filename std::string incube = (*this)["incube"].filename(); diff --git a/src/ctexpcube/ctexpcube.hpp b/src/ctexpcube/ctexpcube.hpp index e0d78f0..beac456 100644 --- a/src/ctexpcube/ctexpcube.hpp +++ b/src/ctexpcube/ctexpcube.hpp @@ -71,7 +71,6 @@ protected: // User parameters GFilename m_outcube; //!< Output exposure cube file - bool m_apply_edisp; //!< Apply energy dispersion? bool m_publish; //!< Publish exposure cube? // Protected members diff --git a/src/ctexpcube/ctexpcube.par b/src/ctexpcube/ctexpcube.par index 987a956..9cbb8f7 100644 --- a/src/ctexpcube/ctexpcube.par +++ b/src/ctexpcube/ctexpcube.par @@ -21,7 +21,6 @@ inobs, f, a, "NONE",,, "Input event list or observation definition XML file" incube, f, a, "NONE",,, "Input counts cube file to extract exposure cube definition" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" outcube, f, a, "expcube.fits",,, "Output exposure cube file" # diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index b352d9a..83da110 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctlike.par # -# Version: 1.0.0 # +# Version: 1.1.0 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -21,6 +21,7 @@ inobs, f, a, "events.fits",,, "Input event list, counts cube or observ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube, f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index acfdb69..0025fe9 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -22,6 +22,7 @@ incube, f, a, "NONE",,, "Input counts cube file to extract model cube definiti inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp index 7684358..6d7b3a7 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -216,18 +216,6 @@ void ctpsfcube::run(void) log << std::endl; } - // Set energy dispersion flag for all CTA observations and save old - // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } - // Write observation(s) into logger if (logTerse()) { log << std::endl; @@ -249,14 +237,6 @@ void ctpsfcube::run(void) // Fill PSF m_psfcube.fill(m_obs, &log); - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } - // Return return; } @@ -308,7 +288,6 @@ void ctpsfcube::init_members(void) { // Initialise members m_outcube.clear(); - m_apply_edisp = false; // Initialise protected members m_obs.clear(); @@ -328,7 +307,6 @@ void ctpsfcube::copy_members(const ctpsfcube& app) { // Copy attributes m_outcube = app.m_outcube; - m_apply_edisp = app.m_apply_edisp; // Copy protected members m_obs = app.m_obs; @@ -401,9 +379,6 @@ void ctpsfcube::get_parameters(void) } // endelse: cube loaded from file - // Read energy dispersion flag - m_apply_edisp = (*this)["edisp"].boolean(); - // Read output filename (if needed) if (read_ahead()) { m_outcube = (*this)["outcube"].filename(); diff --git a/src/ctpsfcube/ctpsfcube.hpp b/src/ctpsfcube/ctpsfcube.hpp index 55d9e49..7b469b3 100644 --- a/src/ctpsfcube/ctpsfcube.hpp +++ b/src/ctpsfcube/ctpsfcube.hpp @@ -70,7 +70,6 @@ protected: // User parameters GFilename m_outcube; //!< Output PSF cube file name - bool m_apply_edisp; //!< Apply energy dispersion? // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctpsfcube/ctpsfcube.par b/src/ctpsfcube/ctpsfcube.par index ab8961c..b9f99eb 100644 --- a/src/ctpsfcube/ctpsfcube.par +++ b/src/ctpsfcube/ctpsfcube.par @@ -21,7 +21,6 @@ inobs, f, a, "NONE",,, "Input event list or observation definition XML file" incube, f, a, "NONE",,, "Input counts cube file to extract PSF cube definition" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" outcube, f, a, "psfcube.fits",,, "Output PSF cube file" # diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 89b784c..6c23137 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -17,35 +17,36 @@ # # General parameters #=================== -inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file" -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" -srcname, s, a, "Crab",,, "Test source name" -expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" -psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" -bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" -caldb, s, a, "prod2",,, "Calibration database" -irf, s, a, "South_0.5h",,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" -outmap, f, a, "tsmap.fits",,, "Output Test Statistic map file" +inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" +srcname, s, a, "Crab",,, "Test source name" +expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" +psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube, f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" +caldb, s, a, "prod2",,, "Calibration database" +irf, s, a, "South_0.5h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" +outmap, f, a, "tsmap.fits",,, "Output Test Statistic map file" # # Spatial binning parameters #=========================== -usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" -nxpix, i, a, 200,,, "Size of the X axis in pixels" -nypix, i, a, 200,,, "Size of the Y axis in pixels" -binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" -xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" -yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" -proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" +nxpix, i, a, 200,,, "Size of the X axis in pixels" +nypix, i, a, 200,,, "Size of the Y axis in pixels" +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" # # Parameters for splitting and speed purpose #=========================== -binmin, i, h, -1,,, "First bin to compute" -binmax, i, h, -1,,, "Last bin to compute" -logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis" +binmin, i, h, -1,,, "First bin to compute" +binmax, i, h, -1,,, "Last bin to compute" +logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis" # # Standard parameters diff --git a/src/ctulimit/ctulimit.par b/src/ctulimit/ctulimit.par index 3e4b3b6..dff2582 100644 --- a/src/ctulimit/ctulimit.par +++ b/src/ctulimit/ctulimit.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctulimit.par # -# Version: 01-00-00 # +# Version: 1.1.0 # # # # Author: Michael Mayer # # michael.mayer@physik.hu-berlin.de # @@ -22,6 +22,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" srcname, s, a, "Crab",,, "Source of interest" expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)" +edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)" bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)" caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index ee6f9fe..5cc444d 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -862,15 +862,18 @@ void ctool::set_obs_response(GCTAObservation* obs) // If the observation contains a counts cube, then first check whether // the expcube and psfcube parameters exist and are not NONE if (dynamic_cast(obs->events()) != NULL) { - if (has_par("expcube") && has_par("psfcube") && has_par("bkgcube")) { + if (has_par("expcube") && has_par("psfcube") && + has_par("edispcube") && has_par("bkgcube")) { // Get filenames - std::string expcube = (*this)["expcube"].filename(); - std::string psfcube = (*this)["psfcube"].filename(); - std::string bkgcube = (*this)["bkgcube"].filename(); + std::string expcube = (*this)["expcube"].filename(); + std::string psfcube = (*this)["psfcube"].filename(); + std::string edispcube = (*this)["edispcube"].filename(); + std::string bkgcube = (*this)["bkgcube"].filename(); // Extract response information if available - if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") && + if ((expcube != "NONE") && (psfcube != "NONE") && + (bkgcube != "NONE") && (gammalib::strip_whitespace(expcube) != "") && (gammalib::strip_whitespace(psfcube) != "") && (gammalib::strip_whitespace(bkgcube) != "")) { @@ -880,8 +883,16 @@ void ctool::set_obs_response(GCTAObservation* obs) GCTACubePsf psf(psfcube); GCTACubeBackground background(bkgcube); - // Set reponse - obs->response(exposure, psf, background); + // Optionally load energy dispersion cube + if ((edispcube != "NONE") && + (gammalib::strip_whitespace(edispcube) != "")) { + GCTACubeEdisp edisp(edispcube); + obs->response(exposure, psf, edisp, background); + } + else { + // Set reponse + obs->response(exposure, psf, background); + } // Signal response availability has_response = true; @@ -904,8 +915,9 @@ void ctool::set_obs_response(GCTAObservation* obs) // to the observation definition XML file, in case an observation // definition XML file is written. std::string parameter = "parameter name=\"Calibration\"" - " database=\""+database+"\"" - " response=\""+irf+"\""; + " database=\""+database+"\"" + " response=\""+irf+"\""; + GXmlElement xml; xml.append(parameter); @@ -915,22 +927,6 @@ void ctool::set_obs_response(GCTAObservation* obs) // Attach response to observation obs->response(response); - // Set calibration database. If "database" is a valid directory then use - // this as the pathname to the calibration database. Otherwise, interpret - // "database" as the instrument name, the mission being "cta" - /* - GCaldb caldb; - if (gammalib::dir_exists(database)) { - caldb.rootdir(database); - } - else { - caldb.open("cta", database); - } - */ - - // Set reponse - //obs->response(irf, caldb); - // Signal response availability has_response = true; @@ -967,11 +963,13 @@ void ctool::set_obs_bounds(GObservations& obs) // parameters and add them if (list->ebounds().is_empty()) { - // If there are "emin" and "emax" parameters then use them + // If there are "emin" and "emax" parameters then use + // them if (has_par("emin") && has_par("emax")) { double emin = (*this)["emin"].real(); double emax = (*this)["emax"].real(); - GEbounds ebounds(GEnergy(emin, "TeV"), GEnergy(emax, "TeV")); + GEbounds ebounds(GEnergy(emin, "TeV"), + GEnergy(emax, "TeV")); list->ebounds(ebounds); } @@ -996,7 +994,6 @@ void ctool::set_obs_bounds(GObservations& obs) } // endfor: looped over observations - // Return return; } @@ -1126,7 +1123,7 @@ void ctool::provide_help(void) const char* ptr = std::getenv("CTOOLS"); if (ptr == NULL) { std::string msg = "CTOOLS environment variable not set, cannot " - "display help file. Please define CTOOLS " + "display help file. Please set the CTOOLS " "environment variable."; throw GException::invalid_value(G_PROVIDE_HELP, msg); } diff --git a/test/data/crab_edispcube.fits.gz b/test/data/crab_edispcube.fits.gz new file mode 100644 index 0000000..370dbe5 --- /dev/null +++ b/test/data/crab_edispcube.fits.gz diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index e367823..673d345 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -405,6 +405,7 @@ fi inobs="data/crab_cntmap_small.fits.gz" \ expcube="data/crab_expcube.fits.gz" \ psfcube="data/crab_psfcube.fits.gz" \ + edispcube="data/crab_edispcube.fits.gz" \ bkgcube="data/crab_bkgcube.fits.gz" \ inmodel="data/crab_bkgcube.xml" \ ebinalg="LOG" \ diff --git a/test/test_ctedispcube.py b/test/test_ctedispcube.py new file mode 100755 index 0000000..e91afea --- /dev/null +++ b/test/test_ctedispcube.py @@ -0,0 +1,111 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the ctedispcube tool. +# +# Copyright (C) 2016 Maria Haupt +# +# 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 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools + + +# ============================= # +# Test class for ctedispcube tool # +# ============================= # +class Test(gammalib.GPythonTestSuite): + """ + Test class for ctedispcube tool. + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Set members + self.events_name = "data/crab_events.fits" + self.caldb = "irf" + self.irf = "cta_dummy_irf" + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("ctedispcube") + + # Append tests + self.append(self.test_functional, "Test ctedispcube functionality") + + # Return + return + + # Test ctedispcube functionnality + def test_functional(self): + """ + Test ctedispcube functionnality. + """ + # Set-up ctedispcube + edispcube = ctools.ctedispcube() + edispcube["inobs"] = self.events_name + edispcube["incube"] = "NONE" + edispcube["outcube"] = "edispcube.fits" + edispcube["caldb"] = self.caldb + edispcube["irf"] = self.irf + edispcube["ebinalg"] = "LOG" + edispcube["emin"] = 0.1 + edispcube["emax"] = 100 + edispcube["enumbins"] = 20 + edispcube["nxpix"] = 10 + edispcube["nypix"] = 10 + edispcube["binsz"] = 0.4 + edispcube["coordsys"] = "CEL" + edispcube["proj"] = "CAR" + edispcube["xref"] = 83.63 + edispcube["yref"] = 22.01 + edispcube["migramax"] = 2.0 + edispcube["migrabins"] = 10 + + # Run tool + self.test_try("Run ctedispcube") + try: + edispcube.run() + self.test_try_success() + #except Exception as e: + # msg = "Exception occured in ctedispcube: %s." % (e,) + except: + msg = "Exception occured in ctedispcube." + self.test_try_failure(msg) + + # Save EDISP cube + self.test_try("Save EDISP cube") + try: + edispcube.save() + self.test_try_success() + #except Exception as e: + # msg = "Exception occured in saving EDISP cube: %s." % (e,) + except: + msg = "Exception occured in saving EDISP cubes." + self.test_try_failure(msg) + + # Return + return diff --git a/test/test_ctools.sh b/test/test_ctools.sh index a730162..5b292eb 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -2,7 +2,7 @@ # ========================================================================== # This script tests all ctools # -# Copyright (C) 2011-2015 Juergen Knoedlseder +# Copyright (C) 2011-2016 Juergen Knoedlseder # # 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 @@ -43,6 +43,7 @@ ctselect=../src/ctselect/ctselect ctskymap=../src/ctskymap/ctskymap ctexpcube=../src/ctexpcube/ctexpcube ctpsfcube=../src/ctpsfcube/ctpsfcube +ctedispcube=../src/ctedispcube/ctedispcube ctbkgcube=../src/ctbkgcube/ctbkgcube ctcubemask=../src/ctcubemask/ctcubemask ctbutterfly=../src/ctbutterfly/ctbutterfly @@ -272,6 +273,57 @@ $ECHO " ok" # +# Test ctedispcube +# ============== +$ECHO -n "Test ctedispcube: " +# +# Run 1 +$ctedispcube inobs="data/crab_events.fits" \ + incube="NONE" \ + outcube="edispcube1.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=10 \ + nypix=10 \ + binsz=0.4 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" \ + migramax=3.0 \ + migrabins=10 +$ECHO -n "." +if [ -s "edispcube1.fits" ] +then +$ECHO -n "." +else +$ECHO " edispcube1.fits file is not found" +exit 1 +fi +# +# Run 2 +$ctedispcube inobs="events.fits" \ + incube="cntmap2.fits" \ + outcube="edispcube2.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + migramax=3.0 \ + migrabins=3 +$ECHO -n "." +if [ -s "edispcube2.fits" ] +then +$ECHO -n "." +else +$ECHO " edispcube2.fits file is not found" +exit 1 +fi +$ECHO " ok" + + +# # Test ctbkgcube # ============== $ECHO -n "Test ctbkgcube: " @@ -336,7 +388,6 @@ fi $ECHO " ok" - # # Test ctmodel # ============ @@ -349,6 +400,7 @@ $ctmodel inobs="cntmap1.fits" \ expcube="NONE" \ psfcube="NONE" \ bkgcube="NONE" \ + edispcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" \ inmodel="data/crab.xml" @@ -368,6 +420,7 @@ $ctmodel inobs="cntmap2.fits" \ outcube="modmap2.fits" \ expcube="expcube2.fits" \ psfcube="psfcube2.fits" \ + edispcube="edispcube2.fits" \ bkgcube="bkgcube2.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ @@ -388,6 +441,7 @@ $ctmodel inobs="NONE" \ outcube="modmap3.fits" \ expcube="NONE" \ psfcube="NONE" \ + edispcube="NONE" \ bkgcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" \ @@ -453,6 +507,7 @@ $ctlike inobs="cntmap1.fits" \ outmodel="results_binned.xml" \ expcube="NONE" \ psfcube="NONE" \ + edispcube="NONE" \ bkgcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" @@ -487,6 +542,7 @@ $ctlike inobs="cntmap2.fits" \ outmodel="results_binned_cube_background.xml" \ expcube="expcube2.fits" \ psfcube="psfcube2.fits" \ + edispcube="NONE" \ bkgcube="bkgcube2.fits" \ caldb="irf" \ irf="cta_dummy_irf" @@ -498,6 +554,28 @@ else $ECHO " results_binned_cube_background.xml file is not found" exit 1 fi + + +# +# Run 4 +$ctlike inobs="cntmap2.fits" \ + inmodel="bkgcube2.xml" \ + outmodel="results_binned_cube_background_edisp.xml" \ + expcube="expcube2.fits" \ + psfcube="psfcube2.fits" \ + edispcube="edispcube2.fits" \ + bkgcube="bkgcube2.fits" \ + edisp="yes" + caldb="irf" \ + irf="cta_dummy_irf" +$ECHO -n "." +if [ -s "results_binned_cube_background_edisp.xml" ] +then + $ECHO -n "." +else + $ECHO " results_binned_cube_background_edisp.xml file is not found" + exit 1 +fi $ECHO " ok" @@ -529,8 +607,6 @@ fi $ECHO " ok" - - # # Test ctcubemask # ===============