Alexander Ziegler / ctools

Commit b930a2584089305409164fa136c1c4f0dc1c1ca6

Authored by Jürgen Knödlseder
2 parents 91a3d3f5 d5049a1a

Merge branch 'hess' into devel 

Conflicts:
	.gitignore
	ChangeLog
Inline Side-by-side
Showing 19 changed files with 1776 additions and 109 deletions
.gitignore
  View file @b930a25
... ... @@ -22,9 +22,13 @@ ltmain.sh
22 22 m4
23 23 missing
24 24 .libs
  25 +.deps
25 26 stamp-h1
  27 +depcomp
26 28  
  29 +src/ctools-setup
27 30 pyext/build/
  31 +pyext/setup.py
28 32 pyext/ctools.py
29 33 pyext/ctools_wrap.cpp
30 34 src/ctbin/ctbin
... ... @@ -34,3 +38,7 @@ src/ctobssim/ctobssim
34 38 src/ctselect/ctselect
35 39 src/ctskymap/ctskymap
36 40 test/pfiles/
  41 +test/*.fits
  42 +test/*.dat
  43 +test/*.xml
  44 +test/test_cscripts.sh
... ...
1 1 Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) - overall design
2 2 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) - openmp support
  3 +Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) - HESS support
  4 +Mathieu De Naurois (denauroi@in2p3.fr) - HESS support
3 5  
... ...
ChangeLog
  View file @b930a25
  1 +2012-12-05 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  2 +
  3 + * Version 0.6.0 released
  4 + ========================
  5 +
  6 + Add ctmodel tool
  7 + Adapt to GammaLib-00-07-00 interface
  8 +
1 9 2012-10-11 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 10  
3 11 * Bug fix version 0.5.2 released
... ...
configure.ac
  View file @b930a25
... ... @@ -32,7 +32,8 @@
32 32 AC_INIT([ctools], [0.6.0], [jurgen.knodlseder@irap.omp.eu], [ctools])
33 33 AC_CONFIG_SRCDIR([src/ctselect/ctselect.cpp])
34 34 AC_CONFIG_MACRO_DIR([m4])
35   -AM_INIT_AUTOMAKE([gnu no-dependencies])
  35 +#AM_INIT_AUTOMAKE([gnu no-dependencies])
  36 +AM_INIT_AUTOMAKE([gnu])
36 37 AM_CONFIG_HEADER([config.h])
37 38  
38 39  
... ... @@ -870,11 +871,12 @@ AC_SUBST(LIB_LINK_OPTION)
870 871 #############################################################################
871 872 AC_CONFIG_FILES([Makefile
872 873 src/Makefile
873   - src/ctselect/Makefile
874 874 src/ctbin/Makefile
875 875 src/ctlike/Makefile
  876 + src/ctmodel/Makefile
876 877 src/ctobssim/Makefile
877 878 src/ctskymap/Makefile
  879 + src/ctselect/Makefile
878 880 src/ctools-setup
879 881 test/Makefile
880 882 test/test_cscripts.sh
... ...
pyext/ctmodel.i 0 → 100644
  View file @b930a25
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctmodel.i
  23 + * @brief CTA counts model tool definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "ctmodel.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class ctmodel
  35 + *
  36 + * @brief CTA counts model tool interface defintion
  37 + ***************************************************************************/
  38 +class ctmodel : public GApplication {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + ctmodel(void);
  43 + explicit ctmodel(GObservations obs);
  44 + ctmodel(int argc, char *argv[]);
  45 + ctmodel(const ctmodel& app);
  46 + virtual ~ctmodel(void);
  47 +
  48 + // Methods
  49 + void clear(void);
  50 + void execute(void);
  51 + void run(void);
  52 + void save(void);
  53 + GObservations& obs(void) { return m_obs; }
  54 + void get_parameters(void);
  55 + void setup_obs(void);
  56 + void model_map(GCTAObservation* obs, const GModels& models);
  57 +};
  58 +
  59 +
  60 +/***********************************************************************//**
  61 + * @brief CTA counts model tool Python extension
  62 + ***************************************************************************/
  63 +%extend ctmodel {
  64 + ctmodel copy() {
  65 + return (*self);
  66 + }
  67 +}
... ...
pyext/ctools.i
  View file @b930a25
1 1 /***************************************************************************
2 2 * ctools - SWIG file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2011 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2012 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,6 +21,11 @@
21 21 * Usage: *
22 22 * swig -c++ -python -Wall ctools.i *
23 23 ***************************************************************************/
  24 +/**
  25 + * @file ctools.i
  26 + * @brief ctools SWIG file
  27 + * @author Juergen Knoedlseder
  28 + */
24 29 %module ctools
25 30 %feature("autodoc", "1");
26 31  
... ... @@ -31,10 +36,12 @@
31 36 %import(module="gammalib.support") "GException.i";
32 37  
33 38 /* __ Inform about base classes __________________________________________ */
34   -%import(module="gammalib.app") "GApplication.i";
  39 +%import(module="gammalib.base") "GBase.i";
  40 +%import(module="gammalib.app") "GApplication.i";
35 41  
36   -/* __ CTA tools __________________________________________________________ */
37   -%include "ctobssim.i"
38   -%include "ctselect.i"
  42 +/* __ ctools _____________________________________________________________ */
39 43 %include "ctbin.i"
  44 +%include "ctobssim.i"
40 45 %include "ctlike.i"
  46 +%include "ctmodel.i"
  47 +%include "ctselect.i"
... ...
scripts/cta_root2caldb.py
  View file @b930a25
... ... @@ -2,11 +2,9 @@
2 2 # ==========================================================================
3 3 # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D
4 4 # performance files.
5   -#
6   -# Requirements:
7   -# - pyroot
8   -#
9   -# Copyright (C) 2011 Jurgen Knodlseder
  5 +# -------------------------------------------------------------------------
  6 +# Copyright (C) 2011-2012 Juergen Knoedlseder
  7 +# -------------------------------------------------------------------------
10 8 #
11 9 # This program is free software: you can redistribute it and/or modify
12 10 # it under the terms of the GNU General Public License as published by
... ... @@ -21,6 +19,10 @@
21 19 # You should have received a copy of the GNU General Public License
22 20 # along with this program. If not, see <http://www.gnu.org/licenses/>.
23 21 #
  22 +# -------------------------------------------------------------------------
  23 +# Requirements:
  24 +# - gammalib
  25 +# - pyroot
24 26 # ==========================================================================
25 27 from ROOT import TFile, TH1F, TH2F
26 28 from gammalib import *
... ... @@ -37,18 +39,27 @@ import os
37 39 class caldb():
38 40 """
39 41 """
40   - def __init__(self, conf, path=""):
  42 + def __init__(self, inst, obsid, path=""):
41 43 """
42   - Constructor.
  44 + CALDB constructor. The constructor will create one CALDB entry for
  45 + an observation.
  46 +
  47 + Parameters:
  48 + inst - Instrument string.
  49 + obsid - Observation ID string.
  50 + Keywords:
  51 + path - CALDB Root path.
43 52 """
44 53 # Store root path
45 54 self.path = path
46 55  
47 56 # Initialise some members
48 57 self.cif = None
  58 + self.irf = None
49 59 self.ea = None
50 60 self.psf = None
51 61 self.edisp = None
  62 + self.bgd = None
52 63 self.hdu_cif = None
53 64  
54 65 # Store UTC date of job execution
... ... @@ -56,7 +67,8 @@ class caldb():
56 67  
57 68 # Set CALDB information
58 69 self.cal_tel = "CTA"
59   - self.cal_inst = conf.upper()
  70 + self.cal_inst = inst.upper()
  71 + self.cal_obsid = obsid
60 72 self.cal_det = "NONE"
61 73 self.cal_flt = "NONE"
62 74 self.cal_class = "BCF"
... ... @@ -84,6 +96,11 @@ class caldb():
84 96 self.edisp_doc = "???"
85 97 self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis]
86 98 self.edisp_desc = "CTA energy dispersion"
  99 + #
  100 + self.bgd_name = "BGD"
  101 + self.bgd_doc = "???"
  102 + self.bgd_bounds = [self.cal_version, self.cal_cut, self.cal_analysis]
  103 + self.bgd_desc = "CTA background"
87 104  
88 105 # Create directory structure
89 106 self.make_dirs()
... ... @@ -93,6 +110,7 @@ class caldb():
93 110  
94 111 def __del__(self):
95 112 """
  113 + Destructor. Close all FITS files.
96 114 """
97 115 # Close calibration
98 116 self.close()
... ... @@ -102,32 +120,50 @@ class caldb():
102 120  
103 121 def make_dirs(self):
104 122 """
  123 + Generate CALDB directory structure for one observation identifier.
  124 + The structure is given by
  125 +
  126 + data/<tel>/<inst>/bcf/<obsid>
  127 +
  128 + where <tel> is "cta" and <inst> is the instrument specified in the
  129 + CALDB constructor (the instrument may be used for different array
  130 + configurations).
  131 +
  132 + Parameters:
  133 + None
  134 + Keywords:
  135 + None
105 136 """
106 137 # Set base path
107 138 self.base_dir = "data"
108 139 self.base_dir += "/"+self.cal_tel.lower()
109 140 self.base_dir += "/"+self.cal_inst.lower()
110 141  
111   -
112 142 # Set directory names
113 143 self.bcf_dir = self.base_dir+"/bcf"
114   - self.ea_dir = self.bcf_dir+"/ea"
115   - self.psf_dir = self.bcf_dir+"/psf"
116   - self.edisp_dir = self.bcf_dir+"/edisp"
  144 + self.obs_dir = self.bcf_dir+"/"+self.cal_obsid
  145 + self.ea_dir = self.obs_dir
  146 + self.psf_dir = self.obs_dir
  147 + self.edisp_dir = self.obs_dir
  148 + self.bgd_dir = self.obs_dir
