diff --git a/.gitignore b/.gitignore index 0c9cb1f..7a0b31c 100644 --- a/.gitignore +++ b/.gitignore @@ -22,9 +22,13 @@ ltmain.sh m4 missing .libs +.deps stamp-h1 +depcomp +src/ctools-setup pyext/build/ +pyext/setup.py pyext/ctools.py pyext/ctools_wrap.cpp src/ctbin/ctbin @@ -34,3 +38,7 @@ src/ctobssim/ctobssim src/ctselect/ctselect src/ctskymap/ctskymap test/pfiles/ +test/*.fits +test/*.dat +test/*.xml +test/test_cscripts.sh diff --git a/AUTHORS b/AUTHORS index 6815846..cef2638 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,3 +1,5 @@ Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) - overall design Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) - openmp support +Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) - HESS support +Mathieu De Naurois (denauroi@in2p3.fr) - HESS support diff --git a/ChangeLog b/ChangeLog index a22cc71..ed8d1bc 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2012-12-05 Juergen Knoedlseder + + * Version 0.6.0 released + ======================== + + Add ctmodel tool + Adapt to GammaLib-00-07-00 interface + 2012-10-11 Juergen Knoedlseder * Bug fix version 0.5.2 released diff --git a/configure.ac b/configure.ac index e2a5eec..3e82b4a 100644 --- a/configure.ac +++ b/configure.ac @@ -32,7 +32,8 @@ AC_INIT([ctools], [0.6.0], [jurgen.knodlseder@irap.omp.eu], [ctools]) AC_CONFIG_SRCDIR([src/ctselect/ctselect.cpp]) AC_CONFIG_MACRO_DIR([m4]) -AM_INIT_AUTOMAKE([gnu no-dependencies]) +#AM_INIT_AUTOMAKE([gnu no-dependencies]) +AM_INIT_AUTOMAKE([gnu]) AM_CONFIG_HEADER([config.h]) @@ -870,11 +871,12 @@ AC_SUBST(LIB_LINK_OPTION) ############################################################################# AC_CONFIG_FILES([Makefile src/Makefile - src/ctselect/Makefile src/ctbin/Makefile src/ctlike/Makefile + src/ctmodel/Makefile src/ctobssim/Makefile src/ctskymap/Makefile + src/ctselect/Makefile src/ctools-setup test/Makefile test/test_cscripts.sh diff --git a/pyext/ctmodel.i b/pyext/ctmodel.i new file mode 100644 index 0000000..c530c29 --- /dev/null +++ b/pyext/ctmodel.i @@ -0,0 +1,67 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by 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 * + * 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 ctmodel.i + * @brief CTA counts model tool definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctmodel.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class ctmodel + * + * @brief CTA counts model tool interface defintion + ***************************************************************************/ +class ctmodel : public GApplication { + +public: + // Constructors and destructors + ctmodel(void); + explicit ctmodel(GObservations obs); + ctmodel(int argc, char *argv[]); + ctmodel(const ctmodel& app); + virtual ~ctmodel(void); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + GObservations& obs(void) { return m_obs; } + void get_parameters(void); + void setup_obs(void); + void model_map(GCTAObservation* obs, const GModels& models); +}; + + +/***********************************************************************//** + * @brief CTA counts model tool Python extension + ***************************************************************************/ +%extend ctmodel { + ctmodel copy() { + return (*self); + } +} diff --git a/pyext/ctools.i b/pyext/ctools.i index 2f145e4..4880856 100755 --- a/pyext/ctools.i +++ b/pyext/ctools.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctools - SWIG file * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2011 by Jurgen Knodlseder * + * copyright (C) 2010-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,6 +21,11 @@ * Usage: * * swig -c++ -python -Wall ctools.i * ***************************************************************************/ +/** + * @file ctools.i + * @brief ctools SWIG file + * @author Juergen Knoedlseder + */ %module ctools %feature("autodoc", "1"); @@ -31,10 +36,12 @@ %import(module="gammalib.support") "GException.i"; /* __ Inform about base classes __________________________________________ */ -%import(module="gammalib.app") "GApplication.i"; +%import(module="gammalib.base") "GBase.i"; +%import(module="gammalib.app") "GApplication.i"; -/* __ CTA tools __________________________________________________________ */ -%include "ctobssim.i" -%include "ctselect.i" +/* __ ctools _____________________________________________________________ */ %include "ctbin.i" +%include "ctobssim.i" %include "ctlike.i" +%include "ctmodel.i" +%include "ctselect.i" diff --git a/scripts/cta_root2caldb.py b/scripts/cta_root2caldb.py index b757081..1d63fb5 100755 --- a/scripts/cta_root2caldb.py +++ b/scripts/cta_root2caldb.py @@ -2,11 +2,9 @@ # ========================================================================== # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D # performance files. -# -# Requirements: -# - pyroot -# -# Copyright (C) 2011 Jurgen Knodlseder +# ------------------------------------------------------------------------- +# Copyright (C) 2011-2012 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 @@ -21,6 +19,10 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # +# ------------------------------------------------------------------------- +# Requirements: +# - gammalib +# - pyroot # ========================================================================== from ROOT import TFile, TH1F, TH2F from gammalib import * @@ -37,18 +39,27 @@ import os class caldb(): """ """ - def __init__(self, conf, path=""): + def __init__(self, inst, obsid, path=""): """ - Constructor. + CALDB constructor. The constructor will create one CALDB entry for + an observation. + + Parameters: + inst - Instrument string. + obsid - Observation ID string. + Keywords: + path - CALDB Root path. """ # Store root path self.path = path # Initialise some members self.cif = None + self.irf = None self.ea = None self.psf = None self.edisp = None + self.bgd = None self.hdu_cif = None # Store UTC date of job execution @@ -56,7 +67,8 @@ class caldb(): # Set CALDB information self.cal_tel = "CTA" - self.cal_inst = conf.upper() + self.cal_inst = inst.upper() + self.cal_obsid = obsid self.cal_det = "NONE" self.cal_flt = "NONE" self.cal_class = "BCF" @@ -84,6 +96,11 @@ class caldb(): self.edisp_doc = "???" self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] self.edisp_desc = "CTA energy dispersion" + # + self.bgd_name = "BGD" + self.bgd_doc = "???" + self.bgd_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_desc = "CTA background" # Create directory structure self.make_dirs() @@ -93,6 +110,7 @@ class caldb(): def __del__(self): """ + Destructor. Close all FITS files. """ # Close calibration self.close() @@ -102,32 +120,50 @@ class caldb(): def make_dirs(self): """ + Generate CALDB directory structure for one observation identifier. + The structure is given by + + data///bcf/ + + where is "cta" and is the instrument specified in the + CALDB constructor (the instrument may be used for different array + configurations). + + Parameters: + None + Keywords: + None """ # Set base path self.base_dir = "data" self.base_dir += "/"+self.cal_tel.lower() self.base_dir += "/"+self.cal_inst.lower() - # Set directory names self.bcf_dir = self.base_dir+"/bcf" - self.ea_dir = self.bcf_dir+"/ea" - self.psf_dir = self.bcf_dir+"/psf" - self.edisp_dir = self.bcf_dir+"/edisp" + self.obs_dir = self.bcf_dir+"/"+self.cal_obsid + self.ea_dir = self.obs_dir + self.psf_dir = self.obs_dir + self.edisp_dir = self.obs_dir + self.bgd_dir = self.obs_dir - # Add path (if required) + # Optionally prefix path if len(self.path) > 0: self.base_path = self.path+"/"+self.base_dir self.bcf_path = self.path+"/"+self.bcf_dir + self.obs_path = self.path+"/"+self.obs_dir self.ea_path = self.path+"/"+self.ea_dir self.psf_path = self.path+"/"+self.psf_dir self.edisp_path = self.path+"/"+self.edisp_dir + self.bgd_path = self.path+"/"+self.bgd_dir else: self.base_path = self.base_dir self.bcf_path = self.bcf_dir + self.obs_path = self.obs_dir self.ea_path = self.ea_dir self.psf_path = self.psf_dir self.edisp_path = self.edisp_dir + self.bgd_path = self.bgd_dir # Create directory structure if not os.path.isdir(self.ea_path): @@ -136,18 +172,36 @@ class caldb(): os.makedirs(self.psf_path) if not os.path.isdir(self.edisp_path): os.makedirs(self.edisp_path) + if not os.path.isdir(self.bgd_path): + os.makedirs(self.bgd_path) # Return return - def open(self, name, clobber=True): + def open(self, version, split=False, clobber=True): """ - Open existing or create new calibration. + Open existing or create new calibration. The actual version will put + all calibrations in the same file, although each part of the response + function will have its own logical name. We can thus easily modify + the script to put each calibration information in a separate file. + + Parameters: + version - Filename version + Keywords: + split - Split IRF over several files? + clobber - Overwrite existing files? """ # Set calibrate file names - self.ea_file = "ea_"+name+".fits" - self.psf_file = "psf_"+name+".fits" - self.edisp_file = "edisp_"+name+".fits" + if split: + self.ea_file = "ea_"+version+".fits" + self.psf_file = "psf_"+version+".fits" + self.edisp_file = "edisp_"+version+".fits" + self.bgd_file = "bgd_"+version+".fits" + else: + self.ea_file = "irf_"+version+".fits" + self.psf_file = "irf_"+version+".fits" + self.edisp_file = "irf_"+version+".fits" + self.bgd_file = "irf_"+version+".fits" # Open calibration database index self.cif = GFits(self.base_path+"/caldb.indx", True) @@ -159,17 +213,42 @@ class caldb(): self.create_cif() self.hdu_cif = self.cif.table("CIF") - # Open calibration files - self.open_ea(self.ea_path+"/"+self.ea_file) - self.open_psf(self.psf_path+"/"+self.psf_file) - self.open_edisp(self.edisp_path+"/"+self.edisp_file) + # Set filenames + ea_filename = self.ea_path+"/"+self.ea_file + psf_filename = self.psf_path+"/"+self.psf_file + edisp_filename = self.edisp_path+"/"+self.edisp_file + bgd_filename = self.bgd_path+"/"+self.bgd_file + + # Open files + if split: + self.ea = GFits(ea_filename, True) + self.psf = GFits(psf_filename, True) + self.edisp = GFits(edisp_filename, True) + self.bgd = GFits(bgd_filename, True) + else: + self.irf = GFits(ea_filename, True) + self.ea = self.irf + self.psf = self.irf + self.edisp = self.irf + self.bgd = self.irf + + # Open HDUs + self.open_ea() + self.open_psf() + self.open_edisp() + self.open_bgd() # Return return def close(self): """ - Close calibration. + Close calibration FITS files. + + Parameters: + None + Keywords: + None """ # Add information to CIF. We do this now as before we did not # necessarily have all the information at hand (in particular @@ -181,31 +260,57 @@ class caldb(): self.cif.save(True) self.cif.close() self.cif = None - - # Close effective area - if self.ea != None: - self.ea.save(True) - self.ea.close() - self.ea = None - - # Close point spread function - if self.psf != None: - self.psf.save(True) - self.psf.close() - self.psf = None - - # Close energy dispersion - if self.edisp != None: - self.edisp.save(True) - self.edisp.close() + + # If IRF file exists then close it now. All file pointers will + # be set to None + if self.irf != None: + self.irf.save(True) + self.irf.close() + self.irf = None + self.ea = None + self.psf = None self.edisp = None + self.bgd = None + + # ... otherwise we have split files, so we have to close them + # all + else: + + # Close effective area file + if self.ea != None: + self.ea.save(True) + self.ea.close() + self.ea = None + + # Close point spread function file + if self.psf != None: + self.psf.save(True) + self.psf.close() + self.psf = None + + # Close energy dispersion file + if self.edisp != None: + self.edisp.save(True) + self.edisp.close() + self.edisp = None + + # Close background file + if self.bgd != None: + self.bgd.save(True) + self.bgd.close() + self.bgd = None # Return return def create_cif(self): """ - Create CIF extension in FITS file. + Create Calibration Index File (CIF) extension in CIF FITS file. + + Parameters: + None + Keywords: + None """ # Create binary table table = GFitsBinTable() @@ -248,6 +353,11 @@ class caldb(): def add_cif_info(self): """ Add information to CIF extension. + + Parameters: + None + Keywords: + None """ # Append 3 rows to CIF extension self.hdu_cif.append_rows(3) @@ -310,20 +420,35 @@ class caldb(): self.hdu_cif["CAL_CBD"][row, n] = "NONE" else: self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] + + # Add background information + row = self.hdu_cif.nrows()-1 + self.hdu_cif["CAL_DIR"][row] = self.bgd_dir + self.hdu_cif["CAL_FILE"][row] = self.bgd_file + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.bgd_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] # Return return - def open_ea(self, filename): + def open_ea(self): """ Open effective area extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.ea = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_ea = self.ea.table("EFFECTIVE AREA") except: @@ -346,16 +471,18 @@ class caldb(): # Return return - def open_psf(self, filename): + def open_psf(self): """ Open point spread function extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.psf = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") except: @@ -378,16 +505,18 @@ class caldb(): # Return return - def open_edisp(self, filename): + def open_edisp(self): """ Open energy dispersion extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.edisp = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") except: @@ -410,9 +539,50 @@ class caldb(): # Return return + def open_bgd(self): + """ + Open background extension. If no extension was found + then create one. We do not yet set any data as we don't + know the dimensions of the table yet. + + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. + try: + self.hdu_bgd = self.bgd.table("BACKGROUND") + except: + # Create binary table + table = GFitsBinTable() + + # Set extension name + table.extname("BACKGROUND") + + # Set OGIP keywords + self.set_ogip_keywords(table, self.bgd_doc, \ + ["RESPONSE", self.bgd_name]) + + # Append table + self.bgd.append(table) + + # Get extension + self.hdu_bgd = self.bgd.table("BACKGROUND") + + # Return + return + def set_ogip_keywords(self, hdu, hdudoc, hduclas): """ Set standard OGIP keywords for extension. + + Parameters: + hdu - Header Data Unit. + hdudoc - Documentation reference string + hduclas - Array of HDUCLAS fields + Keywords: + None """ # Set keywords hdu.card("ORIGIN", "IRAP", "Name of organization making this file") @@ -433,12 +603,17 @@ class caldb(): def make_2D(self, array, hdu, name, unit, scale=1.0): """ Make 2D matrix as function of energy and offset angle from a - ROOT 2D histogram. + ROOT 2D histogram. If the HDU has already the energy and + offset angle columns, this method will simply add another data + column. + If name==None, the method will not append any data column. Parameters: - array - ROOT 2D histgram. + array - ROOT 2D histogram. hdu - FITS HDU. - Keyword: + name - Data column name. + unit - Data unit. + Keywords: scale - Scaling factor for histogram values. """ # Extract energy and offset angle vectors @@ -485,7 +660,7 @@ class caldb(): hdu["THETA_HI"][0,ioff] = o_hi # # "NAME" - if not hdu.hascolumn(name): + if name != None and not hdu.hascolumn(name): hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset)) hdu[name].unit(unit) hdu[name].dim([neng, noffset]) @@ -515,44 +690,217 @@ class caldb(): Parameters: filename - ROOT 2D performance filename. """ - # Open performance file + # Open ROOT performance file file = TFile(filename) - # Allocate ROOT 2D array - array = TH2F() + # Create effective area + self.root2ea(file) + + # Create point spread function + self.root2psf(file) + + # Create energy dispersion + self.root2edisp(file) + + # Create background + self.root2bgd(file) + + # Return + return - # Make 2D effective 80% area + def root2ea(self, file): + """ + Translate ROOT to CALDB effective area extension. The following ROOT + histograms are used: + + EffectiveAreaEtrue_offaxis -> EFFAREA + EffectiveArea_offaxis -> EFFAREA_RECO + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if effective area HDU has been opened if self.hdu_ea != None: - file.GetObject("EffectiveArea80_offaxis", array) - bounds = self.make_2D(array, self.hdu_ea, "EFFAREA", "m2", scale=1.0/0.8) + + # Allocate ROOT 2D arrays + etrue = TH2F() + ereco = TH2F() + file.GetObject("EffectiveAreaEtrue_offaxis", etrue) + file.GetObject("EffectiveArea_offaxis", ereco) + #print etrue.GetXaxis().GetBinLowEdge(1), etrue.GetXaxis().GetBinUpEdge(etrue.GetXaxis().GetNbins()) + #print ereco.GetXaxis().GetBinLowEdge(1), ereco.GetXaxis().GetBinUpEdge(ereco.GetXaxis().GetNbins()) + + # Rebin etrue histogram + etrue.RebinX(10) + neng = etrue.GetXaxis().GetNbins() + noffset = etrue.GetYaxis().GetNbins() + for ioff in range(noffset): + for ieng in range(neng): + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0 + etrue.SetBinContent(ieng+1,ioff+1,value) + + # Write boundaries (use Ereco boundaries) + bounds = self.make_2D(ereco, self.hdu_ea, None, "m2") for b in bounds: self.ea_bounds.append(b) self.set_cif_keywords(self.hdu_ea, self.ea_name, \ self.ea_bounds, self.ea_desc) - # Make 2D 68% point spread function + # EFFAREA + self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2") + + # EFFAREA_RECO + self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2") + + # Return + return + + def root2psf(self, file): + """ + Translate ROOT to CALDB point spread function extension. The following + ROOT histograms are used: + + 1/(2*pi*SIGMA_1) -> SCALE + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8) + 0.0 -> AMPL_2 + 0.0 -> SIGMA_2 + 0.0 -> AMPL_3 + 0.0 -> SIGMA_3 + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if point spread function HDU has been opened if self.hdu_psf != None: - file.GetObject("AngRes_offaxis", array) - bounds = self.make_2D(array, self.hdu_psf, "R68", "deg") + + # Allocate ROOT 2D array + r68 = TH2F() + + # Read 68% containment histogram + file.GetObject("AngRes_offaxis", r68) + neng = r68.GetXaxis().GetNbins() + noffset = r68.GetYaxis().GetNbins() + + # Compute scale histogram + scale = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1) + if integral > 0.0: + value = 1.0 / integral + else: + value = 0.0 + scale.SetBinContent(ieng+1,ioff+1,value) + + # Set zero histogram + zero = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + zero.SetBinContent(ieng+1,ioff+1,0.0) + + # Set boundaries + bounds = self.make_2D(r68, self.hdu_psf, None, "deg") for b in bounds: self.psf_bounds.append(b) self.set_cif_keywords(self.hdu_psf, self.psf_name, \ self.psf_bounds, self.psf_desc) - # Make 2D 80% point spread function - if self.hdu_psf != None: - file.GetObject("AngRes80_offaxis", array) - self.make_2D(array, self.hdu_psf, "R80", "deg") + # SCALE + self.make_2D(scale, self.hdu_psf, "SCALE", "") + + # SIGMA_1 + self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg") - # Make 2D energy resolution + # AMPL_2 + self.make_2D(zero, self.hdu_psf, "AMPL_2", "") + + # SIGMA_2 + self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg") + + # AMPL_3 + self.make_2D(zero, self.hdu_psf, "AMPL_3", "") + + # SIGMA_3 + self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg") + + # Return + return + + def root2edisp(self, file): + """ + Translate ROOT to CALDB energy dispersion extension. The following ROOT + histograms are used: + + ERes_offaxis -> ERES + Ebias_offaxis -> EBIAS + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if energy dispersion HDU has been opened if self.hdu_edisp != None: - file.GetObject("ERes", array) - bounds = self.make_2D(array, self.hdu_edisp, "ERES68", "") + + # Allocate ROOT 2D array + eres = TH2F() + ebias = TH2F() + file.GetObject("ERes_offaxis", eres) + file.GetObject("Ebias_offaxis", ebias) + + # Set boundaries + bounds = self.make_2D(eres, self.hdu_edisp, None, "deg") for b in bounds: self.edisp_bounds.append(b) self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ self.edisp_bounds, self.edisp_desc) - + + # ERES + self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue") + + # EBIAS + self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue") + + # Return + return + + def root2bgd(self, file): + """ + Translate ROOT to CALDB background extension. The following ROOT + histograms are used: + + BGRatePerSqDeg_offaxis -> BGD + BGRatePerSqDeg_offaxis -> BGD_RECO + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if background HDU has been opened + if self.hdu_bgd != None: + + # Allocate ROOT 2D array + array = TH2F() + file.GetObject("BGRatePerSqDeg_offaxis", array) + + # Set boundaries + bounds = self.make_2D(array, self.hdu_bgd, None, "deg") + for b in bounds: + self.bgd_bounds.append(b) + self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \ + self.bgd_bounds, self.bgd_desc) + + # BGD + self.make_2D(array, self.hdu_bgd, "BGD", "") + + # BGD_RECO + self.make_2D(array, self.hdu_bgd, "BGD_RECO", "") + # Return return @@ -589,22 +937,41 @@ class caldb(): #==========================# if __name__ == '__main__': """ - This script translates a ROOT 2D performance file in a CALDB compliant - CTA response. + This script translates a CTA ROOT 2D performance files in a CALDB compliant + CTA response. The script will create the response for 3 observations to + illustrate the run-wise file structure. + + The CTA ROOT 2D performance files contain the following histograms: + - DiffSens_offaxis (2D) - Differential sensitivity + - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?) + - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment + - AngRes_offaxis (2D) - Angular resolution 68% containment + - AngRes80_offaxis (2D) - Angular resolution 80% containment + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree + - BGRate_offaxis (2D) - Background rate + - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins) + - MigMatrix_offaxis (3D) + - EestOverEtrue_offaxis (3D) + - ERes_offaxis (2D) - Energy Resolution + - Ebias_offaxis (2D) - Energy bias """ # Set ROOT filename - #path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxis" - #irf = "SubarrayE_offaxis.root" - path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_Nov2011" - irf = "SubarrayE_IFAE_50hours_20111121_offaxis.root" + path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_May2012" + irf = "SubarrayE_IFAE_50hours_20120510_offaxis.root" filename = path+"/"+irf - # Allocate caldb - db = caldb("E") + # Set observation identifiers + obsids = ["000001", "000002", "000003"] + + # Generate CALDB entries, one for each observation (or run) + for obsid in obsids: + + # Allocate caldb + db = caldb("E", obsid) - # Open calibration - db.open("test") + # Open calibration + db.open("test") - # Translate ROOT to CALDB information - db.root2caldb(filename) + # Translate ROOT to CALDB information + db.root2caldb(filename) \ No newline at end of file diff --git a/scripts/obsutils.py b/scripts/obsutils.py index ba1a37d..02c75ae 100755 --- a/scripts/obsutils.py +++ b/scripts/obsutils.py @@ -183,7 +183,7 @@ def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ event = cast_GCTAEventAtom(event) # Determine sky pixel - skydir = event.dir().skydir() + skydir = event.dir().dir() pixel = map.dir2pix(skydir) # Set pixel @@ -237,7 +237,7 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ # Get sky direction skydir = map.pix2dir(pixel) - instdir.skydir(skydir) + instdir.dir(skydir) # Create event atom for map pixel atom = GCTAEventAtom() @@ -295,7 +295,7 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ # Set ROI roi = GCTARoi() instdir = GCTAInstDir() - instdir.skydir(pntdir) + instdir.dir(pntdir) roi.centre(instdir) roi.radius(rad) diff --git a/src/Makefile.am b/src/Makefile.am index f7f1cd5..8a05e37 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -5,10 +5,11 @@ MACOSX_DEPLOYMENT_TARGET=@CONFIGURE_MACOSX_DEPLOYMENT_TARGET@ @EXPORT_MACOSX_DEPLOYMENT_TARGET@export MACOSX_DEPLOYMENT_TARGET # Set sub directories -SUBDIRS = ctselect \ - ctbin \ +SUBDIRS = ctbin \ ctlike \ + ctmodel \ ctobssim \ + ctselect \ ctskymap # Scripts to be distributed @@ -22,9 +23,10 @@ scripts = ctools-setup \ lib_LTLIBRARIES = libctools.la libctools_la_SOURCES = libctools_la_LIBADD = ctbin/libctbin.la \ - ctselect/libctselect.la \ ctlike/libctlike.la \ + ctmodel/libctmodel.la \ ctobssim/libctobssim.la \ + ctselect/libctselect.la \ ctskymap/libctskymap.la # Distribute scripts diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 15bdecb..cb431d1 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -534,7 +534,7 @@ void ctbin::bin_events(GCTAObservation* obs) // Determine sky pixel GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); - GSkyDir dir = inst->skydir(); + GSkyDir dir = inst->dir(); GSkyPixel pixel = map.dir2xy(dir); // Skip if pixel is out of range diff --git a/src/ctmodel/Makefile.am b/src/ctmodel/Makefile.am new file mode 100644 index 0000000..a968db2 --- /dev/null +++ b/src/ctmodel/Makefile.am @@ -0,0 +1,21 @@ +## Process this file with automake to produce Makefile.in + +# Program name +bin_PROGRAMS = ctmodel + +# Sources +sources = ctmodel.cpp + +# Libtool library +noinst_LTLIBRARIES = libctmodel.la +libctmodel_la_SOURCES = $(sources) + +# Program definition +ctmodel_SOURCES = main.cpp +ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctmodel_CPPFLAGS = @CPPFLAGS@ +ctmodel_LDADD = libctmodel.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctmodel.par diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp new file mode 100644 index 0000000..7fca3d6 --- /dev/null +++ b/src/ctmodel/ctmodel.cpp @@ -0,0 +1,873 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by 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 * + * 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 ctmodel.cpp + * @brief CTA counts model tool implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctmodel.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_SETUP_OBS "ctmodel::setup_obs()" +#define G_MODEL_MAP "ctmodel::model_map(GCTAObservation*)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * This method creates an instance of the class by copying an existing + * observations container. + ***************************************************************************/ +ctmodel::ctmodel(GObservations obs) : GApplication(CTMODEL_NAME, CTMODEL_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // 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. + ***************************************************************************/ +ctmodel::ctmodel(int argc, char *argv[]) : + GApplication(CTMODEL_NAME, CTMODEL_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctmodel::ctmodel(const ctmodel& app) : GApplication(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctmodel::~ctmodel(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app ctmodel application. + * @return Returns ctmodel application. + ***************************************************************************/ +ctmodel& ctmodel::operator= (const ctmodel& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->GApplication::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 instance + ***************************************************************************/ +void ctmodel::clear(void) +{ + // Free members + free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Execute application + * + * This is the main execution method of the ctmodel class. It is invoked + * when the executable is called from command line. The method generates + * the model maps and saves the results. + ***************************************************************************/ +void ctmodel::execute(void) +{ + // Signal that some parameters should be read ahead + m_read_ahead = true; + + // Create the model map(s) + run(); + + // Save the model map(s) into FITS file + save(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate the model map(s) + * + * This method reads the task parameters from the parfile, sets up the + * observation container, loops over all CTA observations in the container + * and generates a model map for each CTA observation. + ***************************************************************************/ +void ctmodel::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(); + + // Setup observation container + setup_obs(); + + // 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; + if (m_obs.size() > 1) { + log.header1("Generate model maps"); + } + else { + log.header1("Generate model map"); + } + } + + // Initialise observation counter + int n_observations = 0; + + // Loop over all observations in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Initialise event input and output filenames + m_infiles.push_back(""); + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // Write header for observation + if (logTerse()) { + if (obs->name().length() > 1) { + log.header3("Observation "+obs->name()); + } + else { + log.header3("Observation"); + } + } + + // Increment number of observations + n_observations++; + + // Save event file name (for possible saving) + m_infiles[i] = obs->eventfile(); + + // Generate model map + model_map(obs, m_obs.models()); + + } // endif: CTA observation found + + } // endfor: looped over observations + + // If more than a single observation has been handled then make sure + // that an XML file will be used for storage + if (n_observations > 1) { + m_use_xml = true; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Observations after model map generation"); + } + else { + log.header1("Observation after model map generation"); + } + log << m_obs << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save model map(s) + * + * This method saves the model map(s) into FITS file(s). There are two + * modes, depending on the m_use_xml flag. + * + * If m_use_xml is true, all model map(s) will be saved into FITS files, + * where the output filenames are constructued from the input filenames by + * prepending the m_prefix string to name. Any path information will be + * stripped form the input name, hence event files will be written into the + * local working directory (unless some path information is present + * in the prefix). In addition, an XML file will be created that gathers + * the filename information for the model map(s). If an XML file was present + * on input, all metadata information will be copied from this input file. + * + * If m_use_xml is false, the model map will be saved into a FITS file. + ***************************************************************************/ +void ctmodel::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Save observations"); + } + else { + log.header1("Save observation"); + } + } + + // Case A: Save model map(s) and XML metadata information + if (m_use_xml) { + save_xml(); + } + + // Case B: Save model map as FITS file + else { + save_fits(); + } + + // 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 ctmodel::get_parameters(void) +{ + // If we do not have any observations in the container then get all + // information to build (at least) one + if (m_obs.size() == 0) { + + // Read general parameters. Only read the output parameters if + // required, and also only read the model XML filename if no + // model is yet in the container. + m_infile = (*this)["infile"].filename(); + if (m_read_ahead) { + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + } + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + if (m_obs.models().size() == 0) { + m_srcmdl = (*this)["srcmdl"].filename(); + } + + // If there is no input filename then read all parameters that + // are required to build a model map from scratch + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) { + m_ra = (*this)["ra"].real(); + m_dec = (*this)["dec"].real(); + m_deadc = (*this)["deadc"].real(); + m_tmin = (*this)["tmin"].real(); + m_tmax = (*this)["tmax"].real(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_enumbins = (*this)["enumbins"].integer(); + m_proj = (*this)["proj"].string(); + m_coordsys = (*this)["coordsys"].string(); + m_xref = (*this)["xref"].real(); + m_yref = (*this)["yref"].real(); + m_binsz = (*this)["binsz"].real(); + m_nxpix = (*this)["nxpix"].integer(); + m_nypix = (*this)["nypix"].integer(); + } + } + + // ... otherwise, read only the parameters that are required for the + // model generation. + else { + if (m_read_ahead) { + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + } + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + if (m_obs.models().size() == 0) { + m_srcmdl = (*this)["srcmdl"].filename(); + } + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Setup observation container + * + * @exception GException::no_cube + * No event cube found in CTA observation. + * + * This method sets up the observation container for processing. There are + * two cases: + * + * If there are no observations in the actual observation container, the + * method will check in "infile" parameter. If this parameter is "NONE" or + * empty, the task parameters will be used to construct a model map. + * Otherwise, the method first tries to interpret the "infile" parameter as + * a counts map, and attemps loading of the file in an event cube. If this + * fails, the method tries to interpret the "infile" parameter as an + * observation definition XML file. If this also fails, an exception will + * be thrown. + * + * If observations exist already in the observation container, the method + * will simply keep them. + * + * Test if all CTA observations contain counts maps. + * + * Finally, if no models exist so far in the observation container, the + * models will be loaded from the model XML file. + ***************************************************************************/ +void ctmodel::setup_obs(void) +{ + // If there are no observations in the container then try to build some + if (m_obs.size() == 0) { + + // If no input filename has been specified, then create a model map + // from the task parameters + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) { + + // Set pointing direction + GCTAPointing pnt; + GSkyDir skydir; + skydir.radec_deg(m_ra, m_dec); + pnt.dir(skydir); + + // Setup energy range covered by model + GEnergy emin; + GEnergy emax; + GEbounds ebds; + emin.TeV(m_emin); + emax.TeV(m_emax); + ebds.setlog(emin, emax, m_enumbins); + + // Setup time interval covered by model + GGti gti; + GTime tmin; + GTime tmax; + tmin.met(m_tmin); + tmax.met(m_tmax); + gti.append(tmin, tmax); + + // Setup skymap + GSkymap map = GSkymap(m_proj, m_coordsys, + m_xref, m_yref, m_binsz, m_binsz, + m_nxpix, m_nypix, m_enumbins); + + // Create model cube from sky map + GCTAEventCube cube(map, ebds, gti); + + // Allocate CTA observation + GCTAObservation obs; + + // Set CTA observation attributes + obs.pointing(pnt); + obs.ontime(gti.ontime()); + obs.livetime(gti.ontime()*m_deadc); + obs.deadc(m_deadc); + + // Set event cube in observation + obs.events(&cube); + + // Append CTA observation to container + m_obs.append(obs); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } // endif: created model map from task parameters + + // ... otherwise try to load information from the file + else { + + // First try to open the file as a counts map + try { + + // Allocate CTA observation + GCTAObservation obs; + + // Load counts map in CTA observation + obs.load_binned(m_infile); + + // Append CTA observation to container + m_obs.append(obs); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise try to open as XML file + catch (GException::fits_open_error &e) { + + // Load observations from XML file. This will throw + // an exception if it fails. + m_obs.load(m_infile); + + // Signal that XML file should be used for storage + m_use_xml = true; + + } + + } // endelse: loaded information from input file + + } // endif: there was no observation in the container + + // If there are no models associated with the observations then + // load now the model definition from the XML file + if (m_obs.models().size() == 0) { + m_obs.models(GModels(m_srcmdl)); + } + + // Check if all CTA observations contain an event cube and setup response + // for all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Is this observation a CTA observation? + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Yes ... + if (obs != NULL) { + + // Throw an exception if this observation does not contain + // an event cube + if (dynamic_cast(obs->events()) == NULL) { + throw GException::no_cube(G_SETUP_OBS); + } + + // Set response + obs->response(m_irf, m_caldb); + + } // endif: observation was a CTA observation + } // endfor: looped over all observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate model map + * + * @param[in] obs CTA observation pointer. + * @param[in] models Model container. + * + * @exception GException::no_cube + * No event cube found in CTA observation. + ***************************************************************************/ +void ctmodel::model_map(GCTAObservation* obs, const GModels& models) +{ + // Continue only if observation pointer is valid + if (obs != NULL) { + + // Get event cube pointer + GCTAEventCube* cube = + const_cast(dynamic_cast(obs->events())); + + // Throw an exception if the observation does not hold and event + // cube + if (cube == NULL) { + throw GException::no_cube(G_MODEL_MAP); + } + + // Initialise statistics + double sum = 0.0; + + // Loop over all events in counts map + for (int i = 0; i < cube->size(); ++i) { + + // Get event bin + GCTAEventBin* bin = (*cube)[i]; + + // Compute model value for event bin + double model = + models.eval(*(const_cast(bin)), *obs) * + bin->size(); + + // Store value + bin->counts(model); + + // Sum all events + sum += model; + } + + // Log results + if (logTerse()) { + log << parformat("Model events in cube"); + log << sum << std::endl; + } + + } // endif: observation pointer was not valid + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctmodel::init_members(void) +{ + // Initialise members + m_infile.clear(); + m_outfile.clear(); + m_prefix.clear(); + m_caldb.clear(); + m_irf.clear(); + m_srcmdl.clear(); + m_proj.clear(); + m_coordsys.clear(); + m_ra = 0.0; + m_dec = 0.0; + m_deadc = 1.0; + m_tmin = 0.0; + m_tmax = 0.0; + m_emin = 0.0; + m_emax = 0.0; + m_enumbins = 0; + m_xref = 0.0; + m_yref = 0.0; + m_binsz = 0.0; + m_nxpix = 0; + m_nypix = 0; + + // Initialise protected members + m_obs.clear(); + m_infiles.clear(); + m_use_xml = false; + m_read_ahead = false; + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctmodel::copy_members(const ctmodel& app) +{ + // Copy attributes + m_infile = app.m_infile; + m_outfile = app.m_outfile; + m_prefix = app.m_prefix; + m_caldb = app.m_caldb; + m_irf = app.m_irf; + m_srcmdl = app.m_srcmdl; + m_proj = app.m_proj; + m_coordsys = app.m_coordsys; + m_ra = app.m_ra; + m_dec = app.m_dec; + m_deadc = app.m_deadc; + m_tmin = app.m_tmin; + m_tmax = app.m_tmax; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_enumbins = app.m_enumbins; + m_xref = app.m_xref; + m_yref = app.m_yref; + m_binsz = app.m_binsz; + m_nxpix = app.m_nxpix; + m_nypix = app.m_nypix; + + // Copy protected members + m_obs = app.m_obs; + m_infiles = app.m_infiles; + m_use_xml = app.m_use_xml; + m_read_ahead = app.m_read_ahead; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctmodel::free_members(void) +{ + // Write separator into logger + if (logTerse()) { + log << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set output file name. + * + * @param[in] filename Input file name. + * + * This method converts an input filename into an output filename by + * prepending the prefix specified by m_prefix to the input filename. Any + * path will be stripped from the input filename. + ***************************************************************************/ +std::string ctmodel::set_outfile_name(const std::string& filename) const +{ + // Split input filename into path elements + std::vector elements = split(filename, "/"); + + // The last path element is the filename + std::string outname = m_prefix + elements[elements.size()-1]; + + // Return output filename + return outname; +} + + +/***********************************************************************//** + * @brief Save model map in FITS format. + * + * Save the model map as a FITS file. The filename of the FITS file is + * specified by the m_outfile member. + ***************************************************************************/ +void ctmodel::save_fits(void) +{ + // Get output filename + m_outfile = (*this)["outfile"].filename(); + + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(&m_obs[0]); + + // Save model map + save_model_map(obs, m_outfile); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save model map(s) in XML format. + * + * Save the model map(s) into FITS files and write the file path information + * into a XML file. The filename of the XML file is specified by the + * m_outfile member, the filename(s) of the model map(s) are built by + * prepending the prefix given by the m_prefix member to the input model + * map(s) filenames. Any path present in the input filename will be stripped, + * i.e. the model map(s) will be written in the local working directory + * (unless a path is specified in the m_prefix member). + ***************************************************************************/ +void ctmodel::save_xml(void) +{ + // Get output filename and prefix + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Handle only CTA observations + if (obs != NULL) { + + // Set event output file name + std::string outfile = set_outfile_name(m_infiles[i]); + + // Store output file name in observation + obs->eventfile(outfile); + + // Save event list + save_model_map(obs, outfile); + + } // endif: observation was a CTA observations + + } // endfor: looped over observations + + // Save observations in XML file + m_obs.save(m_outfile); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save a single model map into a FITS file + * + * @param[in] obs Pointer to CTA observation. + * @param[in] outfile Output file name. + * + * This method saves a single model map into a FITS file. The method does + * nothing if the observation pointer is not valid. + ***************************************************************************/ +void ctmodel::save_model_map(const GCTAObservation* obs, + const std::string& outfile) const +{ + // Save only if observation is valid + if (obs != NULL) { + + // Save observation into FITS file + obs->save(outfile, clobber()); + + } // endif: observation was valid + + // Return + return; +} diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp new file mode 100644 index 0000000..7afcd43 --- /dev/null +++ b/src/ctmodel/ctmodel.hpp @@ -0,0 +1,117 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by 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 * + * 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 ctmodel.hpp + * @brief CTA counts model tool definition + * @author Juergen Knoedlseder + */ + +#ifndef CTMODEL_HPP +#define CTMODEL_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTMODEL_NAME "ctmodel" +#define CTMODEL_VERSION "00-01-00" + + +/***********************************************************************//** + * @class ctmodel + * + * @brief CTA counts model tool interface defintion + * + * This class creates counts model map(s). The definition of the counts model + * can be taken from a predefined observation container, from the counts maps + * found in an observation definition XML file, from an individual counts + * map, or from the parameter definition. + * + * Results are stored in an observation container that can be written to disk + * in form of FITS files (model maps) and an updated observation definition + * XML file. + ***************************************************************************/ +class ctmodel : public GApplication { + +public: + // Constructors and destructors + ctmodel(void); + explicit ctmodel(GObservations obs); + ctmodel(int argc, char *argv[]); + ctmodel(const ctmodel& app); + virtual ~ctmodel(void); + + // Operators + ctmodel& operator= (const ctmodel& app); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + GObservations& obs(void) { return m_obs; } + void get_parameters(void); + void setup_obs(void); + void model_map(GCTAObservation* obs, const GModels& models); + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctmodel& app); + void free_members(void); + std::string set_outfile_name(const std::string& filename) const; + void save_fits(void); + void save_xml(void); + void save_model_map(const GCTAObservation* obs, + const std::string& outfile) const; + + // User parameters + std::string m_infile; //!< Input counts map or XML file + std::string m_outfile; //!< Output model map or XML file + std::string m_prefix; //!< Prefix for multiple model maps + std::string m_caldb; //!< Calibration database + std::string m_irf; //!< Instrument response functions + std::string m_srcmdl; //!< Source model + double m_ra; //!< RA of pointing direction + double m_dec; //!< DEC of pointing direction + double m_deadc; //!< Deadtime correction factor + double m_tmin; //!< Start time + double m_tmax; //!< Stop time + double m_emin; //!< Lower energy + double m_emax; //!< Upper energy + int m_enumbins; //!< Number of energy bins + std::string m_proj; //!< WCS projection + std::string m_coordsys; //!< Coordinate system + double m_xref; //!< Longitude reference coordinate + double m_yref; //!< Latitude reference coordinate + double m_binsz; //!< Pixel size + int m_nxpix; //!< Number of pixels in longitude + int m_nypix; //!< Number of pixels in latitude + + // Protected members + GObservations m_obs; //!< Observation container + std::vector m_infiles; //!< Input event filenames + bool m_use_xml; //!< Use XML file instead of FITS file + bool m_read_ahead; //!< Read ahead parameters +}; + +#endif /* CTMODEL_HPP */ diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par new file mode 100644 index 0000000..3388172 --- /dev/null +++ b/src/ctmodel/ctmodel.par @@ -0,0 +1,62 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctmodel.par # +# Version: 00-01-00 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the ctmodel application # +# # +############################################################################# +# +# General parameters +#=================== +infile, f, a, "NONE",,, "Input counts map or observation definition file" +outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file" +prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file" +caldb, s, a, "test/irf",,, "Calibration database" +irf, s, a, "kb_E_50h_v3",,, "Instrument response function" +srcmdl, f, a, "test/data/crab.xml",,, "Source model" +# +# Instrument parameters +#====================== +ra, r, a, 83.63, 0,360, "RA of pointing (degrees)" +dec, r, a, 22.01, -90, 90, "Dec of pointing (degrees)" +deadc, r, h, 0.95, 0, 1, "Average deadtime correction factor" +# +# Time range parameters +#====================== +tmin, r, a, 0.0, 0, , "Start time (MET in s)" +tmax, r, a, 1800.0, 0, , "End time (MET in s)" +# +# Energy binning parameters +#========================== +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" +emin, r, a, 0.1,,,"Start value for first energy bin in TeV" +emax, r, a, 100.0,,,"Stop value for last energy bin in 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 +#=========================== +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,,, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)" +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees" +proj, s, a, "CAR",,, "Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN" +# +# Standard parameters +#==================== +chatter, i, h, 1,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/ctmodel/main.cpp b/src/ctmodel/main.cpp new file mode 100644 index 0000000..7240904 --- /dev/null +++ b/src/ctmodel/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by 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 * + * 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 main.cpp + * @brief CTA counts model tool main code + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctmodel.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 gtbin application. It allocates a + * ctmodel object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctmodel 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" + " ctmodel. 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/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index e493742..d77379f 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -621,7 +621,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, GCTAEventList* events = static_cast(const_cast(obs->events())); // Extract ROI - GSkyDir dir = events->roi().centre().skydir(); + GSkyDir dir = events->roi().centre().dir(); double rad = events->roi().radius(); // Dump simulation cone information diff --git a/src/ctools-setup.in b/src/ctools-setup.in index e99df45..99a22aa 100644 --- a/src/ctools-setup.in +++ b/src/ctools-setup.in @@ -29,7 +29,7 @@ fi prefix="@prefix@" exec_prefix="@exec_prefix@" ctools_bin="@bindir@" -ctools_lib="@libdir@@TEST_ENV_DIR@" +ctools_lib="@libdir@" uname_system=`uname` python_dir="@pythondir@" python_exec_dir="@pyexecdir@" diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index b30f6a4..51bf608 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -492,7 +492,7 @@ void ctskymap::map_events(GCTAObservation* obs) // Determine sky pixel GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); - GSkyDir dir = inst->skydir(); + GSkyDir dir = inst->dir(); GSkyPixel pixel = m_skymap.dir2xy(dir); // Skip if pixel is out of range diff --git a/test/test_ctools.sh b/test/test_ctools.sh index c35058a..cdee21f 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -34,11 +34,12 @@ $ECHO "***************" # # Define executables # ================== -ctobssim=../src/ctobssim/ctobssim -ctskymap=../src/ctskymap/ctskymap ctbin=../src/ctbin/ctbin -ctselect=../src/ctselect/ctselect ctlike=../src/ctlike/ctlike +ctmodel=../src/ctmodel/ctmodel +ctobssim=../src/ctobssim/ctobssim +ctselect=../src/ctselect/ctselect +ctskymap=../src/ctskymap/ctskymap # @@ -136,6 +137,57 @@ $ECHO " ok" # +# Test ctmodel +# ============ +$ECHO -n "Test ctmodel: " +# +# Run 1 +$ctmodel infile="cntmap.fits" \ + outfile="modmap1.fits" \ + caldb="irf" \ + irf="kb_E_50h_v3" \ + srcmdl="data/crab.xml" +$ECHO -n "." +if [ -s "modmap1.fits" ] +then + $ECHO -n "." +else + $ECHO " modmap1.fits file is not valid" + exit 1 +fi +# +# Run 2 +$ctmodel infile="NONE" \ + outfile="modmap2.fits" \ + caldb="irf" \ + irf="kb_E_50h_v3" \ + srcmdl="data/crab.xml" \ + ra=83.63 \ + dec=22.01 \ + tmin=0.0 \ + tmax=1800.0 \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=200 \ + nypix=200 \ + binsz=0.02 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" +$ECHO -n "." +if [ -s "modmap2.fits" ] +then + $ECHO -n "." +else + $ECHO " modmap2.fits file is not valid" + exit 1 +fi +$ECHO " ok" + + +# # Test ctselect # ============= $ECHO -n "Test ctselect: "