diff --git a/ChangeLog b/ChangeLog index ab12784..60f5aee 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,18 +1,19 @@ -2014-10-10 Juergen Knoedlseder +2014-10-15 Juergen Knoedlseder * Version 0.8.0 released ======================== Clean-up public interface - Adapt to GammaLib 0.9.0 interface - Introduce ctool base class - ctbin now always merges event lists into a single counts cube - ctmodel now always creates a single model cube + Adapt to GammaLib 0.9.0 interface Add ctexpcube tool (generates exposure cube) Add ctpsfcube tool (generates PSF cube) Add ctbkgcube tool (generates background cube) Add ctcubemask tool (masks event bins in event cube) Add cttsmap tool (generates TS map) + Add csresmap script (generates residual map) + Introduce ctool base class + ctbin now always merges event lists into a single counts cube + ctmodel now always creates a single model cube Support energy dispersion in ctobssim, ctlike and ctmodel Correct ctobssim for thread safety Correctly apply time references in ctselect diff --git a/NEWS b/NEWS index 4b576b1..e08cb1c 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -10 October 2014 +15 October 2014 1. Introduction @@ -192,6 +192,9 @@ The following functions have been added to obsutils: * spectrum() generates spectral points by fitting a source model in specified energy bins. +Add the csresmap script that uses ctbin and ctmodel to generate a +residual counts cube for a given model. + 18. caldb --------- diff --git a/pyext/ctmodel.i b/pyext/ctmodel.i index dd43f3d..27632f3 100644 --- a/pyext/ctmodel.i +++ b/pyext/ctmodel.i @@ -49,6 +49,7 @@ public: void save(void); const GObservations& obs(void) const; const GCTAEventCube& cube(void) const; + void cube(const GCTAEventCube& cube); }; diff --git a/scripts/csresmap.py b/scripts/csresmap.py index 70fbbcd..104d744 100755 --- a/scripts/csresmap.py +++ b/scripts/csresmap.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This script generates a residual map. +# Residual map generation script. # # Copyright (C) 2014 Michael Mayer # @@ -18,18 +18,17 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * +import ctools +import gammalib import sys -import csv -import math import tempfile +import os # ============== # # cstsdist class # # ============== # -class csresmap(GApplication): +class csresmap(gammalib.GApplication): """ This class implements the creation of a residual map. It derives from the GammaLib::GApplication class which provides support for parameter @@ -45,14 +44,17 @@ class csresmap(GApplication): self.version = "0.1.0" # Initialise some members - self.obs = None - self.algorithm = "SUB" + self.obs = None + self.algorithm = "SUB" + self.resmap = None # Initialise some members - if isinstance(argv[0],GObservations): + if isinstance(argv[0],gammalib.GObservations): self.obs = argv[0] + argv = argv[1:] else: - self.obs = GObservations() + self.obs = gammalib.GObservations() + self.obs.clear() self.outfile = "" # Make sure that parfile exists @@ -60,9 +62,9 @@ class csresmap(GApplication): # Initialise application if len(argv) == 0: - GApplication.__init__(self, self.name, self.version) + gammalib.GApplication.__init__(self, self.name, self.version) elif len(argv) ==1: - GApplication.__init__(self, self.name, self.version, argv[1:]) + gammalib.GApplication.__init__(self, self.name, self.version, *argv) else: raise TypeError("Invalid number of arguments given.") @@ -93,32 +95,32 @@ class csresmap(GApplication): parfile = self.name+".par" try: - pars = GApplicationPars(parfile) + pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") # Create default parfile - pars = GApplicationPars() - pars.append(GApplicationPar("infile","f","h","obs.xml","","","Observation definition")) - pars.append(GApplicationPar("srcmdl","f","h","crab.xml","","","model file name")) - pars.append(GApplicationPar("outfile","f","h","resmap.fits","","","Output file name")) - pars.append(GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name")) - pars.append(GApplicationPar("modfile","f","h","model.fits","","","Output model file name")) - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)")) - pars.append(GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)")) - pars.append(GApplicationPar("enumbins","i","h","20","","","Number of energy bins")) - pars.append(GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm")) - pars.append(GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System")) - pars.append(GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN")) - pars.append(GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)")) - pars.append(GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)")) - pars.append(GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels")) - pars.append(GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels")) - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)")) - pars.append(GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm")) + pars = gammalib.GApplicationPars() + pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition")) + pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name")) + pars.append(gammalib.GApplicationPar("outfile","f","a","resmap.fits","","","Output file name")) + pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name")) + pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name")) + pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)")) + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)")) + pars.append(gammalib.GApplicationPar("enumbins","i","h","20","","","Number of energy bins")) + pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm")) + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System")) + pars.append(gammalib.GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN")) + pars.append(gammalib.GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)")) + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)")) + pars.append(gammalib.GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels")) + pars.append(gammalib.GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels")) + pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)")) + pars.append(gammalib.GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm")) pars.append_standard() pars.save(parfile) @@ -133,9 +135,27 @@ class csresmap(GApplication): # Set observation if not done before if self.obs.size() == 0: - self.obs = GObservations(self["infile"].filename) - self.m_irf = self["irf"].string() - self.m_caldb = self["caldb"].string() + + obsfile = self["infile"].filename() + try: + + self.obs = gammalib.GObservations(obsfile) + except: + self.obs.clear() + observation = gammalib.GCTAObservation() + observation.load(obsfile) + + self.m_irf = self["irf"].string() + self.m_caldb = self["caldb"].string() + + caldb = gammalib.GCaldb() + if os.path.isdir(self.m_caldb): + caldb.rootdir(self.m_caldb) + else: + caldb.open("cta",self. m_caldb) + + observation.response(self.m_irf, caldb); + self.obs.append(observation) # Check for models in the container # read models from file if there were none @@ -185,6 +205,9 @@ class csresmap(GApplication): """ # Run the script self.run() + + # Save residual map + self.resmap.save(self.m_outfile, self.m_clobber) # Return return @@ -218,7 +241,8 @@ class csresmap(GApplication): self.log.header1("Generate binned map") # Create countsmap - bin = ctbin(self.obs) + bin = ctools.ctbin(self.obs) + bin["debug"].boolean(True) bin["nxpix"].integer(self.m_nxpix) bin["nypix"].integer(self.m_nypix) bin["proj"].string(self.m_proj) @@ -231,75 +255,59 @@ class csresmap(GApplication): bin["emax"].real(self.m_emax) bin["binsz"].real(self.m_binsz) bin["clobber"].boolean(self.m_clobber) - if not self.m_binfile == "NONE": - bin["outfile"].filename(self.m_binfile) - else: - os,filename = tempfile.mkstemp() - bin["outfile"].filename(filename) - bin.execute() + bin.run() + + # 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 + # calling bin.cube() + #residualmap = bin.cube().map() + counts_cube = bin.cube() + self.resmap = counts_cube.map().copy() + self.resmap.stack_maps() # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate model map") - # create model map - model = ctmodel(self.obs) - model["infile"].filename(bin["outfile"].filename()) + # Create model map + model = ctools.ctmodel(self.obs) + model["debug"].boolean(True) + model.cube(bin.cube()) model["clobber"].boolean(self.m_clobber) - if not self.m_modfile == "NONE": - model["outfile"].filename(self.m_modfile) - else: - os,filename = tempfile.mkstemp() - model["outfile"].filename(filename) - model.execute() + model.run() # Get model map into GSkymap object - modelmap = GSkymap(model["outfile"].filename()) + modelmap = model.cube().map() modelmap.stack_maps() - - # Initialise resiudalmap from countmap - residualmap = GSkymap(bin["outfile"].filename()) - residualmap.stack_maps() - - - - # continue calculations depending on given algorithm + # Continue calculations depending on given algorithm if self.m_algorithm == "SUB": # Subtract maps - residualmap -= modelmap + self.resmap -= modelmap elif self.m_algorithm == "SUBDIV": # Subtract and divide by model map - residualmap -= modelmap - residualmap /= modelmap + self.resmap -= modelmap + self.resmap /= modelmap elif self.m_algorithm == "SUBDIVSQRT": # subtract and divide by sqrt of model map - residualmap -= modelmap - residualmap /= modelmap.sqrt() + self.resmap -= modelmap + self.resmap /= modelmap.sqrt() else: - # Rais error if algorithm is unkown + # Raise error if algorithm is unkown raise TypeError("Algorithm \""+self.m_algorithm+"\" not known") - # Save files - residualmap.save(self.m_outfile,self.m_clobber) - - # Remove binned and model files in case they shouldnt be saved - if self.m_binfile == "NONE": - os.remove(bin["outfile"].filename()) - if self.m_modfile == "NONE": - os.remove(model["outfile"].filename()) - # Return return - + # ======================== # # Main routine entry point # diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 43e5019..ac273d2 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -246,8 +246,10 @@ void ctmodel::run(void) log.header1("Generate model cube"); } - // Initialise counts cube - init_cube(); + // Initialise counts cube if cube has not yet been set + if (!has_cube()) { + init_cube(); + } // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { @@ -320,6 +322,31 @@ void ctmodel::save(void) } +/***********************************************************************//** + * @brief Set model cube + * + * @param[in] cube Model cube. + * + * Set model cube and set all cube bins to zero. + ***************************************************************************/ +void ctmodel::cube(const GCTAEventCube& cube) +{ + // Set cube + m_cube = cube; + + // Set all cube bins to zero + for (int i = 0; i < m_cube.size(); ++i) { + m_cube[i]->counts(0.0); + } + + // Signal that cube has been set + m_has_cube = true; + + // Return + return; +} + + /*========================================================================== = = = Private methods = @@ -362,6 +389,7 @@ void ctmodel::init_members(void) m_cube.clear(); m_ebounds.clear(); m_gti.clear(); + m_has_cube = false; // Return return; @@ -406,6 +434,7 @@ void ctmodel::copy_members(const ctmodel& app) m_cube = app.m_cube; m_ebounds = app.m_ebounds; m_gti = app.m_gti; + m_has_cube = app.m_has_cube; // Return return; @@ -431,7 +460,9 @@ void ctmodel::free_members(void) void ctmodel::get_parameters(void) { // Read input and optionally output cube filenames - m_infile = (*this)["infile"].filename(); + if (!has_cube()) { + m_infile = (*this)["infile"].filename(); + } if (read_ahead()) { m_outfile = (*this)["outfile"].filename(); } @@ -482,7 +513,8 @@ void ctmodel::get_parameters(void) // If no cube file has been specified then read all parameters that // are necessary to create the cube from scratch (see method init_cube) - if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { + if ((!has_cube()) && + ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == ""))) { m_ebinalg = (*this)["ebinalg"].string(); if (m_ebinalg == "FILE") { m_ebinfile = (*this)["ebinfile"].filename(); @@ -670,6 +702,9 @@ void ctmodel::init_cube(void) } + // Signal that cube has been set + m_has_cube = true; + // Return return; } @@ -738,7 +773,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs) if (logTerse()) { log << gammalib::parformat("Model events in cube"); log << sum << std::endl; - log << gammalib::parformat("Cube bins outside energy range"); + log << gammalib::parformat("Bins outside energy range"); log << num_outside_ebds << std::endl; } diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 52f71c7..7c43083 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -70,6 +70,7 @@ public: void save(void); const GObservations& obs(void) const; const GCTAEventCube& cube(void) const; + void cube(const GCTAEventCube& cube); protected: // Protected methods @@ -80,6 +81,7 @@ protected: void setup_obs(void); void init_cube(void); void fill_cube(const GCTAObservation* obs); + bool has_cube(void) const; // User parameters std::string m_infile; //!< Input counts cube @@ -112,13 +114,14 @@ protected: GCTAEventCube m_cube; //!< Model cube GEbounds m_ebounds; //!< Energy boundaries GGti m_gti; //!< Model cube GTIs + bool m_has_cube; //!< Signal if cube has been set or loaded }; /***********************************************************************//** * @brief Return observation container * - * @return Reference to observation container + * @return Reference to observation container. ***************************************************************************/ inline const GObservations& ctmodel::obs(void) const @@ -130,7 +133,7 @@ const GObservations& ctmodel::obs(void) const /***********************************************************************//** * @brief Return model cube * - * @return Reference to model cube + * @return Reference to model cube. ***************************************************************************/ inline const GCTAEventCube& ctmodel::cube(void) const @@ -138,4 +141,16 @@ const GCTAEventCube& ctmodel::cube(void) const return m_cube; } + +/***********************************************************************//** + * @brief Signal if cube has been set or loaded + * + * @return True if cube has been set or loaded. + ***************************************************************************/ +inline +bool ctmodel::has_cube(void) const +{ + return m_has_cube; +} + #endif /* CTMODEL_HPP */ diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 643bbea..9b257fc 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -131,3 +131,32 @@ else exit 1 fi $ECHO " ok" + +# +# Test csresmap +# ============= +$ECHO @ECHO_N@ "Test csresmap: " +chmod +x @abs_top_srcdir@/scripts/csresmap.py +@abs_top_srcdir@/scripts/csresmap.py \ + infile="data/crab_events.fits" \ + outfile="resmap.fits" \ + srcmdl="data/crab.xml" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + nxpix=50 \ + nypix=50 \ + binsz=0.02 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" \ + binsz=0.05 > csresmap.log +$ECHO @ECHO_N@ "." +if [ -s "resmap.fits" ] +then + $ECHO @ECHO_N@ "." +else + $ECHO " resmap.fits file is not valid" + exit 1 +fi +$ECHO " ok"