117 149  
118   - # Add path (if required)
  150 + # Optionally prefix path
119 151 if len(self.path) > 0:
120 152 self.base_path = self.path+"/"+self.base_dir
121 153 self.bcf_path = self.path+"/"+self.bcf_dir
  154 + self.obs_path = self.path+"/"+self.obs_dir
122 155 self.ea_path = self.path+"/"+self.ea_dir
123 156 self.psf_path = self.path+"/"+self.psf_dir
124 157 self.edisp_path = self.path+"/"+self.edisp_dir
  158 + self.bgd_path = self.path+"/"+self.bgd_dir
125 159 else:
126 160 self.base_path = self.base_dir
127 161 self.bcf_path = self.bcf_dir
  162 + self.obs_path = self.obs_dir
128 163 self.ea_path = self.ea_dir
129 164 self.psf_path = self.psf_dir
130 165 self.edisp_path = self.edisp_dir
  166 + self.bgd_path = self.bgd_dir
131 167  
132 168 # Create directory structure
133 169 if not os.path.isdir(self.ea_path):
... ... @@ -136,18 +172,36 @@ class caldb():
136 172 os.makedirs(self.psf_path)
137 173 if not os.path.isdir(self.edisp_path):
138 174 os.makedirs(self.edisp_path)
  175 + if not os.path.isdir(self.bgd_path):
  176 + os.makedirs(self.bgd_path)
139 177  
140 178 # Return
141 179 return
142 180  
143   - def open(self, name, clobber=True):
  181 + def open(self, version, split=False, clobber=True):
144 182 """
145   - Open existing or create new calibration.
  183 + Open existing or create new calibration. The actual version will put
  184 + all calibrations in the same file, although each part of the response
  185 + function will have its own logical name. We can thus easily modify
  186 + the script to put each calibration information in a separate file.
  187 +
  188 + Parameters:
  189 + version - Filename version
  190 + Keywords:
  191 + split - Split IRF over several files?
  192 + clobber - Overwrite existing files?
146 193 """
147 194 # Set calibrate file names
148   - self.ea_file = "ea_"+name+".fits"
149   - self.psf_file = "psf_"+name+".fits"
150   - self.edisp_file = "edisp_"+name+".fits"
  195 + if split:
  196 + self.ea_file = "ea_"+version+".fits"
  197 + self.psf_file = "psf_"+version+".fits"
  198 + self.edisp_file = "edisp_"+version+".fits"
  199 + self.bgd_file = "bgd_"+version+".fits"
  200 + else:
  201 + self.ea_file = "irf_"+version+".fits"
  202 + self.psf_file = "irf_"+version+".fits"
  203 + self.edisp_file = "irf_"+version+".fits"
  204 + self.bgd_file = "irf_"+version+".fits"
151 205  
152 206 # Open calibration database index
153 207 self.cif = GFits(self.base_path+"/caldb.indx", True)
... ... @@ -159,17 +213,42 @@ class caldb():
159 213 self.create_cif()
160 214 self.hdu_cif = self.cif.table("CIF")
161 215  
162   - # Open calibration files
163   - self.open_ea(self.ea_path+"/"+self.ea_file)
164   - self.open_psf(self.psf_path+"/"+self.psf_file)
165   - self.open_edisp(self.edisp_path+"/"+self.edisp_file)
  216 + # Set filenames
  217 + ea_filename = self.ea_path+"/"+self.ea_file
  218 + psf_filename = self.psf_path+"/"+self.psf_file
  219 + edisp_filename = self.edisp_path+"/"+self.edisp_file
  220 + bgd_filename = self.bgd_path+"/"+self.bgd_file
  221 +
  222 + # Open files
  223 + if split:
  224 + self.ea = GFits(ea_filename, True)
  225 + self.psf = GFits(psf_filename, True)
  226 + self.edisp = GFits(edisp_filename, True)
  227 + self.bgd = GFits(bgd_filename, True)
  228 + else:
  229 + self.irf = GFits(ea_filename, True)
  230 + self.ea = self.irf
  231 + self.psf = self.irf
  232 + self.edisp = self.irf
  233 + self.bgd = self.irf
  234 +
  235 + # Open HDUs
  236 + self.open_ea()
  237 + self.open_psf()
  238 + self.open_edisp()
  239 + self.open_bgd()
166 240  
167 241 # Return
168 242 return
169 243  
170 244 def close(self):
171 245 """
172   - Close calibration.
  246 + Close calibration FITS files.
  247 +
  248 + Parameters:
  249 + None
  250 + Keywords:
  251 + None
173 252 """
174 253 # Add information to CIF. We do this now as before we did not
175 254 # necessarily have all the information at hand (in particular
... ... @@ -181,31 +260,57 @@ class caldb():
181 260 self.cif.save(True)
182 261 self.cif.close()
183 262 self.cif = None
184   -
185   - # Close effective area
186   - if self.ea != None:
187   - self.ea.save(True)
188   - self.ea.close()
189   - self.ea = None
190   -
191   - # Close point spread function
192   - if self.psf != None:
193   - self.psf.save(True)
194   - self.psf.close()
195   - self.psf = None
196   -
197   - # Close energy dispersion
198   - if self.edisp != None:
199   - self.edisp.save(True)
200   - self.edisp.close()
  263 +
  264 + # If IRF file exists then close it now. All file pointers will
  265 + # be set to None
  266 + if self.irf != None:
  267 + self.irf.save(True)
  268 + self.irf.close()
  269 + self.irf = None
  270 + self.ea = None
  271 + self.psf = None
201 272 self.edisp = None
  273 + self.bgd = None
  274 +
  275 + # ... otherwise we have split files, so we have to close them
  276 + # all
  277 + else:
  278 +
  279 + # Close effective area file
  280 + if self.ea != None:
  281 + self.ea.save(True)
  282 + self.ea.close()
  283 + self.ea = None
  284 +
  285 + # Close point spread function file
  286 + if self.psf != None:
  287 + self.psf.save(True)
  288 + self.psf.close()
  289 + self.psf = None
  290 +
  291 + # Close energy dispersion file
  292 + if self.edisp != None:
  293 + self.edisp.save(True)
  294 + self.edisp.close()
  295 + self.edisp = None
  296 +
  297 + # Close background file
  298 + if self.bgd != None:
  299 + self.bgd.save(True)
  300 + self.bgd.close()
  301 + self.bgd = None
202 302  
203 303 # Return
204 304 return
205 305  
206 306 def create_cif(self):
207 307 """
208   - Create CIF extension in FITS file.
  308 + Create Calibration Index File (CIF) extension in CIF FITS file.
  309 +
  310 + Parameters:
  311 + None
  312 + Keywords:
  313 + None
209 314 """
210 315 # Create binary table
211 316 table = GFitsBinTable()
... ... @@ -248,6 +353,11 @@ class caldb():
248 353 def add_cif_info(self):
249 354 """
250 355 Add information to CIF extension.
  356 +
  357 + Parameters:
  358 + None
  359 + Keywords:
  360 + None
251 361 """
252 362 # Append 3 rows to CIF extension
253 363 self.hdu_cif.append_rows(3)
... ... @@ -310,20 +420,35 @@ class caldb():
310 420 self.hdu_cif["CAL_CBD"][row, n] = "NONE"
311 421 else:
312 422 self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n]
  423 +
  424 + # Add background information
  425 + row = self.hdu_cif.nrows()-1
  426 + self.hdu_cif["CAL_DIR"][row] = self.bgd_dir
  427 + self.hdu_cif["CAL_FILE"][row] = self.bgd_file
  428 + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name
  429 + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc
  430 + self.hdu_cif["CAL_XNO"][row] = 1
  431 + for n in range(9):
  432 + if n >= len(self.bgd_bounds):
  433 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  434 + else:
  435 + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
313 436  
314 437 # Return
315 438 return
316 439  
317   - def open_ea(self, filename):
  440 + def open_ea(self):
318 441 """
319 442 Open effective area extension. If no extension was found
320 443 then create one. We do not yet set any data as we don't
321 444 know the dimensions of the table yet.
322   - """
323   - # Open FITS file
324   - self.ea = GFits(filename, True)
325 445  
326   - # Get extensions. If they do not exist then create them now
  446 + Parameters:
  447 + None
  448 + Keywords:
  449 + None
  450 + """
  451 + # Get extension. If it does not exist then create it now.
327 452 try:
328 453 self.hdu_ea = self.ea.table("EFFECTIVE AREA")
329 454 except:
... ... @@ -346,16 +471,18 @@ class caldb():
346 471 # Return
347 472 return
348 473  
349   - def open_psf(self, filename):
  474 + def open_psf(self):
350 475 """
351 476 Open point spread function extension. If no extension was
352 477 found then create one. We do not yet set any data as we don't
353 478 know the dimensions of the table yet.
354   - """
355   - # Open FITS file
356   - self.psf = GFits(filename, True)
357 479  
358   - # Get extensions. If they do not exist then create them now
  480 + Parameters:
  481 + None
  482 + Keywords:
  483 + None
  484 + """
  485 + # Get extension. If it does not exist then create it now.
359 486 try:
360 487 self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
361 488 except:
... ... @@ -378,16 +505,18 @@ class caldb():
378 505 # Return
379 506 return
380 507  
381   - def open_edisp(self, filename):
  508 + def open_edisp(self):
382 509 """
383 510 Open energy dispersion extension. If no extension was found
384 511 then create one. We do not yet set any data as we don't
385 512 know the dimensions of the table yet.
386   - """
387   - # Open FITS file
388   - self.edisp = GFits(filename, True)
389 513  
390   - # Get extensions. If they do not exist then create them now
  514 + Parameters:
  515 + None
  516 + Keywords:
  517 + None
  518 + """
  519 + # Get extension. If it does not exist then create it now.
391 520 try:
392 521 self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
393 522 except:
... ... @@ -410,9 +539,50 @@ class caldb():
410 539 # Return
411 540 return
412 541  
  542 + def open_bgd(self):
  543 + """
  544 + Open background extension. If no extension was found
  545 + then create one. We do not yet set any data as we don't
  546 + know the dimensions of the table yet.
  547 +
  548 + Parameters:
  549 + None
  550 + Keywords:
  551 + None
  552 + """
  553 + # Get extension. If it does not exist then create it now.
  554 + try:
  555 + self.hdu_bgd = self.bgd.table("BACKGROUND")
  556 + except:
  557 + # Create binary table
  558 + table = GFitsBinTable()
  559 +
  560 + # Set extension name
  561 + table.extname("BACKGROUND")
  562 +
  563 + # Set OGIP keywords
  564 + self.set_ogip_keywords(table, self.bgd_doc, \
  565 + ["RESPONSE", self.bgd_name])
  566 +
  567 + # Append table
  568 + self.bgd.append(table)
  569 +
  570 + # Get extension
  571 + self.hdu_bgd = self.bgd.table("BACKGROUND")
  572 +
  573 + # Return
  574 + return
  575 +
413 576 def set_ogip_keywords(self, hdu, hdudoc, hduclas):
414 577 """
415 578 Set standard OGIP keywords for extension.
  579 +
  580 + Parameters:
  581 + hdu - Header Data Unit.
  582 + hdudoc - Documentation reference string
  583 + hduclas - Array of HDUCLAS fields
  584 + Keywords:
  585 + None
416 586 """
417 587 # Set keywords
418 588 hdu.card("ORIGIN", "IRAP", "Name of organization making this file")
... ... @@ -433,12 +603,17 @@ class caldb():
433 603 def make_2D(self, array, hdu, name, unit, scale=1.0):
434 604 """
435 605 Make 2D matrix as function of energy and offset angle from a
436   - ROOT 2D histogram.
  606 + ROOT 2D histogram. If the HDU has already the energy and
  607 + offset angle columns, this method will simply add another data
  608 + column.
  609 + If name==None, the method will not append any data column.
437 610  
438 611 Parameters:
439   - array - ROOT 2D histgram.
  612 + array - ROOT 2D histogram.
440 613 hdu - FITS HDU.
441   - Keyword:
  614 + name - Data column name.
  615 + unit - Data unit.
  616 + Keywords:
442 617 scale - Scaling factor for histogram values.
443 618 """
444 619 # Extract energy and offset angle vectors
... ... @@ -485,7 +660,7 @@ class caldb():
485 660 hdu["THETA_HI"][0,ioff] = o_hi
486 661 #
487 662 # "NAME"
488   - if not hdu.hascolumn(name):
  663 + if name != None and not hdu.hascolumn(name):
489 664 hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset))
490 665 hdu[name].unit(unit)
491 666 hdu[name].dim([neng, noffset])
... ... @@ -515,44 +690,217 @@ class caldb():
515 690 Parameters:
516 691 filename - ROOT 2D performance filename.
517 692 """
518   - # Open performance file
  693 + # Open ROOT performance file
519 694 file = TFile(filename)
520 695  
521   - # Allocate ROOT 2D array
522   - array = TH2F()
  696 + # Create effective area
  697 + self.root2ea(file)
  698 +
  699 + # Create point spread function
  700 + self.root2psf(file)
  701 +
  702 + # Create energy dispersion
  703 + self.root2edisp(file)
  704 +
  705 + # Create background
  706 + self.root2bgd(file)
  707 +
  708 + # Return
  709 + return
523 710  
524   - # Make 2D effective 80% area
  711 + def root2ea(self, file):
  712 + """
  713 + Translate ROOT to CALDB effective area extension. The following ROOT
  714 + histograms are used:
  715 +
  716 + EffectiveAreaEtrue_offaxis -> EFFAREA
  717 + EffectiveArea_offaxis -> EFFAREA_RECO
  718 +
  719 + Parameters:
  720 + file - ROOT file.
  721 + Keywords:
  722 + None
  723 + """
  724 + # Continue only if effective area HDU has been opened
525 725 if self.hdu_ea != None:
526   - file.GetObject("EffectiveArea80_offaxis", array)
527   - bounds = self.make_2D(array, self.hdu_ea, "EFFAREA", "m2", scale=1.0/0.8)
  726 +
  727 + # Allocate ROOT 2D arrays
  728 + etrue = TH2F()
  729 + ereco = TH2F()
  730 + file.GetObject("EffectiveAreaEtrue_offaxis", etrue)
  731 + file.GetObject("EffectiveArea_offaxis", ereco)
  732 + #print etrue.GetXaxis().GetBinLowEdge(1), etrue.GetXaxis().GetBinUpEdge(etrue.GetXaxis().GetNbins())
  733 + #print ereco.GetXaxis().GetBinLowEdge(1), ereco.GetXaxis().GetBinUpEdge(ereco.GetXaxis().GetNbins())
  734 +
  735 + # Rebin etrue histogram
  736 + etrue.RebinX(10)
  737 + neng = etrue.GetXaxis().GetNbins()
  738 + noffset = etrue.GetYaxis().GetNbins()
  739 + for ioff in range(noffset):
  740 + for ieng in range(neng):
  741 + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
  742 + etrue.SetBinContent(ieng+1,ioff+1,value)
  743 +
  744 + # Write boundaries (use Ereco boundaries)
  745 + bounds = self.make_2D(ereco, self.hdu_ea, None, "m2")
528 746 for b in bounds:
529 747 self.ea_bounds.append(b)
530 748 self.set_cif_keywords(self.hdu_ea, self.ea_name, \
531 749 self.ea_bounds, self.ea_desc)
532 750  
533   - # Make 2D 68% point spread function
  751 + # EFFAREA
  752 + self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2")
  753 +
  754 + # EFFAREA_RECO
  755 + self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2")
  756 +
  757 + # Return
  758 + return
  759 +
  760 + def root2psf(self, file):
  761 + """
  762 + Translate ROOT to CALDB point spread function extension. The following
  763 + ROOT histograms are used:
  764 +
  765 + 1/(2*pi*SIGMA_1) -> SCALE
  766 + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
  767 + 0.0 -> AMPL_2
  768 + 0.0 -> SIGMA_2
  769 + 0.0 -> AMPL_3
  770 + 0.0 -> SIGMA_3
  771 +
  772 + Parameters:
  773 + file - ROOT file.
  774 + Keywords:
  775 + None
  776 + """
  777 + # Continue only if point spread function HDU has been opened
534 778 if self.hdu_psf != None:
535   - file.GetObject("AngRes_offaxis", array)
536   - bounds = self.make_2D(array, self.hdu_psf, "R68", "deg")
  779 +
  780 + # Allocate ROOT 2D array
  781 + r68 = TH2F()
  782 +
  783 + # Read 68% containment histogram
  784 + file.GetObject("AngRes_offaxis", r68)
  785 + neng = r68.GetXaxis().GetNbins()
  786 + noffset = r68.GetYaxis().GetNbins()
  787 +
  788 + # Compute scale histogram
  789 + scale = r68.Clone()
  790 + for ioff in range(noffset):
  791 + for ieng in range(neng):
  792 + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
  793 + if integral > 0.0:
  794 + value = 1.0 / integral
  795 + else:
  796 + value = 0.0
  797 + scale.SetBinContent(ieng+1,ioff+1,value)
  798 +
  799 + # Set zero histogram
  800 + zero = r68.Clone()
  801 + for ioff in range(noffset):
  802 + for ieng in range(neng):
  803 + zero.SetBinContent(ieng+1,ioff+1,0.0)
  804 +
  805 + # Set boundaries
  806 + bounds = self.make_2D(r68, self.hdu_psf, None, "deg")
537 807 for b in bounds:
538 808 self.psf_bounds.append(b)
539 809 self.set_cif_keywords(self.hdu_psf, self.psf_name, \
540 810 self.psf_bounds, self.psf_desc)
541 811  
542   - # Make 2D 80% point spread function
543   - if self.hdu_psf != None:
544   - file.GetObject("AngRes80_offaxis", array)
545   - self.make_2D(array, self.hdu_psf, "R80", "deg")
  812 + # SCALE
  813 + self.make_2D(scale, self.hdu_psf, "SCALE", "")
  814 +
  815 + # SIGMA_1
  816 + self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg")
546 817  
547   - # Make 2D energy resolution
  818 + # AMPL_2
  819 + self.make_2D(zero, self.hdu_psf, "AMPL_2", "")
  820 +
  821 + # SIGMA_2
  822 + self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg")
  823 +
  824 + # AMPL_3
  825 + self.make_2D(zero, self.hdu_psf, "AMPL_3", "")
  826 +
  827 + # SIGMA_3
  828 + self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg")
  829 +
  830 + # Return
  831 + return
  832 +
  833 + def root2edisp(self, file):
  834 + """
  835 + Translate ROOT to CALDB energy dispersion extension. The following ROOT
  836 + histograms are used:
  837 +
  838 + ERes_offaxis -> ERES
  839 + Ebias_offaxis -> EBIAS
  840 +
  841 + Parameters:
  842 + file - ROOT file.
  843 + Keywords:
  844 + None
  845 + """
  846 + # Continue only if energy dispersion HDU has been opened
548 847 if self.hdu_edisp != None:
549   - file.GetObject("ERes", array)
550   - bounds = self.make_2D(array, self.hdu_edisp, "ERES68", "")
  848 +
  849 + # Allocate ROOT 2D array
  850 + eres = TH2F()
  851 + ebias = TH2F()
  852 + file.GetObject("ERes_offaxis", eres)
  853 + file.GetObject("Ebias_offaxis", ebias)
  854 +
  855 + # Set boundaries
  856 + bounds = self.make_2D(eres, self.hdu_edisp, None, "deg")
551 857 for b in bounds:
552 858 self.edisp_bounds.append(b)
553 859 self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \
554 860 self.edisp_bounds, self.edisp_desc)
555   -
  861 +
  862 + # ERES
  863 + self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue")
  864 +
  865 + # EBIAS
  866 + self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue")
  867 +
  868 + # Return
  869 + return
  870 +
  871 + def root2bgd(self, file):
  872 + """
  873 + Translate ROOT to CALDB background extension. The following ROOT
  874 + histograms are used:
  875 +
  876 + BGRatePerSqDeg_offaxis -> BGD
  877 + BGRatePerSqDeg_offaxis -> BGD_RECO
  878 +
  879 + Parameters:
  880 + file - ROOT file.
  881 + Keywords:
  882 + None
  883 + """
  884 + # Continue only if background HDU has been opened
  885 + if self.hdu_bgd != None:
  886 +
  887 + # Allocate ROOT 2D array
  888 + array = TH2F()
  889 + file.GetObject("BGRatePerSqDeg_offaxis", array)
  890 +
  891 + # Set boundaries
  892 + bounds = self.make_2D(array, self.hdu_bgd, None, "deg")
  893 + for b in bounds:
  894 + self.bgd_bounds.append(b)
  895 + self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
  896 + self.bgd_bounds, self.bgd_desc)
  897 +
  898 + # BGD
  899 + self.make_2D(array, self.hdu_bgd, "BGD", "")
  900 +
  901 + # BGD_RECO
  902 + self.make_2D(array, self.hdu_bgd, "BGD_RECO", "")
  903 +
556 904 # Return
557 905 return
558 906  
... ... @@ -589,22 +937,41 @@ class caldb():
589 937 #==========================#
590 938 if __name__ == '__main__':
591 939 """
592   - This script translates a ROOT 2D performance file in a CALDB compliant
593   - CTA response.
  940 + This script translates a CTA ROOT 2D performance files in a CALDB compliant
  941 + CTA response. The script will create the response for 3 observations to
  942 + illustrate the run-wise file structure.
  943 +
  944 + The CTA ROOT 2D performance files contain the following histograms:
  945 + - DiffSens_offaxis (2D) - Differential sensitivity
  946 + - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?)
  947 + - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment
  948 + - AngRes_offaxis (2D) - Angular resolution 68% containment
  949 + - AngRes80_offaxis (2D) - Angular resolution 80% containment
  950 + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree
  951 + - BGRate_offaxis (2D) - Background rate
  952 + - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins)
  953 + - MigMatrix_offaxis (3D)
  954 + - EestOverEtrue_offaxis (3D)
  955 + - ERes_offaxis (2D) - Energy Resolution
  956 + - Ebias_offaxis (2D) - Energy bias
594 957 """
595 958 # Set ROOT filename
596   - #path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxis"
597   - #irf = "SubarrayE_offaxis.root"
598   - path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_Nov2011"
599   - irf = "SubarrayE_IFAE_50hours_20111121_offaxis.root"
  959 + path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_May2012"
  960 + irf = "SubarrayE_IFAE_50hours_20120510_offaxis.root"
600 961 filename = path+"/"+irf
601 962  
602   - # Allocate caldb
603   - db = caldb("E")
  963 + # Set observation identifiers
  964 + obsids = ["000001", "000002", "000003"]
  965 +
  966 + # Generate CALDB entries, one for each observation (or run)
  967 + for obsid in obsids:
  968 +
  969 + # Allocate caldb
  970 + db = caldb("E", obsid)
604 971  
605   - # Open calibration
606   - db.open("test")
  972 + # Open calibration
  973 + db.open("test")
607 974  
608   - # Translate ROOT to CALDB information
609   - db.root2caldb(filename)
  975 + # Translate ROOT to CALDB information
  976 + db.root2caldb(filename)
610 977  
611 978 \ No newline at end of file
... ...
scripts/obsutils.py
  View file @b930a25
... ... @@ -183,7 +183,7 @@ def cntmap(obs, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
183 183 event = cast_GCTAEventAtom(event)
184 184  
185 185 # Determine sky pixel
186   - skydir = event.dir().skydir()
  186 + skydir = event.dir().dir()
187 187 pixel = map.dir2pix(skydir)
188 188  
189 189 # Set pixel
... ... @@ -237,7 +237,7 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
237 237  
238 238 # Get sky direction
239 239 skydir = map.pix2dir(pixel)
240   - instdir.skydir(skydir)
  240 + instdir.dir(skydir)
241 241  
242 242 # Create event atom for map pixel
243 243 atom = GCTAEventAtom()
... ... @@ -295,7 +295,7 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
295 295 # Set ROI
296 296 roi = GCTARoi()
297 297 instdir = GCTAInstDir()
298   - instdir.skydir(pntdir)
  298 + instdir.dir(pntdir)
299 299 roi.centre(instdir)
300 300 roi.radius(rad)
301 301  
... ...
src/Makefile.am
  View file @b930a25
... ... @@ -5,10 +5,11 @@ MACOSX_DEPLOYMENT_TARGET=@CONFIGURE_MACOSX_DEPLOYMENT_TARGET@
5 5 @EXPORT_MACOSX_DEPLOYMENT_TARGET@export MACOSX_DEPLOYMENT_TARGET
6 6  
7 7 # Set sub directories
8   -SUBDIRS = ctselect \
9   - ctbin \
  8 +SUBDIRS = ctbin \
10 9 ctlike \
  10 + ctmodel \
11 11 ctobssim \
  12 + ctselect \
12 13 ctskymap
13 14  
14 15 # Scripts to be distributed
... ... @@ -22,9 +23,10 @@ scripts = ctools-setup \
22 23 lib_LTLIBRARIES = libctools.la
23 24 libctools_la_SOURCES =
24 25 libctools_la_LIBADD = ctbin/libctbin.la \
25   - ctselect/libctselect.la \
26 26 ctlike/libctlike.la \
  27 + ctmodel/libctmodel.la \
27 28 ctobssim/libctobssim.la \
  29 + ctselect/libctselect.la \
28 30 ctskymap/libctskymap.la
29 31  
30 32 # Distribute scripts
... ...
src/ctbin/ctbin.cpp
  View file @b930a25
... ... @@ -534,7 +534,7 @@ void ctbin::bin_events(GCTAObservation* obs)
534 534  
535 535 // Determine sky pixel
536 536 GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
537   - GSkyDir dir = inst->skydir();
  537 + GSkyDir dir = inst->dir();
538 538 GSkyPixel pixel = map.dir2xy(dir);
539 539  
540 540 // Skip if pixel is out of range
... ...
src/ctmodel/Makefile.am 0 → 100644
  View file @b930a25
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Program name
  4 +bin_PROGRAMS = ctmodel
  5 +
  6 +# Sources
  7 +sources = ctmodel.cpp
  8 +
  9 +# Libtool library
  10 +noinst_LTLIBRARIES = libctmodel.la
  11 +libctmodel_la_SOURCES = $(sources)
  12 +
  13 +# Program definition
  14 +ctmodel_SOURCES = main.cpp
  15 +ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  16 +ctmodel_CPPFLAGS = @CPPFLAGS@
  17 +ctmodel_LDADD = libctmodel.la
  18 +
  19 +# Parameter files
  20 +pfilesdir = $(prefix)/syspfiles
  21 +pfiles_DATA = ctmodel.par
... ...
src/ctmodel/ctmodel.cpp 0 → 100644
  View file @b930a25
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctmodel.cpp
  23 + * @brief CTA counts model tool implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctmodel.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_SETUP_OBS "ctmodel::setup_obs()"
  37 +#define G_MODEL_MAP "ctmodel::model_map(GCTAObservation*)"
  38 +
  39 +/* __ Debug definitions __________________________________________________ */
  40 +
  41 +/* __ Coding definitions _________________________________________________ */
  42 +
  43 +
  44 +/*==========================================================================
  45 + = =
  46 + = Constructors/destructors =
  47 + = =
  48 + ==========================================================================*/
  49 +
  50 +/***********************************************************************//**
  51 + * @brief Void constructor
  52 + ***************************************************************************/
  53 +ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
  54 +{
  55 + // Initialise members
  56 + init_members();
  57 +
  58 + // Write header into logger
  59 + log_header();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Observations constructor
  68 + *
  69 + * This method creates an instance of the class by copying an existing
  70 + * observations container.
  71 + ***************************************************************************/
  72 +ctmodel::ctmodel(GObservations obs) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Set observations
  78 + m_obs = obs;
  79 +
  80 + // Write header into logger
  81 + log_header();
  82 +
  83 + // Return
  84 + return;
  85 +}
  86 +
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Command line constructor
  91 + *
  92 + * @param[in] argc Number of arguments in command line.
  93 + * @param[in] argv Array of command line arguments.
  94 + ***************************************************************************/
  95 +ctmodel::ctmodel(int argc, char *argv[]) :
  96 + GApplication(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
  97 +{
  98 + // Initialise members
  99 + init_members();
  100 +
  101 + // Write header into logger
  102 + log_header();
  103 +
  104 + // Return
  105 + return;
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Copy constructor
  111 + *
  112 + * @param[in] app Application.
  113 + ***************************************************************************/
  114 +ctmodel::ctmodel(const ctmodel& app) : GApplication(app)
  115 +{
  116 + // Initialise members
  117 + init_members();
  118 +
  119 + // Copy members
  120 + copy_members(app);
  121 +
  122 + // Return
  123 + return;
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Destructor
  129 + ***************************************************************************/
  130 +ctmodel::~ctmodel(void)
  131 +{
  132 + // Free members
  133 + free_members();
  134 +
  135 + // Return
  136 + return;
  137 +}
  138 +
  139 +
  140 +/*==========================================================================
  141 + = =
  142 + = Operators =
  143 + = =
  144 + ==========================================================================*/
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Assignment operator
  148 + *
  149 + * @param[in] app ctmodel application.
  150 + * @return Returns ctmodel application.
  151 + ***************************************************************************/
  152 +ctmodel& ctmodel::operator= (const ctmodel& app)
  153 +{
  154 + // Execute only if object is not identical
  155 + if (this != &app) {
  156 +
  157 + // Copy base class members
  158 + this->GApplication::operator=(app);
  159 +
  160 + // Free members
  161 + free_members();
  162 +
  163 + // Initialise members
  164 + init_members();
  165 +
  166 + // Copy members
  167 + copy_members(app);
  168 +
  169 + } // endif: object was not identical
  170 +
  171 + // Return this object
  172 + return *this;
  173 +}
  174 +
  175 +
  176 +/*==========================================================================
  177 + = =
  178 + = Public methods =
  179 + = =
  180 + ==========================================================================*/
  181 +
  182 +/***********************************************************************//**
  183 + * @brief Clear instance
  184 + ***************************************************************************/
  185 +void ctmodel::clear(void)
  186 +{
  187 + // Free members
  188 + free_members();
  189 + this->GApplication::free_members();
  190 +
  191 + // Initialise members
  192 + this->GApplication::init_members();
  193 + init_members();
  194 +
  195 + // Return
  196 + return;
  197 +}
  198 +
  199 +
  200 +/***********************************************************************//**
  201 + * @brief Execute application
  202 + *
  203 + * This is the main execution method of the ctmodel class. It is invoked
  204 + * when the executable is called from command line. The method generates
  205 + * the model maps and saves the results.
  206 + ***************************************************************************/
  207 +void ctmodel::execute(void)
  208 +{
  209 + // Signal that some parameters should be read ahead
  210 + m_read_ahead = true;
  211 +
  212 + // Create the model map(s)
  213 + run();
  214 +
  215 + // Save the model map(s) into FITS file
  216 + save();
  217 +
  218 + // Return
  219 + return;
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
  224 + * @brief Generate the model map(s)
  225 + *
  226 + * This method reads the task parameters from the parfile, sets up the
  227 + * observation container, loops over all CTA observations in the container
  228 + * and generates a model map for each CTA observation.
  229 + ***************************************************************************/
  230 +void ctmodel::run(void)
  231 +{
  232 + // If we're in debug mode then all output is also dumped on the screen
  233 + if (logDebug()) {
  234 + log.cout(true);
  235 + }
  236 +
  237 + // Get task parameters
  238 + get_parameters();
  239 +
  240 + // Setup observation container
  241 + setup_obs();
  242 +
  243 + // Write parameters into logger
  244 + if (logTerse()) {
  245 + log_parameters();
  246 + log << std::endl;
  247 + }
  248 +
  249 + // Write observation(s) into logger
  250 + if (logTerse()) {
  251 + log << std::endl;
  252 + if (m_obs.size() > 1) {
  253 + log.header1("Observations");
  254 + }
  255 + else {
  256 + log.header1("Observation");
  257 + }
  258 + log << m_obs << std::endl;
  259 + }
  260 +
  261 + // Write header
  262 + if (logTerse()) {
  263 + log << std::endl;
  264 + if (m_obs.size() > 1) {
  265 + log.header1("Generate model maps");
  266 + }
  267 + else {
  268 + log.header1("Generate model map");
  269 + }
  270 + }
  271 +
  272 + // Initialise observation counter
  273 + int n_observations = 0;
  274 +
  275 + // Loop over all observations in the container
  276 + for (int i = 0; i < m_obs.size(); ++i) {
  277 +
  278 + // Initialise event input and output filenames
  279 + m_infiles.push_back("");
  280 +
  281 + // Get CTA observation
  282 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  283 +
  284 + // Continue only if observation is a CTA observation
  285 + if (obs != NULL) {
  286 +
  287 + // Write header for observation
  288 + if (logTerse()) {
  289 + if (obs->name().length() > 1) {
  290 + log.header3("Observation "+obs->name());
  291 + }
  292 + else {
  293 + log.header3("Observation");
  294 + }
  295 + }
  296 +
  297 + // Increment number of observations
  298 + n_observations++;
  299 +
  300 + // Save event file name (for possible saving)
  301 + m_infiles[i] = obs->eventfile();
  302 +
  303 + // Generate model map
  304 + model_map(obs, m_obs.models());
  305 +
  306 + } // endif: CTA observation found
  307 +
  308 + } // endfor: looped over observations
  309 +
  310 + // If more than a single observation has been handled then make sure
  311 + // that an XML file will be used for storage
  312 + if (n_observations > 1) {
  313 + m_use_xml = true;
  314 + }
  315 +
  316 + // Write observation(s) into logger
  317 + if (logTerse()) {
  318 + log << std::endl;
  319 + if (m_obs.size() > 1) {
  320 + log.header1("Observations after model map generation");
  321 + }
  322 + else {
  323 + log.header1("Observation after model map generation");
  324 + }
  325 + log << m_obs << std::endl;
  326 + }
  327 +
  328 + // Return
  329 + return;
  330 +}
  331 +
  332 +
  333 +/***********************************************************************//**
  334 + * @brief Save model map(s)
  335 + *
  336 + * This method saves the model map(s) into FITS file(s). There are two
  337 + * modes, depending on the m_use_xml flag.
  338 + *
  339 + * If m_use_xml is true, all model map(s) will be saved into FITS files,
  340 + * where the output filenames are constructued from the input filenames by
  341 + * prepending the m_prefix string to name. Any path information will be
  342 + * stripped form the input name, hence event files will be written into the
  343 + * local working directory (unless some path information is present
  344 + * in the prefix). In addition, an XML file will be created that gathers
  345 + * the filename information for the model map(s). If an XML file was present
  346 + * on input, all metadata information will be copied from this input file.
  347 + *
  348 + * If m_use_xml is false, the model map will be saved into a FITS file.
  349 + ***************************************************************************/
  350 +void ctmodel::save(void)
  351 +{
  352 + // Write header
  353 + if (logTerse()) {
  354 + log << std::endl;
  355 + if (m_obs.size() > 1) {
  356 + log.header1("Save observations");
  357 + }
  358 + else {
  359 + log.header1("Save observation");
  360 + }
  361 + }
  362 +
  363 + // Case A: Save model map(s) and XML metadata information
  364 + if (m_use_xml) {
  365 + save_xml();
  366 + }
  367 +
  368 + // Case B: Save model map as FITS file
  369 + else {
  370 + save_fits();
  371 + }
  372 +
  373 + // Return
  374 + return;
  375 +}
  376 +
  377 +
  378 +/***********************************************************************//**
  379 + * @brief Get application parameters
  380 + *
  381 + * Get all task parameters from parameter file or (if required) by querying
  382 + * the user. The parameters are read in the correct order.
  383 + ***************************************************************************/
  384 +void ctmodel::get_parameters(void)
  385 +{
  386 + // If we do not have any observations in the container then get all
  387 + // information to build (at least) one
  388 + if (m_obs.size() == 0) {
  389 +
  390 + // Read general parameters. Only read the output parameters if
  391 + // required, and also only read the model XML filename if no
  392 + // model is yet in the container.
  393 + m_infile = (*this)["infile"].filename();
  394 + if (m_read_ahead) {
  395 + m_outfile = (*this)["outfile"].filename();
  396 + m_prefix = (*this)["prefix"].string();
  397 + }
  398 + m_caldb = (*this)["caldb"].string();
  399 + m_irf = (*this)["irf"].string();
  400 + if (m_obs.models().size() == 0) {
  401 + m_srcmdl = (*this)["srcmdl"].filename();
  402 + }
  403 +
  404 + // If there is no input filename then read all parameters that
  405 + // are required to build a model map from scratch
  406 + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) {
  407 + m_ra = (*this)["ra"].real();
  408 + m_dec = (*this)["dec"].real();
  409 + m_deadc = (*this)["deadc"].real();
  410 + m_tmin = (*this)["tmin"].real();
  411 + m_tmax = (*this)["tmax"].real();
  412 + m_emin = (*this)["emin"].real();
  413 + m_emax = (*this)["emax"].real();
  414 + m_enumbins = (*this)["enumbins"].integer();
  415 + m_proj = (*this)["proj"].string();
  416 + m_coordsys = (*this)["coordsys"].string();
  417 + m_xref = (*this)["xref"].real();
  418 + m_yref = (*this)["yref"].real();
  419 + m_binsz = (*this)["binsz"].real();
  420 + m_nxpix = (*this)["nxpix"].integer();
  421 + m_nypix = (*this)["nypix"].integer();
  422 + }
  423 + }
  424 +
  425 + // ... otherwise, read only the parameters that are required for the
  426 + // model generation.
  427 + else {
  428 + if (m_read_ahead) {
  429 + m_outfile = (*this)["outfile"].filename();
  430 + m_prefix = (*this)["prefix"].string();
  431 + }
  432 + m_caldb = (*this)["caldb"].string();
  433 + m_irf = (*this)["irf"].string();
  434 + if (m_obs.models().size() == 0) {
  435 + m_srcmdl = (*this)["srcmdl"].filename();
  436 + }
  437 + }
  438 +
  439 + // Return
  440 + return;
  441 +}
  442 +
  443 +
  444 +/***********************************************************************//**
  445 + * @brief Setup observation container
  446 + *
  447 + * @exception GException::no_cube
  448 + * No event cube found in CTA observation.
  449 + *
  450 + * This method sets up the observation container for processing. There are
  451 + * two cases:
  452 + *
  453 + * If there are no observations in the actual observation container, the
  454 + * method will check in "infile" parameter. If this parameter is "NONE" or
  455 + * empty, the task parameters will be used to construct a model map.
  456 + * Otherwise, the method first tries to interpret the "infile" parameter as
  457 + * a counts map, and attemps loading of the file in an event cube. If this
  458 + * fails, the method tries to interpret the "infile" parameter as an
  459 + * observation definition XML file. If this also fails, an exception will
  460 + * be thrown.
  461 + *
  462 + * If observations exist already in the observation container, the method
  463 + * will simply keep them.
  464 + *
  465 + * Test if all CTA observations contain counts maps.
  466 + *
  467 + * Finally, if no models exist so far in the observation container, the
  468 + * models will be loaded from the model XML file.
  469 + ***************************************************************************/
  470 +void ctmodel::setup_obs(void)
  471 +{
  472 + // If there are no observations in the container then try to build some
  473 + if (m_obs.size() == 0) {
  474 +
  475 + // If no input filename has been specified, then create a model map
  476 + // from the task parameters
  477 + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) {
  478 +
  479 + // Set pointing direction
  480 + GCTAPointing pnt;
  481 + GSkyDir skydir;
  482 + skydir.radec_deg(m_ra, m_dec);
  483 + pnt.dir(skydir);
  484 +
  485 + // Setup energy range covered by model
  486 + GEnergy emin;
  487 + GEnergy emax;
  488 + GEbounds ebds;
  489 + emin.TeV(m_emin);
  490 + emax.TeV(m_emax);
  491 + ebds.setlog(emin, emax, m_enumbins);
  492 +
  493 + // Setup time interval covered by model
  494 + GGti gti;
  495 + GTime tmin;
  496 + GTime tmax;
  497 + tmin.met(m_tmin);
  498 + tmax.met(m_tmax);
  499 + gti.append(tmin, tmax);
  500 +
  501 + // Setup skymap
  502 + GSkymap map = GSkymap(m_proj, m_coordsys,
  503 + m_xref, m_yref, m_binsz, m_binsz,
  504 + m_nxpix, m_nypix, m_enumbins);
  505 +
  506 + // Create model cube from sky map
  507 + GCTAEventCube cube(map, ebds, gti);
  508 +
  509 + // Allocate CTA observation
  510 + GCTAObservation obs;
  511 +
  512 + // Set CTA observation attributes
  513 + obs.pointing(pnt);
  514 + obs.ontime(gti.ontime());
  515 + obs.livetime(gti.ontime()*m_deadc);
  516 + obs.deadc(m_deadc);
  517 +
  518 + // Set event cube in observation
  519 + obs.events(&cube);
  520 +
  521 + // Append CTA observation to container
  522 + m_obs.append(obs);
  523 +
  524 + // Signal that no XML file should be used for storage
  525 + m_use_xml = false;
  526 +
  527 + } // endif: created model map from task parameters
  528 +
  529 + // ... otherwise try to load information from the file
  530 + else {
  531 +
  532 + // First try to open the file as a counts map
  533 + try {
  534 +
  535 + // Allocate CTA observation
  536 + GCTAObservation obs;
  537 +
  538 + // Load counts map in CTA observation
  539 + obs.load_binned(m_infile);
  540 +
  541 + // Append CTA observation to container
  542 + m_obs.append(obs);
  543 +
  544 + // Signal that no XML file should be used for storage
  545 + m_use_xml = false;
  546 +
  547 + }
  548 +
  549 + // ... otherwise try to open as XML file
  550 + catch (GException::fits_open_error &e) {
  551 +
  552 + // Load observations from XML file. This will throw
  553 + // an exception if it fails.
  554 + m_obs.load(m_infile);
  555 +
  556 + // Signal that XML file should be used for storage
  557 + m_use_xml = true;
  558 +
  559 + }
  560 +
  561 + } // endelse: loaded information from input file
  562 +
  563 + } // endif: there was no observation in the container
  564 +
  565 + // If there are no models associated with the observations then
  566 + // load now the model definition from the XML file
  567 + if (m_obs.models().size() == 0) {
  568 + m_obs.models(GModels(m_srcmdl));
  569 + }
  570 +
  571 + // Check if all CTA observations contain an event cube and setup response
  572 + // for all observations
  573 + for (int i = 0; i < m_obs.size(); ++i) {
  574 +
  575 + // Is this observation a CTA observation?
  576 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  577 +
  578 + // Yes ...
  579 + if (obs != NULL) {
  580 +
  581 + // Throw an exception if this observation does not contain
  582 + // an event cube
  583 + if (dynamic_cast<const GCTAEventCube*>(obs->events()) == NULL) {
  584 + throw GException::no_cube(G_SETUP_OBS);
  585 + }
  586 +
  587 + // Set response
  588 + obs->response(m_irf, m_caldb);
  589 +
  590 + } // endif: observation was a CTA observation
  591 + } // endfor: looped over all observations
  592 +
  593 + // Return
  594 + return;
  595 +}
  596 +
  597 +
  598 +/***********************************************************************//**
  599 + * @brief Generate model map
  600 + *
  601 + * @param[in] obs CTA observation pointer.
  602 + * @param[in] models Model container.
  603 + *
  604 + * @exception GException::no_cube
  605 + * No event cube found in CTA observation.
  606 + ***************************************************************************/
  607 +void ctmodel::model_map(GCTAObservation* obs, const GModels& models)
  608 +{
  609 + // Continue only if observation pointer is valid
  610 + if (obs != NULL) {
  611 +
  612 + // Get event cube pointer
  613 + GCTAEventCube* cube =
  614 + const_cast<GCTAEventCube*>(dynamic_cast<const GCTAEventCube*>(obs->events()));
  615 +
  616 + // Throw an exception if the observation does not hold and event
  617 + // cube
  618 + if (cube == NULL) {
  619 + throw GException::no_cube(G_MODEL_MAP);
  620 + }
  621 +
  622 + // Initialise statistics
  623 + double sum = 0.0;
  624 +
  625 + // Loop over all events in counts map
  626 + for (int i = 0; i < cube->size(); ++i) {
  627 +
  628 + // Get event bin
  629 + GCTAEventBin* bin = (*cube)[i];
  630 +
  631 + // Compute model value for event bin
  632 + double model =
  633 + models.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
  634 + bin->size();
  635 +
  636 + // Store value
  637 + bin->counts(model);
  638 +
  639 + // Sum all events
  640 + sum += model;
  641 + }
  642 +
  643 + // Log results
  644 + if (logTerse()) {
  645 + log << parformat("Model events in cube");
  646 + log << sum << std::endl;
  647 + }
  648 +
  649 + } // endif: observation pointer was not valid
  650 +
  651 + // Return
  652 + return;
  653 +}
  654 +
  655 +
  656 +/*==========================================================================
  657 + = =
  658 + = Private methods =
  659 + = =
  660 + ==========================================================================*/
  661 +
  662 +/***********************************************************************//**
  663 + * @brief Initialise class members
  664 + ***************************************************************************/
  665 +void ctmodel::init_members(void)
  666 +{
  667 + // Initialise members
  668 + m_infile.clear();
  669 + m_outfile.clear();
  670 + m_prefix.clear();
  671 + m_caldb.clear();
  672 + m_irf.clear();
  673 + m_srcmdl.clear();
  674 + m_proj.clear();
  675 + m_coordsys.clear();
  676 + m_ra = 0.0;
  677 + m_dec = 0.0;
  678 + m_deadc = 1.0;
  679 + m_tmin = 0.0;
  680 + m_tmax = 0.0;
  681 + m_emin = 0.0;
  682 + m_emax = 0.0;
  683 + m_enumbins = 0;
  684 + m_xref = 0.0;
  685 + m_yref = 0.0;
  686 + m_binsz = 0.0;
  687 + m_nxpix = 0;
  688 + m_nypix = 0;
  689 +
  690 + // Initialise protected members
  691 + m_obs.clear();
  692 + m_infiles.clear();
  693 + m_use_xml = false;
  694 + m_read_ahead = false;
  695 +
  696 + // Set logger properties
  697 + log.date(true);
  698 +
  699 + // Return
  700 + return;
  701 +}
  702 +
  703 +
  704 +/***********************************************************************//**
  705 + * @brief Copy class members
  706 + *
  707 + * @param[in] app Application.
  708 + ***************************************************************************/
  709 +void ctmodel::copy_members(const ctmodel& app)
  710 +{
  711 + // Copy attributes
  712 + m_infile = app.m_infile;
  713 + m_outfile = app.m_outfile;
  714 + m_prefix = app.m_prefix;
  715 + m_caldb = app.m_caldb;
  716 + m_irf = app.m_irf;
  717 + m_srcmdl = app.m_srcmdl;
  718 + m_proj = app.m_proj;
  719 + m_coordsys = app.m_coordsys;
  720 + m_ra = app.m_ra;
  721 + m_dec = app.m_dec;
  722 + m_deadc = app.m_deadc;
  723 + m_tmin = app.m_tmin;
  724 + m_tmax = app.m_tmax;
  725 + m_emin = app.m_emin;
  726 + m_emax = app.m_emax;
  727 + m_enumbins = app.m_enumbins;
  728 + m_xref = app.m_xref;
  729 + m_yref = app.m_yref;
  730 + m_binsz = app.m_binsz;
  731 + m_nxpix = app.m_nxpix;
  732 + m_nypix = app.m_nypix;
  733 +
  734 + // Copy protected members
  735 + m_obs = app.m_obs;
  736 + m_infiles = app.m_infiles;
  737 + m_use_xml = app.m_use_xml;
  738 + m_read_ahead = app.m_read_ahead;
  739 +
  740 + // Return
  741 + return;
  742 +}
  743 +
  744 +
  745 +/***********************************************************************//**
  746 + * @brief Delete class members
  747 + ***************************************************************************/
  748 +void ctmodel::free_members(void)
  749 +{
  750 + // Write separator into logger
  751 + if (logTerse()) {
  752 + log << std::endl;
  753 + }
  754 +
  755 + // Return
  756 + return;
  757 +}
  758 +
  759 +
  760 +/***********************************************************************//**
  761 + * @brief Set output file name.
  762 + *
  763 + * @param[in] filename Input file name.
  764 + *
  765 + * This method converts an input filename into an output filename by
  766 + * prepending the prefix specified by m_prefix to the input filename. Any
  767 + * path will be stripped from the input filename.
  768 + ***************************************************************************/
  769 +std::string ctmodel::set_outfile_name(const std::string& filename) const
  770 +{
  771 + // Split input filename into path elements
  772 + std::vector<std::string> elements = split(filename, "/");
  773 +
  774 + // The last path element is the filename
  775 + std::string outname = m_prefix + elements[elements.size()-1];
  776 +
  777 + // Return output filename
  778 + return outname;
  779 +}
  780 +
  781 +
  782 +/***********************************************************************//**
  783 + * @brief Save model map in FITS format.
  784 + *
  785 + * Save the model map as a FITS file. The filename of the FITS file is
  786 + * specified by the m_outfile member.
  787 + ***************************************************************************/
  788 +void ctmodel::save_fits(void)
  789 +{
  790 + // Get output filename
  791 + m_outfile = (*this)["outfile"].filename();
  792 +
  793 + // Get CTA observation from observation container
  794 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[0]);
  795 +
  796 + // Save model map
  797 + save_model_map(obs, m_outfile);
  798 +
  799 + // Return
  800 + return;
  801 +}
  802 +
  803 +
  804 +/***********************************************************************//**
  805 + * @brief Save model map(s) in XML format.
  806 + *
  807 + * Save the model map(s) into FITS files and write the file path information
  808 + * into a XML file. The filename of the XML file is specified by the
  809 + * m_outfile member, the filename(s) of the model map(s) are built by
  810 + * prepending the prefix given by the m_prefix member to the input model
  811 + * map(s) filenames. Any path present in the input filename will be stripped,
  812 + * i.e. the model map(s) will be written in the local working directory
  813 + * (unless a path is specified in the m_prefix member).
  814 + ***************************************************************************/
  815 +void ctmodel::save_xml(void)
  816 +{
  817 + // Get output filename and prefix
  818 + m_outfile = (*this)["outfile"].filename();
  819 + m_prefix = (*this)["prefix"].string();
  820 +
  821 + // Loop over all observation in the container
  822 + for (int i = 0; i < m_obs.size(); ++i) {
  823 +
  824 + // Get CTA observation
  825 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(&m_obs[i]);
  826 +
  827 + // Handle only CTA observations
  828 + if (obs != NULL) {
  829 +
  830 + // Set event output file name
  831 + std::string outfile = set_outfile_name(m_infiles[i]);
  832 +
  833 + // Store output file name in observation
  834 + obs->eventfile(outfile);
  835 +
  836 + // Save event list
  837 + save_model_map(obs, outfile);
  838 +
  839 + } // endif: observation was a CTA observations
  840 +
  841 + } // endfor: looped over observations
  842 +
  843 + // Save observations in XML file
  844 + m_obs.save(m_outfile);
  845 +
  846 + // Return
  847 + return;
  848 +}
  849 +
  850 +
  851 +/***********************************************************************//**
  852 + * @brief Save a single model map into a FITS file
  853 + *
  854 + * @param[in] obs Pointer to CTA observation.
  855 + * @param[in] outfile Output file name.
  856 + *
  857 + * This method saves a single model map into a FITS file. The method does
  858 + * nothing if the observation pointer is not valid.
  859 + ***************************************************************************/
  860 +void ctmodel::save_model_map(const GCTAObservation* obs,
  861 + const std::string& outfile) const
  862 +{
  863 + // Save only if observation is valid
  864 + if (obs != NULL) {
  865 +
  866 + // Save observation into FITS file
  867 + obs->save(outfile, clobber());
  868 +
  869 + } // endif: observation was valid
  870 +
  871 + // Return
  872 + return;
  873 +}
... ...
src/ctmodel/ctmodel.hpp 0 → 100644
  View file @b930a25
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctmodel.hpp
  23 + * @brief CTA counts model tool definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef CTMODEL_HPP
  28 +#define CTMODEL_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +
  34 +/* __Definitions _________________________________________________________ */
  35 +#define CTMODEL_NAME "ctmodel"
  36 +#define CTMODEL_VERSION "00-01-00"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class ctmodel
  41 + *
  42 + * @brief CTA counts model tool interface defintion
  43 + *
  44 + * This class creates counts model map(s). The definition of the counts model
  45 + * can be taken from a predefined observation container, from the counts maps
  46 + * found in an observation definition XML file, from an individual counts
  47 + * map, or from the parameter definition.
  48 + *
  49 + * Results are stored in an observation container that can be written to disk
  50 + * in form of FITS files (model maps) and an updated observation definition
  51 + * XML file.
  52 + ***************************************************************************/
  53 +class ctmodel : public GApplication {
  54 +
  55 +public:
  56 + // Constructors and destructors
  57 + ctmodel(void);
  58 + explicit ctmodel(GObservations obs);
  59 + ctmodel(int argc, char *argv[]);
  60 + ctmodel(const ctmodel& app);
  61 + virtual ~ctmodel(void);
  62 +
  63 + // Operators
  64 + ctmodel& operator= (const ctmodel& app);
  65 +
  66 + // Methods
  67 + void clear(void);
  68 + void execute(void);
  69 + void run(void);
  70 + void save(void);
  71 + GObservations& obs(void) { return m_obs; }
  72 + void get_parameters(void);
  73 + void setup_obs(void);
  74 + void model_map(GCTAObservation* obs, const GModels& models);
  75 +
  76 +protected:
  77 + // Protected methods
  78 + void init_members(void);
  79 + void copy_members(const ctmodel& app);
  80 + void free_members(void);
  81 + std::string set_outfile_name(const std::string& filename) const;
  82 + void save_fits(void);
  83 + void save_xml(void);
  84 + void save_model_map(const GCTAObservation* obs,
  85 + const std::string& outfile) const;
  86 +
  87 + // User parameters
  88 + std::string m_infile; //!< Input counts map or XML file
  89 + std::string m_outfile; //!< Output model map or XML file
  90 + std::string m_prefix; //!< Prefix for multiple model maps
  91 + std::string m_caldb; //!< Calibration database
  92 + std::string m_irf; //!< Instrument response functions
  93 + std::string m_srcmdl; //!< Source model
  94 + double m_ra; //!< RA of pointing direction
  95 + double m_dec; //!< DEC of pointing direction
  96 + double m_deadc; //!< Deadtime correction factor
  97 + double m_tmin; //!< Start time
  98 + double m_tmax; //!< Stop time
  99 + double m_emin; //!< Lower energy
  100 + double m_emax; //!< Upper energy
  101 + int m_enumbins; //!< Number of energy bins
  102 + std::string m_proj; //!< WCS projection
  103 + std::string m_coordsys; //!< Coordinate system
  104 + double m_xref; //!< Longitude reference coordinate
  105 + double m_yref; //!< Latitude reference coordinate
  106 + double m_binsz; //!< Pixel size
  107 + int m_nxpix; //!< Number of pixels in longitude
  108 + int m_nypix; //!< Number of pixels in latitude
  109 +
  110 + // Protected members
  111 + GObservations m_obs; //!< Observation container
  112 + std::vector<std::string> m_infiles; //!< Input event filenames
  113 + bool m_use_xml; //!< Use XML file instead of FITS file
  114 + bool m_read_ahead; //!< Read ahead parameters
  115 +};
  116 +
  117 +#endif /* CTMODEL_HPP */
... ...
src/ctmodel/ctmodel.par 0 → 100644
  View file @b930a25
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctmodel.par #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Juergen Knoedlseder #
  11 +# Jurgen.Knodlseder@irap.omp.eu #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the ctmodel application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +infile, f, a, "NONE",,, "Input counts map or observation definition file"
  21 +outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file"
  22 +prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file"
  23 +caldb, s, a, "test/irf",,, "Calibration database"
  24 +irf, s, a, "kb_E_50h_v3",,, "Instrument response function"
  25 +srcmdl, f, a, "test/data/crab.xml",,, "Source model"
  26 +#
  27 +# Instrument parameters
  28 +#======================
  29 +ra, r, a, 83.63, 0,360, "RA of pointing (degrees)"
  30 +dec, r, a, 22.01, -90, 90, "Dec of pointing (degrees)"
  31 +deadc, r, h, 0.95, 0, 1, "Average deadtime correction factor"
  32 +#
  33 +# Time range parameters
  34 +#======================
  35 +tmin, r, a, 0.0, 0, , "Start time (MET in s)"
  36 +tmax, r, a, 1800.0, 0, , "End time (MET in s)"
  37 +#
  38 +# Energy binning parameters
  39 +#==========================
  40 +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins"
  41 +emin, r, a, 0.1,,,"Start value for first energy bin in TeV"
  42 +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV"
  43 +enumbins, i, a, 20,,,"Number of energy bins"
  44 +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
  45 +#
  46 +# Spatial binning parameters
  47 +#===========================
  48 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  49 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  50 +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
  51 +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  52 +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)"
  53 +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)"
  54 +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees"
  55 +proj, s, a, "CAR",,, "Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN"
  56 +#
  57 +# Standard parameters
  58 +#====================
  59 +chatter, i, h, 1,0,4, "Chattiness of output"
  60 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  61 +debug, b, h, no,,, "Debugging mode activated"
  62 +mode, s, h, "ql",,, "Mode of automatic parameters"
... ...
src/ctmodel/main.cpp 0 → 100644
  View file @b930a25
  1 +/***************************************************************************
  2 + * ctmodel - CTA counts model tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file main.cpp
  23 + * @brief CTA counts model tool main code
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctmodel.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point of application
  37 + *
  38 + * @param[in] argc Number of command line arguments.
  39 + * @param[in] argv Command line arguments.
  40 + *
  41 + * This is the main entry point of the gtbin application. It allocates a
  42 + * ctmodel object and runs the application.
  43 + ***************************************************************************/
  44 +int main (int argc, char *argv[])
  45 +{
  46 + // Initialise return code
  47 + int rc = 1;
  48 +
  49 + // Create instance of application
  50 + ctmodel application(argc, argv);
  51 +
  52 + // Run application
  53 + try {
  54 + // Execute application
  55 + application.execute();
  56 +
  57 + // Signal success
  58 + rc = 0;
  59 + }
  60 + catch (std::exception &e) {
  61 +
  62 + // Extract error message
  63 + std::string message = e.what();
  64 + std::string signal = "*** ERROR encounterted in the execution of"
  65 + " ctmodel. Run aborted ...";
  66 +
  67 + // Write error in logger
  68 + application.log << signal << std::endl;
  69 + application.log << message << std::endl;
  70 +
  71 + // Write error on standard output
  72 + std::cout << signal << std::endl;
  73 + std::cout << message << std::endl;
  74 +
  75 + } // endcatch: catched any application error
  76 +
  77 + // Return
  78 + return rc;
  79 +}
... ...
src/ctobssim/ctobssim.cpp
  View file @b930a25
... ... @@ -621,7 +621,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
621 621 GCTAEventList* events = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
622 622  
623 623 // Extract ROI
624   - GSkyDir dir = events->roi().centre().skydir();
  624 + GSkyDir dir = events->roi().centre().dir();
625 625 double rad = events->roi().radius();
626 626  
627 627 // Dump simulation cone information
... ...
src/ctools-setup.in
  View file @b930a25
... ... @@ -29,7 +29,7 @@ fi
29 29 prefix="@prefix@"
30 30 exec_prefix="@exec_prefix@"
31 31 ctools_bin="@bindir@"
32   -ctools_lib="@libdir@@TEST_ENV_DIR@"
  32 +ctools_lib="@libdir@"
33 33 uname_system=`uname`
34 34 python_dir="@pythondir@"
35 35 python_exec_dir="@pyexecdir@"
... ...
src/ctskymap/ctskymap.cpp
  View file @b930a25
... ... @@ -492,7 +492,7 @@ void ctskymap::map_events(GCTAObservation* obs)
492 492  
493 493 // Determine sky pixel
494 494 GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
495   - GSkyDir dir = inst->skydir();
  495 + GSkyDir dir = inst->dir();
496 496 GSkyPixel pixel = m_skymap.dir2xy(dir);
497 497  
498 498 // Skip if pixel is out of range
... ...
test/test_ctools.sh
  View file @b930a25
... ... @@ -34,11 +34,12 @@ $ECHO &quot;***************&quot;
34 34 #
35 35 # Define executables
36 36 # ==================
37   -ctobssim=../src/ctobssim/ctobssim
38   -ctskymap=../src/ctskymap/ctskymap
39 37 ctbin=../src/ctbin/ctbin
40   -ctselect=../src/ctselect/ctselect
41 38 ctlike=../src/ctlike/ctlike
  39 +ctmodel=../src/ctmodel/ctmodel
  40 +ctobssim=../src/ctobssim/ctobssim
  41 +ctselect=../src/ctselect/ctselect
  42 +ctskymap=../src/ctskymap/ctskymap
42 43  
43 44  
44 45 #
... ... @@ -136,6 +137,57 @@ $ECHO &quot; ok&quot;
136 137  
137 138  
138 139 #
  140 +# Test ctmodel
  141 +# ============
  142 +$ECHO -n "Test ctmodel: "
  143 +#
  144 +# Run 1
  145 +$ctmodel infile="cntmap.fits" \
  146 + outfile="modmap1.fits" \
  147 + caldb="irf" \
  148 + irf="kb_E_50h_v3" \
  149 + srcmdl="data/crab.xml"
  150 +$ECHO -n "."
  151 +if [ -s "modmap1.fits" ]
  152 +then
  153 + $ECHO -n "."
  154 +else
  155 + $ECHO " modmap1.fits file is not valid"
  156 + exit 1
  157 +fi
  158 +#
  159 +# Run 2
  160 +$ctmodel infile="NONE" \
  161 + outfile="modmap2.fits" \
  162 + caldb="irf" \
  163 + irf="kb_E_50h_v3" \
  164 + srcmdl="data/crab.xml" \
  165 + ra=83.63 \
  166 + dec=22.01 \
  167 + tmin=0.0 \
  168 + tmax=1800.0 \
  169 + emin=0.1 \
  170 + emax=100.0 \
  171 + enumbins=20 \
  172 + nxpix=200 \
  173 + nypix=200 \
  174 + binsz=0.02 \
  175 + coordsys="CEL" \
  176 + xref=83.63 \
  177 + yref=22.01 \
  178 + proj="CAR"
  179 +$ECHO -n "."
  180 +if [ -s "modmap2.fits" ]
  181 +then
  182 + $ECHO -n "."
  183 +else
  184 + $ECHO " modmap2.fits file is not valid"
  185 + exit 1
  186 +fi
  187 +$ECHO " ok"
  188 +
  189 +
  190 +#
139 191 # Test ctselect
140 192 # =============
141 193 $ECHO -n "Test ctselect: "
... ...